§ ABSTRACT

Fracture mechanics is crucial for many fields of engineering applications, as precisely predicting failure of structures and parts is required for efficient designs.
The simulation of failure processes is, both from a mechanical and a numerical point of view, challenging, especially for inhomogeneous materials, where the microstructure influences crack initiation and propagation and might lead to very complex crack patterns.

The phase field method for fracture is a promising approach to encounter such materials, since it is able to describe complex fracture phenomena like crack kinking, branching and coalescence.
Moreover, it is a largely mesh independent approach, given that the mesh is homogenous in the area of the crack.
However, the broadly used formulation of the phase field method is limited to isotropic materials and does not account for preferable fracture planes defined through the material's microstructure.

In this work, the method is expanded to take orthotropic constitutive behavior and preferable directions of crack propagation into account.

We show that by using a stress-based split and multiple phase field variables with preferable fracture planes, in combination with a hybrid phase field approach, a general framework can be found for simulating anisotropic, inhomogeneous materials.

The stress-based split is based on fictitious crack faces and is, herein, expanded to support anisotropic materials.
Furthermore, a novel hybrid approach is used, where the degradation of the sound material is performed based on a smooth traction free crack boundary condition, which proves to be the main driving factor for recovering commonly observed crack patterns.

This is shown by means of a detailed analysis of two examples: a wooden single edge notched plate and a wood board with a single knot and complex fiber directions.
In both cases, the proposed novel hybrid phase field approach is able to realistically reproduce complex failure modes.

§ INTRODUCTION

The use of fracture mechanical modeling approaches is crucial in areas of applied engineering.
Being able to describe failure processes of structures and components realistically and allows for their optimization with regard to design and material consumption, while ensuring high reliability standards.
Simulation of failure processes has been the topic of numerous publications over the recent decades and is, due the complexity of possible material failure mechanisms, still an area that is being researched very intensively.

The foundation of most studies on brittle fracture processes is the work of Griffith1921, which describes fracture in terms of a critical fracture energy release rate required for crack propagation.
Based on Griffith's work, Irwin1958 introduced the so-called stress intensity factors to characterize stress fields around the crack tip, depending on geometry and load.
While those linear elastic fracture mechanics theories can be applied to describe crack propagation of existing cracks, effects like crack kinking, branching, coalescence, initiation and cohesive behavior are not covered.
With the broad establishment of finite element methods for problems related to continuum mechanics, also new methods for simulating fracture processes emerged:
Approaches based on remeshing and usage of special crack tip elements barsoum1977,shahani2009, the node split method peng2012, cohesive elements and cohesive zones models Barenblatt1962,ponnusami2015 and XFEM moes2002.

Those methods allow overcoming some previously mentioned limitations of theories rooted in Griffiths's work, however, each approach comes with its own weaknesses.

One of the most recent and promising method is the so-called phase field method for fracture.
This method was initially proposed by Francfort1998 and is also rooted in Griffith's theory of brittle fracture, however, formulated by a variational approach, through which the total energy of the system is minimized.
The main advantage is that no predefined crack paths are needed and branching as well as coalescence of cracks is naturally included in this approach.
However, finding a solution to the proposed framework turned out to be very difficult.
Therefore, in Bourdin2000,Bourdin2008 a regularization method was developed that allows the minimization problem to be solved numerically efficient.
By introduction of an auxiliary field $d(x) \in [0,1]$ – the so-called crack phase field – the crack discontinuity is modeled by including a smooth transition zone from intact ($d = 0$) to cracked ($d = 1$) solid.
The width of this transition zone is controlled by a regularization or length scale parameter.
As this parameter approaches zero (recovering the discontinuous transition from solid to crack), the solution gamma-converges to Griffith's theory.

By expressing cracks in form of a field variable, the mentioned complex fracture phenomena like kinking, branching and coalescence, naturally arise from the defining system of differential equations.
Thus, the phase field method theoretically allows crack topologies of arbitrary complexity, only limited by the mesh size and mesh structure.
This motivates usage of the phase field method for materials with a complex micro- or macrostructure, like concrete, fiber-reinforced composites, polycrystalline structures and wood.
In those materials, the micro- and macrostructure strongly affects both the elastic behavior, in the sense of having anisotropic constitutive relationships, and the crack topology, by introduction of favorable fracture planes due to weak principal material directions.

§.§ Fracture phenomena of wood

For wood, this results in crack topologies driven by both the direction of least resistance, orthogonal to the wood fiber direction, and the maximum principal stress Smith2003.
This causes the often observed zig-zag pattern, where cracks jump from one growth layer to another, representing a combination of Mode-I, Mode-II and Mode-III failure modes.

In this way, materialspecific microstructural features might influence the fracture behavior and crack propagation during and after crack formation in complex materials.
In wood, crack growth is mainly triggered by defects in the cell wall material at a microscopic level.
This induces the often observed decrease in macroscopic stiffness, evidently visible in load-deflection plots, the nonlinear behavior before reaching the peak load.
Close to the actual peak load the microscopic cracks localize and the actual macroscopic crack and fracture process zone forms.
Failure processes in wood after crack initiation, like cohesive behavior, where studied by vasic2002 and are also observable in the experimental study of Dourado2008.
Their research concluded with identifying so-called fiber bridging as the main cause of toughening effects at the crack tip.

§.§ The phase field method for anisotropic materials

The phase field models proposed by Amor2009, Miehe2010b contain the assumptions of isotropic constitutive behavior and ideal-brittle fracture.
However, in recent years, three approaches to allow for consideration of anisotropic behavior were published:
BLEYER2018213 proposed a method that introduces additional phase field DOF, which are uncoupled in the geometrical terms of the phase field equation and thus allow for different fracture energy release rates and length scale parameters.
The actual coupling is introduced on a constitutive level.
TEICHTMEISTER20171 presented a different approach for anisotropic fracture by adding preferable directions through a second order tensor – the so-called structural tensor.
This tensor scales the gradient of the crack phase field and imposes an orientation on the geometrical terms of the phase field equation.
A similar method was pursued by NGUYEN2017279.
Here, a model with additional phase field DOF is used where each of the phase variables is linked to its own second order structural tensor that invalidates crack growth in certain direction.
Effectively, by suppressing one direction, favorable fracture planes, defined by the microstructure of the material, can be considered.
The tensor used in their work is similar to the transversely isotropic formulation found in the work of TEICHTMEISTER20171.

§.§ Cohesive behavior in the phase field method

Already the original publication Bourdin2000, covering the numerical implementation of the model from Francfort1998, discusses the role of the length scale parameter.
This was further elaborated by Amor2009.
Originally, the length scale parameter was seen as a purely numerical value in the context of the regularization scheme.
However, due to the similarities of the phase field method with gradient damage approaches, considering it a material parameter is not far-fetched.
As already outlined, with this parameter tending towards zero, the solution gamma-converges to Griffith's theory, describing ideal brittle failure.
By increasing this parameter, more complex softening behavior can be described, and ductile effects of quasi-brittle material failure can be implemented
An alternative approach is proposed by Wu2018a in form of the so-called non-standard phase field method, where the phase field's characteristic functions are tuned to a cohesive zone model with a predefined traction separation constitutive relation.
In their work, the length scale parameter is still regarded as a purely numerical value.
However, including the material's tensile strength into the formulation allows controlling the ductile behavior.

§.§ Scope of this work

In this paper, a cohesive phase field implementation is presented based on the above-mentioned non-standard model from Wu2018a.
Furthermore, the outlined methods for considering anisotropic behavior are evaluated and their possible combinations and limitations are discussed.
Based on the work of hu2020PhasefieldModelFracture, a novel hybrid approach Ambati2014, employing a crack boundary condition for the degradation part of the phase field equation, is formulated.

The paper is organized as follows:
<Ref> provides the theoretical background and introduces methods and concepts from recent literature.
Additionally, a generalization of the stress split from steinke2018 for anisotropic materials is given and the hybrid approach, utilizing a crack boundary condition, is described.
In <Ref> the introduced phase field models are evaluated using two numerical examples with different levels of complexity.
First, a simple single edge notch plate with a predefined material direction is looked at.
The simulation results serve as a basis for the discussion of the hybrid approach described in this work.
Secondly, a more complex realistic model of a wooden board with a single knot is simulated to show the effects of spatially varying principal material directions, resulting from different wood fiber orientations.
The paper closes with a short summary of the introduced concepts for applying the phase field model to orthotropic materials showing different complex damage mechanisms, and gives a brief outlook and concluding remarks.

§ FUNDAMENTALS AND METHODS

§.§ Fundamentals

(a) sharp and (b) diffuse representation of a crack topology.
Within the diffuse representation, the sharp crack $\Gamma$ in the body $\Omega$ is approximated using multiple diffuse crack fields, which depend on the crack phase fields $\bm{d}$.

To verify the various methods able to take anisotropic and quasi-brittle behavior into account, the unified phase field model from Wu2018 is generalized to include $n$ different phase field variables BLEYER2018213 and a second order structural tensor $\bm{A}_i$ TEICHTMEISTER20171,NGUYEN2017279, where $i$ denotes the $i$-th phase field.
Excluding body forces and surface tractions, the regularized form of the total energy $\Pi$ of the system, defined on the domain $\Omega$, thus reads:

\begin{align}
\Pi(\bm{u},\bm{d}) = \int_\Omega\left[ \psi^+(\bm{u},\bm{d}) + \psi^-(\bm{u}) \right]d\Omega + \sum_i^n G_{c,i} \int_\Omega \gamma_i(d_i)d\Omega,
\label{eq:phasefieldeq}
\end{align}

where $\bm{u}$ is the displacement field, $\bm{d}$ the phase field of dimensionality $n$, $G_{c,i}$ the critical energy release read of phase field $i$ and $\gamma_i$ is the regularized crack surface density functional that approximates the sharp crack surface.
As the driving force for crack propagation is energy based, a criterion is needed for preventing fracture under pure compressive stress states modes.
Thus, the strain energy density is separated into $\psi^+(\bm{u},\bm{d})$ which contributes to fracture and $\psi^-(\bm{u})$ which does not.
Only $\psi^+$ depends on the crack phase field $\bm{d}$, in such a way that the so-called degradation function $\omega(d_i)$, which damages the solid $\Omega$, degrades the strain energy density.

Within the unified phase field theory, the crack surface density functional is defined as

\begin{align}
\gamma_i(d_i) = \frac{1}{c_{0,i}}\left[ \frac{1}{l_i} \alpha(d_i) + l_i\nabla{}d_i\bm{A}_i\nabla{}d_i \right]\;\text{with}\; c_{0,i} = 4\int_0^1\sqrt{\alpha(\zeta)}d\zeta.
\end{align}

As outlined in Kuhn2015, Pham2010, the function $\alpha(d_i)$ defines the local part of the dissipated fracture energy and determines the ultimate crack phase field.
It satisfies the properties $\alpha(d=0) = 0$ and $\alpha(d=1) = 1$ Wu2018.
In order to recover the actual surface measure of the crack set for $l_i \to 0$, the normalization constant $c_{0,i}$ is needed.
$\bm{A}_i$ is the so-called structural tensor, which scales the gradient of the crack phase field to define preferable or invalid crack propagation directions.
In a multi-phase field setting, the definition

\begin{align}
\bm{A}_i = \bm{I} + \beta_i\left(\bm{I} - \bm{a}_i\otimes\bm{a}_i \right)\label{eq:structural_tensor}
\end{align}

from NGUYEN2017279 can be used, which allows assigning material directions $\bm{a}_i$ and penalty factors $\beta_i$ for penalizing planes not orthogonal to the material directions for each phase field variable $d_i$.
For $\beta_i = 0$, the standard isotropic formulation of the crack surface density functional is recovered.

The two functions that mainly influence the fracture process are the degradation function $\omega(d_i)$ and $\alpha(d_i)$.
As described by BLEYER2018213, while multiple phase field variables are geometrically uncoupled ($\alpha(d_i)$ is defined for each DOF), they are coupled in the constitutive relation ($\psi^+(\bm{u},\bm{d})$ is defined for all DOF).

A general expression for those two functions is given by Wu2017, in the following way

\begin{align}
\alpha_i(d_i) &= \xi d_i +(1-\xi) d_i^2 \quad \forall d_i \in [0,1]\;\xi \in [0,2] \;\text{and} \label{eq:localpart}\\
\omega_i(d_i) &= \frac{(1-d_i)^p}{(1-d_i)^p+Q(d_i)} \quad p \geq 2 \label{eq:degrad} \\
Q_i(d_i) &= a_{1,i}d_i + a_{1,i}a_{2,i}d_i^2 + a_{1,i}a_{2,i}a_{3,i}d_i^3, \label{eq:Q_a1_a2_a3}
\end{align}

where $a_{1,i}$, $a_{2,i}$ and $a_{3,i}$ are coefficients that can be calibrated to model a certain cohesive behavior related to the $i$-th phase field.

The general expression for the local part of the dissipated fracture energy resembles the well-known monotonous model.
Alternative formulations are of type of a double well function.
Such a formulation has an interesting property, because it provides an energy barrier between the undamaged and the fully damaged state, and thus, naturally models irreversibility of the crack phase field Kuhn2015.
However, those functions are barely used because of the two states being energetically equivalent, which can lead to phase fields expanding perpendicular to the crack path.
<Ref> can be specialized for the two commonly used models $\alpha(d_i) = d_i$ for $\xi = 1$ (AT-1) and $\alpha(d_i) = d_i^2$ for $\xi = 0$ (AT-2) ambrosio1990.

As outlined in Miehe2010b, the degradation function must satisfy $\omega(d_i = 0)=1$ and $\omega(d_i = 1)=0$.
Furthermore, as the first derivative of the degradation function with respect to the phase field variable controls the amount of energetic driving force, $\omega'(d_i = 1)=0$ is needed, in order to eliminate this elastic driving term once full damage is reached Miehe2010b,Kuhn2015.
This ultimately stops further crack growth in regions characterized by $d_i = 1$.
steinke2018 discuss the additional soft requirement of $\omega(d_i = 0) \neq 0$, which, if not satisfied, as well leads to the elimination of the elastic driving term.
For $d_i = 0$, the undamaged state, crack growth is hindered, and no phase field evolution can take place.
Therefore, such a model requires additional treatment in form of a numerical perturbation of the initial state, such that the energetic driving forces become unequal to zero.
One commonly used definition is $\omega(d_i) = (1-d_i)^2$ Miehe2010.
The general expression in <Ref> contains this simple case for $p=2$, $a_1=2$, $a_2=-\sfrac{1}{2}$ and $a_3=0$.

The remaining part left to be specified from <Ref> is the strain energy density split.
In order to properly discuss the various kinds of methods to approach this separation, at first, two commonly applied methods for dealing with coupled equations within the variation framework are discussed.

§.§ Isotropic, anisotropic and hybrid formulation

As described by Ambati2014, there are two basic formulations originating from the regularized variational framework Bourdin2000,Bourdin2008, the isotropic formulation and the anisotropic formulation Miehe2010b, Miehe2010[The terms isotropic and anisotropic are not to be interpreted in terms of the local material behavior, but they refer to the decomposition of the strain energy density].
The isotropic formulation does not contain the additive decomposition of the strain energy density and, thus, gives a linear relation in $\bm{u}$, reading

\begin{align}
\bm{\sigma}(\bm{u},\bm{d}) = \frac{\partial\psi(\bm{u},\bm{d})}{\partial\bm{\varepsilon}}.
\label{eq:constiso}
\end{align}

The anisotropic formulation

\begin{align}
\bm{\sigma}(\bm{u},\bm{d}) = \frac{\partial \psi^+(\bm{u},\bm{d}) + \psi^-(\bm{u}) }{\partial\bm{\varepsilon}},
\label{eq:constaniso}
\end{align}

however, contains the split and is thus non-linear in $\bm{u}$.
This property simplifies the solution process when the so-called staggered approach is used (see <Ref>), as the deformation subproblem can be treated as an uncoupled linear problem.

One obvious downside of the isotropic formulation is that every deformation state is degraded equally, thus unphysical behavior like interpenetration of crack faces or crack growth under pure compressive stress states can occur.
Each effect is related to a different aspect of the evolution equations.
Crack growth is rooted in the energetic driving force and interpenetration of crack faces in the constitutive relation.
The so-called hybrid formulation Ambati2014 combines the advantages of both formulations, the linear behavior of the isotropic formulation and the physical, more appropriate, modeling of the anisotropic one.
This is achieved by using the constitutive relation from <Ref> with the additional constraint that for all $\bm{u}$ with $\psi^+$ smaller than $\psi^-$ (the passive energy parts outweighs the crack driving one), the material is treated as undamaged, $d_i = 0$.
For the energetic driving force, the anisotropic formulation is used, leading to the proper crack propagation behavior.
As shown in Ambati2014, the hybrid formulation manages to produce results qualitatively and quantitatively similar to the anisotropic formulation, with a computation effort close to the isotropic formulation.

§.§ Fracture contributing and passive parts

The energy split into a part that drives and is affected by fracture and into a part that is neither affected by nor drives fracture is one key ingredient of the phase field formulation for ensuring a realistic fracture behavior.
As outlined in <Ref>, by ignoring this split, using the isotropic formulation, effects like interpenetration of crack faces or crack growth under pure compressive stress states can occur.
This unphysical behavior can be eliminated by properly splitting the energy, such that the energetic driving force and also the damaged part in the constitutive relation is only related to the contributing part $\psi^+$.

Besides this physical motivation, there is also a simple conceptional one: The additive decomposition $\psi = \psi^+ + \psi^-$, must retain its validity.
In dijk2020, this is very well explained based on the two common methods for splitting the strain energy density: the spectral decomposition by Miehe2010b and the volumetric-deviatoric decomposition by Amor2009.
Starting with the well-known formulation of the strain energy density, the fundamental idea is to separate stresses and strains into strictly fracture-contributing and passive parts.
Hence, the strain energy density reads

\begin{align}
\psi &= \frac{1}{2} (\bm{\sigma}^+ + \bm{\sigma}^-) : (\bm{\varepsilon}^+ + \bm{\varepsilon}^-)\label{eq:standard_decompo}\\
&= \frac{1}{2} ( \underbrace{\bm{\sigma}^+ : \bm{\varepsilon}^+}_{\psi^+} + \bm{\sigma}^+ : \bm{\varepsilon}^- + \bm{\sigma}^- : \bm{\varepsilon}^+ + \underbrace{\bm{\sigma}^- : \bm{\varepsilon}^-}_{\psi^-}).
\end{align}

So, for the additive decomposition to hold, the terms $\bm{\sigma}^+ : \bm{\varepsilon}^-$ and $\bm{\sigma}^- : \bm{\varepsilon}^+$, consisting of contributing and passive parts, must vanish.

The spectral decomposition is conceptually rooted in fracture mechanics, so tensile principal stresses (and strains) result in crack opening, while compressive principal stresses (and strains) do not.
This idea motivates the decomposition to be performed in the principal strain space.
There, for an isotropic material, where the directions of the principal strains and principal stresses coincide, also the contributing strains are orthogonal to the passive stresses and the passive strains are orthogonal to the contributing stresses.
Therefore, as required, the terms containing both contributing and passive parts indeed vanish.
By the same reasoning, the volumetric-deviatoric decomposition is also valid for an isotropic material, as the volumetric and deviatoric parts are orthogonal.
So both models can be applied to isotropic materials.
However, for non-isotropic materials, this orthogonality between principle stresses and strains no longer exists and, thus, neither of the above models can be used, as $\psi \neq \psi^+ + \psi^-$.

The phase field formulation consists of two essential parts, one being the formulation of the energetic driving force and the other being the actual constitutive behavior.
On those two parts, the above described decomposition has a very different impact.
Regarding the first part, Miehe2015 introduced the concept of a dimensionless crack driving function, which is a generalization of the history function Miehe2010b and allows replacing the energetic driving force by an arbitrary failure function, like a maximum stress, a maximum strain or a Tsai-Wu GULTEKIN201823 criterion.
While this opens the phase field formulation to a variety of materials, the restrictions on the constitutive relation remain.
However, the current literature finally provides three approaches of a valid strain energy decomposition for non-isotropic materials:

* A proper split for anisotropic materials dijk2020, based on generalizations of decomposition from Miehe2010b and Amor2009,  which is not yet validated,

* a stress-based split steinke2018, hu2020PhasefieldModelFracture that is based on a different formulation of <Ref> and

* the hybrid approach described in <Ref>.

The hybrid approach is a simple method to circumvent the validity of the split, as for the constitutive part no split is required.
In a sense, the hybrid approach represents a generalized framework for applying arbitrary crack driving functions, as long as an additional constraint ensures that crack faces cannot interpenetrate.
The additional constraint formulation from Ambati2014, where for $\psi^- > \psi^+$ contact of the crack faces is assumed, however, might not be suitable for any crack driving function, thus, requiring an alternative formulation of this constraint, suited for ensuring a physical, sound fracture behavior.

A stress-based split represents an alternative approach for separating the strain energy density into a contributing and a passive part.
The method, as described by steinke2018, is based on the idea that only stresses are additively decomposed.
So <Ref> changes to

\begin{align}
\psi &= \frac{1}{2} (\bm{\sigma}^+ + \bm{\sigma}^-) : (\bm{\varepsilon}) = \frac{1}{2} ( \underbrace{\bm{\sigma}^+ : \bm{\varepsilon}}_{\psi^+} + \underbrace{\bm{\sigma}^- : \bm{\varepsilon}}_{\psi^-}).
\end{align}

Obviously, this approach is applicable to any constitutive behavior, as there are no terms consisting of both contributing and passive parts.
Nevertheless, steinke2018 specialized the model for isotropic materials, as some observations on the inherent properties of the deformations in presence of a crack, discussed below,  require changes in the initial formulation of $\bm{\sigma}^+$.

The stress based split is performed in a local crack coordinate system, as shown in <Ref>.
This coordinate system related to a fictitious crack surface allows identification of crack driving forces for Mode-I, Mode-II and Mode-III fracture.
By representing the stress tensor in the crack coordinate system, the crack driving and passive stresses can be identified:

\begin{align}
\bm{\sigma}^+ &= \underbrace{\left<\sigma_{rr}\right>_{+}(\bm{r}\otimes\bm{r})}_{\text{Mode-I}} + \underbrace{\bm{\sigma}_{rt} + \bm{\sigma}_{tr}}_{\text{Mode-II}} + \underbrace{\bm{\sigma}_{rs} + \bm{\sigma}_{sr}}_{\text{Mode-III}}\label{eq:stress_modes_1}\\
\bm{\sigma}^- &= \left<\sigma_{rr}\right>_{-}(\bm{r}\otimes\bm{r}) + \bm{\sigma}_{tt} + \bm{\sigma}_{ss} + \bm{\sigma}_{ts} + \bm{\sigma}_{st},\label{eq:stress_modes_2}
\end{align}

where $\bm{\sigma}_{ij} = (\bm{\sigma}:(\bm{i}\otimes\bm{j}))(\bm{i}\otimes\bm{j})$, the contributions to the stress tensor related to the $i$ and $j$ direction in the crack coordinate system and $\left<\bullet\right>_+ = \sfrac{(\bullet + |\bullet|)}{2}$ and $\left<\bullet\right>_- = \sfrac{(\bullet - |\bullet|)}{2}$.
At this point, the decomposition is still applicable to any constitutive law, however, it can lead to physically inconsistent results, as in the fully damaged state the essential crack boundary conditions – no positive normal stress perpendicular to the crack and no shear stresses along a frictionless crack surface strobl2015NovelTreatmentCrack – are not recovered for certain strain states.

Crack coordinate system of a fictitious crack surface, required to identify Mode-I, Mode-II and Mode-III crack driving forces.

This unphysical behavior can easily be shown by picturing a state of pure crack normal strain ($\varepsilon_{rr} > 0$ and all other strain components equal to zero).
As a generalization of the isotropic model from steinke2018, following a linear elastic constitutive law, such a strain state would result in non-zero stress components:

\begin{align}
\sigma_{rr} &= C_{rrrr}\varepsilon_{rr}\text{ and }\label{eq:normal_strain_1}\\
\sigma_{ij} &= C_{ijrr}\varepsilon_{rr}.\label{eq:normal_strain_2}
\end{align}

A fully developed crack state, however, should be stress-free due to this very strain state, as the two crack surfaces should be able to move freely along the crack's normal vector.
Additionally, the stresses related to Poisson's effect (<Ref>), must be considered in the fracture contributing stresses $\bm{\sigma}^+$, by expressing them in terms of the crack normal stresses from <Ref> as

\begin{equation}
\sigma_{ij} = \frac{C_{ijrr}}{C_{rrrr}}\sigma_{rr},
\end{equation}

which leads to the following formulation of the crack contributing and passive stresses:

\begin{align}
\begin{split}
\bm{\sigma}^+ &= \left<\sigma_{rr}\right>_{+}(\bm{r}\otimes\bm{r}) + \bm{\sigma}_{rt} + \bm{\sigma}_{tr} + \bm{\sigma}_{rs} + \bm{\sigma}_{sr} + \frac{\left<\sigma_{rr}\right>_{+}}{C_{rrrr}} \big[
\end{split}\label{eq:total_active_stress}\\
\begin{split}
\bm{\sigma}^- &= \left<\sigma_{rr}\right>_{-}(\bm{r}\otimes\bm{r}) + \bm{\sigma}_{tt} + \bm{\sigma}_{ss} + \bm{\sigma}_{ts} + \bm{\sigma}_{st} - \frac{\left<\sigma_{rr}\right>_{+}}{C_{rrrr}} \big[
\end{split}\label{eq:total_passive_stress}
\end{align}

<Ref> is based on <Ref>, which contains crack driving stresses identified using classic fracture mechanics failure modes.
In addition, <Ref> also considers the stresses related to Poisson's effect, which must vanish for a fully developed phase field.
Having this formulation for $\bm{\sigma}^+$ allows finding $\bm{\sigma}^-$ from $\bm{\sigma}=\bm{\sigma}^+ + \bm{\sigma}^-$.
It can be shown that this generalized formulation of the stress split contains the isotropic formulation proposed by steinke2018.
In the isotropic case, only the entries $C_{rrrr} = \lambda + 2 \mu$, $C_{ssrr} = C_{ttrr} = \lambda$, required for the additional decomposed stresses, are non-zero, leading to the expression

\begin{equation}
\bm{\sigma}^+ = \left<\sigma_{rr}\right>_{+}(\bm{r}\otimes\bm{r}) + \bm{\sigma}_{rt} + \bm{\sigma}_{rs} + \bm{\sigma}_{tr} + \bm{\sigma}_{ts} + \frac{\lambda}{\lambda + 2 \mu}\left<\sigma_{rr}\right>_{+}(\bm{s}\otimes\bm{s} + \bm{t}\otimes\bm{t}),
\end{equation}

where $\lambda$ and $\mu$ are the two Lamé constants.
This expression matches the one from steinke2018.

§.§.§ A novel hybrid approach based on a crack boundary condition

Referring to the key requirement defined in strobl2015NovelTreatmentCrack that in a fully damaged state, tensile crack normal stresses and shear stresses along a frictionless crack surface should be zero, hu2020PhasefieldModelFracture developed a stress-based decomposition approach based on a smooth traction-free crack boundary condition.
This approach is similar to the stress-based split by steinke2018, however, instead of considering degradation from the perspective of crack driving stresses in Mode-I, Mode-II and Mode-III, they view degradation as a contact problem.
With $\bm{r}$ as the crack face's normal vector, this results in the following decomposition of the stress tensor:

\begin{align}
\bm{\sigma}^+ &= \underbrace{\left<\sigma_{rr}\right>_{+}(\bm{r}\otimes\bm{r})}_{\text{tensile normal stress}} + \underbrace{\bm{\sigma} - \sigma_{rr}(\bm{r}\otimes\bm{r})}_{\text{tangential stress}}\\
\bm{\sigma}^- &= \underbrace{\left<\sigma_{rr}\right>_{-}(\bm{r}\otimes\bm{r})}_\text{compressive normal stress}
\end{align}

The main difference compared to <Ref> is in the treatment of tangential components, which in the case of the crack boundary condition always result in an energetic driving force.
This results in an unrealistic overestimation of the driving strain energy density ($\bm{\sigma}_{tt}$ and $\bm{\sigma}_{ss}$ are considered crack driving).
Therefore, hu2020PhasefieldModelFracture highlight that this decomposition approach should not be used for crack initiation and crack propagation, but should only serve as a boundary condition, which is activated after the phase field variable reaches a certain threshold.

An alternative approach for dealing with the problem of having unphysical crack driving forces is utilizing the properties of the hybrid approach (<Ref>).
As discussed in <Ref>, this method allows arbitrary combinations of energetic driving forces and definitions of the constitutive behavior, given that the solver is based on the staggered approach (see <Ref>).
Commonly, the isotropic formulation from <Ref> is used for defining the constitutive behavior, which, however, requires an additional constraint for preventing interpenetration of crack faces.
In order to circumvent this additional constraint, we therefore propose using

\begin{align}
\bm{\sigma}^+ &= \left<\sigma_{rr}\right>_{+}(\bm{r}\otimes\bm{r}) + \bm{\sigma}_{rt} + \bm{\sigma}_{tr} + \bm{\sigma}_{rs} + \bm{\sigma}_{sr} + \bm{\sigma}_{tt} + \bm{\sigma}_{ss} + \bm{\sigma}_{ts} + \bm{\sigma}_{st} \;\text{and}\label{eq:crack_boundary_active}\\
\bm{\sigma}^- &= \left<\sigma_{rr}\right>_{-}(\bm{r}\otimes\bm{r})\label{eq:crack_boundary_passive},
\end{align}

which matches the formulation from hu2020PhasefieldModelFracture, for defining the constitutive behavior, where only $\bm{\sigma}^+ $ is degraded.
The crack driving part of the coupled system remains based on the orthotropic stress split derived in <Ref>.
This additionally bypasses the need for considering the crack boundary condition only after the phase field variable reaches a certain threshold, as the driving parts are now rooted in the classic fracture mechanic failure modes.
To conclude, we propose using a stress-based decomposition, where crack driving stresses are identified by fracture mechanics failure modes (<Ref>), in a hybrid-approach, combined with a degradation function, where stresses are degraded, such that the resulting constitutive behavior matches the one of a traction-free crack surface (<Ref>).
This allows physical, sound estimation of crack driving forces for orthotropic materials and proper modeling of crack faces.

§.§.§ Application of the stress split in a multi-phase field theory

In highly orthotropic materials, fracture is driven based on two principles Smith2003:

* Cracks following the direction of least resistance, defined by the microstructure of the material.

* Cracks opening perpendicular to the largest principal stress, thus, leading to a maximum reduction of the total energy.

For wood this results in the often observed zig-zag pattern (<Ref>), where cracks follow the path of maximum total energy reduction, until reaching a growth ring.
At the growth ring, which is essentially a weak interface, the crack direction changes to the direction of least resistance, which for wood is always parallel to its fiber direction.
So cracks are likely to follow the material's structure, for wood the longitudinal, radial or tangential direction.

Typical crack pattern observed for wood.
A crack switches from following the direction of least resistance along the fiber, to the direction of maximum total energy reduction perpendicular to the largest principal stress.
It changes its orientation again, when it hits the next growth ring. Smith2003

The fundamental part of the stress-based decomposition is the definition of the crack coordinate system.
Due to the material characteristics and often observed crack pattern of wood, instead of considering an arbitrary crack face orientation, the identification of three crack coordinate systems defined by $\bm{r}_i$, $\bm{s}_i$ and $\bm{t}_i$ (see <Ref>) is plausible, where $\bm{r}_1$ is the longitudinal, $\bm{r}_2$ the radial and $\bm{r}_3$ the tangential direction.
This results in three different crack driving energy terms, where according to the principle of maximum dissipation, the failure mode with the highest energy release determines the main cause of failure.

In order to consider the very different fracture toughness and strength of each of the possible types of cracks (longitudinal, radial, tangential), each material axis is related to only a single phase field variable.
This is similar to approaches taken for crystalline materials for considering the influence of cleavage planes NGUYEN2017279.
With this multi-phase field description, only the phase field variable related to driving failure mode is activated for degradation.
In comparison to other multi-phase field models, BLEYER2018213, NGUYEN2017279, this approach results in no coupling of the phase field variables in the constitutive relation, as a single phase field variable is already sufficient to describe the state of a fully developed crack.
Thus, the strain energy terms $\psi^+(\bm{u},\bm{d})$ and $\psi^-(\bm{u},\bm{d})$ in <Ref>, are replaced by

\begin{align}
\psi_i^+(\bm{u},\bm{d}) = \omega_i(d_i) \psi^+_i(\bm{u})\text{ and } \psi_i^-(\bm{u}) = \psi_i(\bm{u}) - \psi^+_i(\bm{u}),
\end{align}

respectively, where $i$ is the index of the defining failure mode.

From a numerical point of view, having no coupling between the phase field DOF is favorable.
However, choosing one driving failure mechanism introduces a strong non-linearity, as the strain energy density function depends on $i$, for which the total energy is not a continous function.
This makes solving the problem very difficult, if not impossible.
Therefore, instead of evaluating strong non-linearities (Heaviside functions or the Macaulay brackets) based on the current state variables, they are computed using the deformation and phase field values of the last converged increment.
Given that the increments are sufficiently small, this vastly improves the convergence rate, while leading to similar results.
Our tests showed no significant influence on the obtained solutions, except that the algorithm is more robust.

§.§ Solver

Solving the posed minimization problem is a challenging task, because of the regularized functionals being non-convex in the state variables Amor2009.
Furthermore, to ensure thermodynamical consistency and consider the fully dissipative nature of crack growth, it is necessary to impose irreversibility constraints Miehe2010b.
Solving a non-convex bound-constrained optimization problem for the global optimum is hard; in this case, due to the requirement of a very fine mesh and the resulting problem size, it is almost impossible.
Therefore, searching for local minima is the only option.
Here, mainly two different approaches are pursued:
Application of an alternate minimization scheme (also referred to as the staggered approach) Bourdin2000,Miehe2010 and the monolithic approach Miehe2010b.
The consideration of the constraints is not part of these two approaches and is carried out separately.

In the monolithic approach, the problem is solved by using a standard Netwon-Raphson procedure on the full set of state variables.
This can lead to convergence issues, due to the non-convexity of the problem in those variables.
Improvements to these algorithms were proposed, by gerasimov2016 in the form of a line search that is performed in case the energy functional increases and more recently by kopanicakova2020 through additional application of a thrust region method.
Still, most of the recent works (as well as this) do not rely on the monolithic approach and instead use the staggered one.

Staggered solution process for a single load increment.
Convergence is assured by imposing a limit on the size of the state variable's increment.
Additionally, the residual of the phase field subproblem must lie below a threshold.
The convergence of the overall problem is given by checking whether the last change in the crack phase field only resulted in a converged state in the displacement field.

The fundamental idea of the staggered approach is to separate the problem into two subproblems, where one is used to determine the displacement field at a constant crack phase field and the other to find the crack phase field at a constant displacement field.
The problems are solved in an alternating manner, such that each is based on the resulting state variable of the other.
One main advantage of the staggered approach is that the subproblems are convex and, thus, the robustness of the solution process is relatively high Ambati2015.
The alternating scheme, however, requires a higher number of iterations compared to the monolithic approach and thus greater computational effort.
Also, convergence criteria are required not only for the overall problem, but for both subproblems.

The algorithmic scheme applied in this work is outlined in the flowchart in <Ref>.
Similar to what was proposed by Amor2009, the following criteria are used as a convergence measure for the subproblems in iteration $n$ and increment $k$,

\begin{align}
\sum_{i \in \text{DOFs}}\frac{|\dFE^k_{n+1,i}-\dFE^k_{n,i}|}{|\dFE^k_{n+1,i}|} \leq \epsilon_d \quad
\sum_{i \in \text{DOFs}}\frac{|\uFE^k_{n+1,i}-\uFE^k_{n,i}|}{|\uFE^k_{n+1,i}|} \leq \epsilon_u \quad
\text{with}\quad
\epsilon_u < \epsilon_d.
\end{align}

Here, $\hat{\bullet}$ refers to the discretized representation of the continuous field $\bullet$.
Additionally, the L${^2}$-norm of the crack phase field's residual $\bm{\hat{R}}^k_{\dFE,n+1}$ must be smaller than or equal to $\epsilon_{R,d}$.
Convergence of the overall problem is assured by requiring convergence of the deformation problem and the phase field problem.
If the set of state variables resulting from the last phase field step ($\dFE^k_{n+1}$ and $\uFE^k_{n}$) result in a converged state for the deformation problem, the newly obtained deformation state $\uFE^k_{n+1}$ is ignored (as the phase field subproblem converged for $\uFE^k_{n}$) and $\dFE^k_{n+1}$ and $\uFE^k_{n}$ are accepted as a solution.
The whole algorithm uses an adaptive time increment stepping scheme, such that a better performance can be achieved by having larger increments in less critical regions (linear elastic regime) and smaller increments in critical ones (close to peak load).

As mentioned above, to ensure thermodynamical consistency, an irreversibility constraint on the phase field variable is required.
Additionally, it is also necessary to ensure the bound constraint $d(x) \in [0,1]$.
Lorenzis2020 give a quite comprehensive overview of current approaches.
There are various methods and considerations that allow formulating those constraints in various places in the implementation.

Regarding the bound constraint, in <Ref>, the requirement on the degradation function, $\omega'(d_i = 1)=0$, was discussed.
Essentially, this disables the energetic driving force for $d_i = 1$, and therefore, a further development of the related phase field $i$.
Thus, the upper bound is naturally satisfied by the choice of such a degradation function.
Given that there is a solution for enforcing irreversibility of the phase field problem in general, the lower bound is implicitly fulfilled as the crack phase field is initially zero.
Thus, for every subsequent increment the condition $0 \leq \dFE^{0} \leq \dFE^{k} \text{ }\forall\text{ } k$ holds.
Conceptually, two different schemes for the inequality constraint are currently used:
Applying constraints on a local level, like the history field Miehe2010 and the penalty method Gerasimov2019, or applying constraints on the global level, by application of Dirichlet-type boundary conditions Bourdin2000 or active set methods Heister2015, Hintermueller2002, which are essentially Dirichlet-type boundary conditions on a subset of $\dFE$.
In this work, the computationally less expensive active set reduced space method yang2016a, which is both qualitatively and quantitatively similar to the primal-dual active set method from Heister2015, is used (for a comparison see <Ref>).

A set is said to be active when the inequality constraint is violated.
In <cit.> a box constraint solver is discussed, hence identifying two sets, a lower bound active set $\mathcal{A}_\phi$ and an upper bound active set $\mathcal{A}_\psi$.
The problem is, subsequently, solved on the inactive set

\begin{align}
\mathcal{I}(\dFE^k_{n+1}) = \mathcal{S} \setminus \left(\mathcal{A}_\phi(\dFE^k_{n+1}) \cup \mathcal{A}_\psi(\dFE^k_{n+1})\right),
\end{align}

where $\mathcal{S}$ is the set of all crack phase field DOF.
The DOF values on the active set are fixed to the boundary values using hard Dirichlet-type boundary conditions.
Applying this method would allow considering alternative degradation functions for which $\omega'(d_i = 1) \neq 0$, however, in this work no functions of this kind are used.
The active sets are computed as follows[In their work, yang2016a compare the current DOF's value with the one from the last converged state by equality. Given that they apply a special operator that cuts off values lower than the lower bound and larger than the upper bound, comparing by lower than or larger than, respectively, leads to the same result. We use those operators instead, to make the comparison to the primal-dual active set algorithm from Heister2015 more clear.]:

\begin{align}
\mathcal{A}_\phi(\dFE^k_{n+1}) &= \left\{  i \in \mathcal{S} | \left(\dFE^k_{n+1}\right)_i \leq \left(\dFE^{k-1}\right)_i \text{ and } (\bm{\hat{R}}^k_{\dFE,n+1})_i > 0 \right\}\label{eq:yang-active-set}\\
\mathcal{A}_\psi(\dFE^k_{n+1}) &= \left\{  i \in \mathcal{S} | \left(\dFE^k_{n+1}\right)_i \geq 1 \text{ and } (\bm{\hat{R}}^k_{\dFE,n+1})_i < 0 \right\},
\end{align}

where the lower bound is the element-wise restriction that the current state variable must be larger than or equal to the state variable from the last converged increment, $k-1$.
These methods also affects the convergence conditions from above, such that the L$^2$-norm of the residual is only computed on the inactive set.
To conclude, we propose using the staggered approach, as it is more robust than the monolithic one, in combination with the active set reduce spaced method, which allows assuring irreversibility without requiring additional terms like penalty functions in the phase field formulation.

The entire code is implemented in Julia bezanson2017JuliaFreshApproach.
For automatically deriving the element stiffness matrices and residual vectors from the energy formulation, the ForwardDiff-Package revels2016ForwardModeAutomaticDifferentiation is used.
Pardiso 6.0 deconinck2016NeedlesLargeScaleGenomic, kourounis2018NextGenerationMultiperiod, verbosio2017EnhancingScalabilitySelected is employed as the sparse linear solver.

§ RESULTS AND DISCUSSION

In the following Section, the stress split described in <Ref> is assessed based on two different models: A simple notched plate (<Ref>) with varying fiber orientation and a more complex example of a wooden board including a knot (<Ref>).
As the envisioned use of the phase field model described in this work is the application to complex three-dimensional geometries, three-dimensional linear tetrahedral elements are used.
Since the hybrid approach strongly alters the phase field formulation, it can be expected to provide different crack topologies, compared to the variationally consistent formulation (the anisotropic formulation from <Ref>).

For all examples, the material stiffness tensor is defined as $\bm{C}_{\text{LLLL}} = 9000.016,\, \bm{C}_{\text{LLRR}} = 269.384,\, \bm{C}_{\text{LLTT}} = 175.104,\, \bm{C}_{\text{RRRR}} = 480.096,\, \bm{C}_{\text{RRTT}} = 118.528,\, \bm{C}_{\text{TTTT}} = 270.6,\, \bm{C}_{\text{RTRT}} = 32,\, \bm{C}_{\text{LRLR}} = 552,\, \bm{C}_{\text{LTLT}} = 552$, all in $\si{MPa}$.
This resembles the elastic properties of so-called clear wood, describing wood areas without defects and knots.
Following <cit.>, to the knot in the wooden board example, a stiffness tensor reduced by a factor of $0.5$, compared to the clear wood stiffness tensor, is assigned.
This reduction takes cracks perpendicular to the grain direction, often observed in knots, into account.
The elastic properties are defined in a local cylindrical coordinate system, as is commonly used for describing wood.
$\bm{a}_1$ defines the longitudinal (L) direction, $\bm{a}_2$ the radial (R) direction and $\bm{a}_3$ the tangential (T) direction.

In order to account for the cohesive behavior of wood, the coefficients $a_{1,i}$, $a_{2,i}$ and $a_{3,i}$ in <Ref> are tuned to match a linear softening law.
Based on the analytical solution of a one-dimensional bar problem, Wu2017 gives the following definitions:

\begin{align}
a_{1,i} = \frac{4}{\pi}\frac{l_{ch,i}}{l_i},\; a_{2,i} = -\frac{1}{2}\text{ and } a_{3,i} = 0,\label{eq:a1_a2_a3_softening}
\end{align}

for $\xi = 2$ and $p = 2$, in <Ref>, respectively.
$l_{ch,i}$ defines Irwin's characteristic length, given as $l_{ch,i} = \sfrac{E_{0,i} G_{c,i}}{f_{t,i}^2}$, for the $i$-th phase field, and $l_{i}$ is the length scale parameter for the $i$-th phase field, which is chosen to be larger than the effective element size $l_{\text{eff}}$ (third root of the average volume of the finite elements in region of the expected phase field crack).

§.§ Single edge notched plate

Geometry of the single edge notched plate.
In plane, the plate is fixed in all directions at the bottom edge and only out of plane across the entire back surface.
The load is applied in form of a prescribed vertical deformation along the upper edge.
The fiber direction $\bm{a}_1$ is changed by setting the fiber angle relative to the horizontal direction.
All measurements are in $\si{mm}$.

The notched plate's geometry is depicted in <Ref>.
In plane, it is supported at the bottom edge and out of plane on the entire back surface.
The load is applied in form of a prescribed vertical deformation along the upper edge.
For considering the orthotropic behavior, the fiber direction ($\bm{a}_1$) is changed by setting the fiber angle relative to the horizontal direction, $0^\circ{}$ meaning $\bm{a}_1$ points into the x-direction and $90^\circ{}$ meaning $\bm{a}_1$ points into the y-direction.
The remaining axes are defined such that the tangential direction ($\bm{a}_3$) always points into the z-direction.
The parameters controlling the phase field problem are given in <Ref>.

Defining parameters for the single edge notch plate problem.

$d_i$     $\beta_i$a      $G_{c,i}$b     $f_{t,i}$c     $\sfrac{l_i}{l_{\text{eff}}}$

$d_1$     5.0     0.05     50.0      4.0

$d_2$     5.0     0.1      14.42     4.0

$d_3$     5.0     0.1       7.21     4.0

a Structural tensor scale in <Ref>

b in $\sfrac{\si{Nmm}}{mm^2}$

c in $\si{MPa}$

In order to reduce the computational effort of such problems, often, the mesh density is increased in regions of a priori known crack paths.
As changing the fiber angle is expected to also change the resulting crack topology, the crack paths cannot precisely be known in advance.
Therefore, all models are consistently meshed with the same effective element size over the entire specimen's geometry.
This greatly reduces the influence of the mesh structure on the resulting crack paths.
Initially, seven different element sizes, ranging from a very coarse mesh with 2018 nodes to a very fine one with 144825 nodes, are tested.
The finest mesh results from a characteristic element size of $0.005$, a value which is also used in other publications, by hu2020PhasefieldModelFracture.

The results of this mesh study are shown in <Ref>.
For both the hybrid and the consistent approach, the total external energy, normalized to the maximum value of the external energy for the specific model and the studied fiber angle, is plotted over the number of nodes.
For all five material directions, with increasing number of nodes, the total external energy shows clear convergence against a value that can already be captured well by the two finest meshes (6 and 7).
This is also reflected in the phase field developments in <Ref> (c), where there is no qualitative difference in the crack topology between mesh 6 and 7, however, a significant change in the failure mode in meshes 1 to 5.
Therefore, only mesh number 7 was used, for the further simulations.
Nevertheless, it should be pointed out that also mesh number 6, though quite coarse compared to common mesh sizes used in literature, could already be used, which would allow a major reduction of the computational effort.
The generally low mesh sensitivity, as also pointed out by Yang2019, is related to the usage of the unified phase field theory from Wu2017, as the regularization parameter is actually considered in the phase field formulation for calibration of the coefficients $a_{1,i}$, $a_{2,i}$ and $a_{3,i}$ in <Ref>, which compensates the size effect resulting from a larger crack phase field.

Results of the mesh study (total external energy) for seven different meshes with increasing mesh density, for 5 different fiber angles and for both the consistent and the hybrid approach.
Additionally, the crack topology for a fiber angle of $72^\circ$ for the hybrid approach is depicted for each of the meshes.

As already pointed out, a main question of this work is to evaluate whether the proposed model is capable of considering effects resulting from the material structure appropriately.
With respect to wood, it is of particular interest, whether the commonly observed zig-zag failure pattern (as shown in <Ref>), arising from cracks that jump between growth layers, can be modeled.
In order to show this effect, various fiber angles ranging from $0^\circ$ to $90^\circ$ are investigated.
<Ref> shows the transition of failure modes, obtained with the hybrid approach, for selected characteristic fiber angles in this range, where the mesh and the load are identical.
In <Ref>, qualitatively similar crack topologies are summarized in a graphic to improve comparability.

The results clearly show that at a certain fiber angle, the failure mode switches from a crack driven by stresses perpendicular to the fiber ($d_2$) to a crack driven by stresses in fiber direction ($d_1$).
The first and foremost observation from those results is that by using the proposed hybrid approach, it is actually possible to recover the zig-zag fracture pattern, even if a completely homogenous mesh and material definition is used.
The main influencing factors are the structural tensor, which forces the geometric phase field evolution to stay on planes perpendicular to the crack normal direction, and considering the driving stresses on the fictitious crack face for each likely crack orientation.
This will further be elaborated at the end of this section, where the influence of the hybrid approach on the appearance of this pattern is discussed.

Comparing Figures <ref> (a) and (c) it becomes obvious that phase field $d_2$ (green) has a stronger tendency to follow the prescribed fiber orientation than phase field $d_1$ (red).
This result suggests that even without explicitly defining a weak interface in-between fibers, cracks perpendicular to the $\bm{a}_2$-direction (<Ref>) are influenced by an effect that has a similar impact as a weak interface.
Nevertheless, the cracks do not exactly follow the fiber direction and the stronger the incline is, the stronger the tendency away from this fracture plane.
Theoretically, by increasing the structural tensor scale, one could increase the affinity to this plane, however, closer to the edge of the plate, mixed fracture modes, for example for the $72^\circ$ and the $80^\circ$ model, cause a rather strong deviation.

Transition of failure modes with increasing fiber incline. Similar failure mechanisms are plotted above each other, (a) shows three different simulations.
(a) shows only phase field $d_2$, a crack that propagates along the fiber,
(b) shows the interaction of phase field $d_1$ and $d_2$, cracks that propagates along and perpendicular to the fiber
(c) shows only phase field $d_1$, a crack that propagates perpendicular to the fiber.

For the fiber angles shown in <Ref> (a) and (c), the results of the hybrid approach and the consistent approach agree.
The main difference appeared for the $72^\circ$ case shown in <Ref> (b).
In case of the consistent approach, no interaction occurred, meaning that the zig-zag pattern could not be reproduced (see the final state of the phase field in <Ref>).
The hybrid and the consistent approach differ in how the solid material is degraded, which influences the stress distribution and, thus, the crack driving forces.

To take a close look at these differences, a state prior to the fully developed phase field is compared in <Ref>.
Of primary interest are the stress components which only contribute to the development of one phase field, either $d_1$ or $d_2$, which for this quasi two-dimensional example are the Mode-I stresses, $\bm{\sigma}_{LL}$ for $d_1$ and $\bm{\sigma}_{RR}$ for $d_2$.
As $\bm{\sigma}_{RR}$ is degraded in both approaches, the main difference is in $\bm{\sigma}_{LL}$.
Looking at the stress plot in <Ref>, there is actually a notable difference between the hybrid and the consistent approach.
In the consistent approach the longitudinal stresses are not degraded, thus peaking at the geometrical crack tip from the single edge notched plate.
For the hybrid approach, with the longitudinal stresses being fully degraded, the stress peak moves with the diffuse phase field crack tip, causing a stress state at this location, which favors the evolution of phase field $d_1$.
This is an important finding, because it clearly shows that the usage of the hybrid approach is required in order to recover a zig-zag pattern.

Comparison of the hybrid and the consistent approach for two states, one prior to the fully developed phase field and one for the fully developed field.
The figure shows the main difference in the Mode-I stresses for phase field $d_1$ ($\bm{\sigma}_{LL}$) which leads to the zig-zag pattern for a $72^{\circ}$ fiber incline in case of the hybrid approach.
Stresses are given in $\si{MPa}$.

§.§ Comparison of the hybrid approach and the consistent approach

While <Ref> gives a strong argument in favor of the hybrid approach, it is still of interest, which of the two approaches is closer to modeling an actual discrete crack.
Therefore, results of both methods are compared with the resulting stress distribution and deformation of a model with a discretely modeled predefined crack.
The two different approaches, including a crack with a kink of $60^\circ$, are shown in <Ref>.
To properly compare the two cases, a discrete crack is modeled, and the same crack is modeled by solving the phase field problem for a Dirichlet-type boundary condition, prescribing $d_2 = 1$ on the same region.
Both cracks follow the fiber incline of $60^\circ$.
Subsequently, a linear elastic simulation for the model with the discrete crack and a simulation of the deformation problem considering the phase field distribution from <Ref>, using the hybrid and the consistent approach, was conducted.
For all three models the vertical deformation at the top edge is set to $u_\text{y} = 0.1$.

Two different approaches to compare a phase field crack with a geometrically modeled, discrete crack.
The phase field crack is computed by prescribing the Dirichlet-type boundary condition $d_2=1$ on the same region as the discrete crack is defined.

First, the deformation state of the three models is compared.
The results are shown in <Ref>, where only the outer edges of each of the notched plates are depicted.
Clearly, as the hybrid approach degrades all elastic components, it resembles the solution of the discrete crack model very well.
In contrast, in the case of the consistent approach, the remaining stresses related to phase field $d_1$ result in a quite large deviation from the discrete crack solution.

Comparison of the deformation of the outer edges of a model with a discrete crack and two phase field models, using the hybrid approach and the consistent approach.

<Ref> shows a comparison of stress distributions for the – essentially two-dimensional – notched plate.
As expected, the main differences between the three models occur for stress components related to the longitudinal direction ($\bm{\sigma}_{LL}$ and $\bm{\sigma}_{LR}$), as those are not fully degraded in the consistent approach.
<Ref>, region 1 shows that while the specimen is cracked at the given location, there are still stresses transferred through the crack.
Neither the hybrid approach, nor the model with the discrete crack show longitudinal stresses in this region.
Furthermore, this behavior also affects regions away from the crack, in region 2.
The influence on the remaining stress components is marginal, as $\bm{\sigma}_{RR}$ and $\bm{\sigma}_{LR}$ are fully degraded in both cases.
Nevertheless, in region 3, the consistent approach shows a peak in shear stresses at the kink which is considerably smaller at the other two models.

Comparison of three stress distributions between the consistent approach and the hybrid approach using a predefined phase field distribution and a discrete crack model.
1 shows that longitudinal stresses are still transferred through the crack when using the consistent approach and 2 shows that this does not only affect the region of the crack but also regions further away.
3 shows a peak in shear stresses at the kink, when using the consistent approach. Stresses are given in $\si{MPa}$.

Given the similarities of the stress distributions and deformation states, and the hybrid approach's ability to reproduce the zig-zag failure patter found in wood, the results give clear support for using the hybrid approach over the consistent one, when modeling complex material failure.

§.§ Wooden board with a knot

<Ref> shows the geometry of the wooden board with a single knot.
The board is supported both at the top and the right face.
The load is applied in form of a prescribed deformation along the bottom left edge.
In order to control crack initiation, the board has a 10mm notch through the specimen.

The fiber orientation $\bm{a}_1$ (longitudinal), $\bm{a}_2$ (radial), $\bm{a}_3$ (tangential), is computed for each individual integration point and defines the local material directions.
In this work, the model from lukacevic20193DModelKnots is used, where a knot is represented by a rotationally symmetric cone and wood fibers are streamlines flowing around an obstacle, which is the knot.
Consequently, the fiber orientation in the LT-plane can be computed using the so-called Rankine oval, which describes the fluid flow around an elliptical object.
Additionally, as fibers are situated on so-called growth surfaces, which motivate the cylindrical coordinate system commonly used to describe the elastic properties of wood, the third dimension of the fiber direction vector, the so-called dive angle, can be computed by restricting it to be orthogonal to the LT-plane.
<Ref> also shows a rendering of the fiber course.
The texture is generated based on the previously computed fiber directions on a $0.10\si{mm}$ grid in the mid-plane of the board.

Similar to the notched plate from <Ref>, the mesh size is homogenous in the left part of the specimen, where the crack will open.
Restricting phase field evolution just to the left part reduces the number of DOF and keeping the mesh homogenous reduces the influence of the mesh structure on the crack paths.
The characteristic length of the elements in the phase field activated region is set to $0.75\si{mm}$.
Relative to the volume of this region, this matches the mesh density of notched plate model number 6 in <Ref>, which is in good agreement with the results from an even smaller characteristic length.
The resulting total number of DOF is $976,200$.
As outlined in <Ref>, the hybrid approach is required for properly modeling fracture processes of wood.
Therefore, for simulation of this more complex example, only the hybrid approach is used.
The parameters controlling the phase field problem are given in <Ref>.

Geometry of the wooden board with a single knot.
The fiber orientation $\bm{a}_1$ (longitudinal), $\bm{a}_2$ (radial), $\bm{a}_3$ (tangential) is prescribed in each integration point and is rendered in the mid-plane of the board.
A notch is placed in order to control crack initiation.
The board is loaded at the bottom left edge by prescribing the deformation.
All measurements are in mm.

Defining parameters for the wooden board with a single knot.

$d_i$     $\beta_i$a      $G_{c,i}$b     $f_{t,i}$c     $\sfrac{l_i}{l_{\text{eff}}}$

$d_1$     2.0     2.0      80.0      2.0

$d_2$     4.0     0.8      5.0     2.0

$d_3$     2.0     0.1       3.6     2.0

a Structural tensor scale in <Ref>

b in $\sfrac{\si{Nmm}}{mm^2}$

c in $\si{MPa}$

<Ref> shows the resulting crack path when the specimen is almost fully cracked.
Phase field $d_2$ is visualized using three-dimensional contour lines ranging from $d_2 = 0.0$ to $d_2 = 1.0$ in increments of $0.1$.
The fiber direction is depicted by plotting $\bm{a}_1$ on a uniformly spatially distributed subset of integration points.
Obviously, the varying fiber directions influence the orientation of the crack face.
<Ref> shows the evolution of the phase field variable $d_2$.
The crack initially starts with a slight decline and changes its orientation in the vicinity of the knot, where the fibers become parallel to the knot's surface.
When reaching the knot, the crack kinks and follows the weak interface.
It stops propagating close to the lower edge of the board, as the compressive Mode-I stresses in this region do not result in crack driving forces, due to the additive decomposition of the strain energy density term.

Fully cracked wooden board.
The fiber orientation vector $\bm{a}_1$ is plotted on a uniformly spatially distributed subset of the integration points.
The crack mostly follows the fiber direction, starting horizontally, tilts and subsequently propagates along the interface region between the clear wood area and the knot.

<Ref> shows the load-deflection plot of the simulation, measured at the lower left edge of the board.
The horizontal axis is split into two differently scaled parts, as the change in the reaction force from 03mm is quite large compared to the change from 1030mm.
Past the initial opening of the crack, the load-deflection plot shows a cohesive behavior during further crack propagation.
The softening effect can be controlled by properly setting $f_{t,i}$ and $G_{c,i}$ in <Ref>.
The reaction force is heavily reduced while the crack propagates along the fiber towards the knot.
With the crack further progressing, tensile and compressive stresses, similar to the bending stresses at the clamped end of a cantilever beam, concentrate at the lower left edge of the knot.
This shift in the stress distribution results in a less stiff response of the system, therefore, larger deformations are required for further crack growth along the weak interface between the clear wood area and the knot.

Load-deflection plot measured at the lower left edge of the board.
Additionally, the evolution of the phase field variable $d_2$ is shown for multiple locations along the path.

The simulation of a more complex model considering a realistic fiber course showed that the discussed phase field model is capable of considering the effects resulting from fiber deviations and that sudden changes in the crack face orientation, the kink when the crack reaches the knot, can be modeled.
Furthermore, using the unified phase field theory, adapted to a linear softening law, allows for a cohesive behavior during crack propagation.

§ CONCLUSION AND OUTLOOK

The present work addresses the formulation of a phase field model for orthotropic non-brittle materials, able to reproduce multiple, very different failure mechanisms.
In order to extend the phase field method for fracture to support cohesive behavior, the so-called unified phase field theory is applied and tuned to a linear softening law, resembling a cohesive zones model.
Subsequently, a stress-based split for anisotropic materials is derived, which is based on considering Mode-I, Mode-II and Mode-III crack driving stresses on a fictitious crack plane.
The plane's orientations are defined accordingly to material specific fracture planes, in this case for wood: a crack plane perpendicular to the longitudinal, the radial and the tangential direction.
The orientation in which crack growth results in the largest energy dissipation is the driving failure mechanism.

This formulation is coupled in form of a hybrid approach, by separating the energetic driving force term and the actual degradation of the solid.
In this novel hybrid approach, a smooth traction-free crack boundary condition is used, which incorporates a contact constraint and, thus, does not require an additional constraint for preventing interpenetration of crack faces.
This concept is then put into a multi-phase field model, which allows defining a different fracture behavior for each phase field variable individually.
Therefore, very different failure mechanisms can be modeled and described realistically.
In order to consider the effect of the material's structure on the crack paths, a second-order tensor is added to the crack density function, which scales the phase field's gradient on the plane perpendicular to the associated crack orientation vector.
Hence, preferable planes for crack propagation can be defined, a crack perpendicular to the radial direction is likely to propagate along the fiber direction (longitudinal) and less likely to propagate in the radial direction, due to the weak interface between the fibers and the matrix.

The outlined method is tested using two numerical examples of different complexity, both being wooden specimens.
By means of the model of a sinlge edge notched plate, it is shown that changing the fiber orientation leads to different crack topologies, where cracks travel along the fiber when the load direction is in an obtuse angle relative to the fiber direction.
At a certain fiber incline, as expected, the crack kinks and jumps to the next fictitious growth layer, rupturing the fibers in between.
At a sharp angle, the failure mode changes to a crack perpendicular to the fiber orientation.
The simulations showed that such common phenomena of wood (the zig-zag pattern) cannot be recovered when a variationally consistent  approach is used, thus motivating the use of the hybrid approach with a smooth traction free crack boundary condition.
Subsequently, a more complex example of a wooden board with a single knot and a spatially varying fiber orientation was tested.
The model shows that the phase field crack actually follows the curvature of the wood fibers and also allows for sudden changes in the crack face orientation, in the vicinity of the knot where the crack kinks.
Furthermore, the influence of the cohesive behavior during crack propagation can be observed.

This allows the conclusion that the phase field method can be used to model wood failure, as crack phenomena like the zig-zag pattern can be modeled, complex crack topologies can be depicted and cohesive behavior can be considered.
An apparent limitation of this work lies in the formulation of the energetic driving force, which, while allowing the definition of a different fracture characteristic on the level of each phase field variable, allows no distinction between Mode-I, Mode-II and Mode-III.
Thus, always mixed mode failure is assumed.
As this study's focus is on the implementation of a phase field model for wood and investigation of commonly found crack patterns, future research on validating the model with experimental data is needed.
Furthermore, examining more complex examples like wooden boards with multiple knots is of interest.

§ ACKNOWLEDGEMENT

This research was funded in whole, or in part, by the Austrian Science Fund (FWF) . For the purpose of open access, the author has applied a CC BY public copyright licence to any Author's Accepted Manuscript version arising from this submission.
The authors also acknowledge gratefully the support by the ForestValue project InnoCrossLam.

§ COMPARISON OF THE ACTIVE SET REDUCED SPACE METHOD AND THE PRIMAL-DUAL ACTIVE SET METHOD

The only difference between the active set reduced space method from yang2016a, and the primal-dual active set method from Heister2015, is in the selection of the active and the inactive set.
Generally speaking, whether a phase field DOF is in the active or the inactive set, is determined by two aspects, the current value of the phase field and the current value of the residual.
Adaption to the notation used in this work and rearranging the parts in the primal-dual active set formulation, leads to the following definition of the primal-dual active set:

\begin{equation}
\mathcal{A}(\dFE^k_{n+1}) = \left\{ i \in \mathcal{S} | c \mathbf{B}_{ii}\left(\dFE^{k-1} - \dFE^k_{n+1} \right)_i  + (\bm{\hat{R}}^k_{\dFE,n+1})_i > 0 \right\},
\label{eq:heister-active-set}
\end{equation}

where $c$ is a constant larger than 0 and $\mathbf{B}_{ii}$ is the entry of the $i$-th DOF in the diagonal mass matrix $\mathbf{B}$.
Given that both $c$ and the mass are strictly larger than 0, their product is as well.
Comparing the two approaches from <Ref> leads to the following observations:

* If both $\left(\dFE^{k-1} - \dFE^k_{n+1}\right)_i$  and $(\bm{\hat{R}}^k_{\dFE,n+1})_i$ are larger than 0, the $i$-th DOF is in both formulations considered to be active.

* If both $\left(\dFE^{k-1} - \dFE^k_{n+1}\right)_i$  and $(\bm{\hat{R}}^k_{\dFE,n+1})_i$ are smaller than 0, the $i$-th DOF is in both formulations considered to be inactive.

* If the sign of the two terms is different, the $i$-th DOF is inactive in the method from yang2016a, however, in the method from Heister2015 it depends on the choice of the constant $c$, which is not further elaborated in their work.

Missing 'biblatex' package
The bibliography requires the 'biblatex' package.

Springer Science and Business Media LLC

journaltitleComputational Mechanics
titlePhase-Field Modeling of Ductile Fracture
spech/Google Drive/Zotero/Ambati_et_al_2015_Phase-field_modeling_of_ductile_fracture.pdf
phase field,phase_field_implementation,phase_field_implementation_paper

Springer Science and Business Media LLC

journaltitleComput. Mech.
titleA Review on Phase-Field Models of Brittle Fracture and a New Fast Hybrid Formulation
spech/Google Drive/Zotero/Ambati_et_al_2014_A_review_on_phase-field_models_of_brittle_fracture_and_a_new_fast_hybrid.pdf
phase field,phase_field_implementation,phase_field_implementation_paper,reading list

abstractWe show how it is possible to approximate the Mumford-Shah (see [29]) image segmentation functional by elliptic functionals defined on Sobolev spaces. The heuristic idea is to consider functionals 𝒢h(u, z) with z ranging between 0 and 1 and related to the set K. The minimizing zh are near to 1 in a neighborhood of the set K, and far from the neighborhood they are very small. The neighborhood shrinks as h → + ∞. For a similar approach to the problem compare Kulkarni; see [25]. The approximation of 𝒢h to 𝒢 takes place in a variational sense, the De Giorgi F-convergence.
annotation_eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/cpa.3160430805
journaltitleCommunications on Pure and Applied Mathematics
titleApproximation of Functional Depending on Jumps by Elliptic Functional via T-Convergence
spech/Google Drive/Zotero/Ambrosio_Tortorelli_1990_Approximation_of_functional_depending_on_jumps_by_elliptic_functional_via.pdf;/Users/spech/Zotero/storage/U8CEWFN7/cpa.html

Elsevier BV

journaltitleJournal of the Mechanics and Physics of Solids
titleRegularized Formulation of the Variational Brittle Fracture with Unilateral Contact: Numerical Experiments
spech/Google Drive/Zotero/Amor_et_al_2009_Regularized_formulation_of_the_variational_brittle_fracture_with_unilateral.pdf

booktitleAdvances in Applied Mechanics
titleThe Mathematical Theory of Equilibrium Cracks in Brittle Fracture
spech/Google Drive/Zotero/Barenblatt_1962_The_mathematical_theory_of_equilibrium_cracks_in_brittle_fracture.pdf

abstractTriangular and prismatic quadratic isoparametric elements, formed by collapsing one side and placing the mid-side node near the crack tip at the quarter point, are shown to embody the (1/√r) singularity of elastic fracture mechanics and the (1/r) singularity of perfect plasticity. The procedure of performing the fracture analysis for the case of small scale yielding is discussed, and the finite element results are compared with theoretical results. The proposed elements have wide application in the fracture analysis of structures where ductile fracture is investigated. They permit a determination of the relationship between crack tip field parameters, loading, and geometry. And for a given fracture criterion can be applied to the prediction of fracture in structures such as pressure vessels under in service conditions.
annotation_eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.1620110109
journaltitleInternational Journal for Numerical Methods in Engineering
titleTriangular Quarter-Point Elements as Elastic and Perfectly-Plastic Crack Tip Elements

journaltitleSIAM Rev.
titleJulia: A Fresh Approach to Numerical Computing
spech/Google Drive/Zotero/Bezanson_et_al_2017_Julia.pdf;/Users/spech/Google Drive/Zotero/Bezanson_et_al_2017_Julia2.pdf
65Y05,68N15,97P40,Julia,numerical,parallel,phase_field_implementation,scientific computing

journaltitleComput. Methods Appl. Mech. Eng.
titlePhase-Field Modeling of Anisotropic Brittle Fracture Including Several Damage Mechanisms
spech/Google Drive/Zotero/Bleyer_Alessi_2018_Phase-field_modeling_of_anisotropic_brittle_fracture_including_several_damage.pdf;/Users/spech/Google Drive/Zotero/Bleyer_Alessi_2018_Phase-field_modeling_of_anisotropic_brittle_fracture_including_several_damage2.pdf
ase field,phase_field_implementation,phase_field_implementation_paper,Phase-field,reading list,Variational model

Elsevier BV

journaltitleJournal of the Mechanics and Physics of Solids
titleNumerical Experiments in Revisited Brittle Fracture
spech/Google Drive/Zotero/Bourdin_et_al_2000_Numerical_experiments_in_revisited_brittle_fracture.pdf
phase field,phase_field_implementation,phase_field_implementation_paper,reading list

Springer-Verlag GmbH

titleThe Variational Approach to Fracture
Brittleness,Calculus of variations,Fracture mechanics,Mathematical models,phase field,phase_field_implementation,phase_field_implementation_paper,reading list

abstractAbstract Genomic prediction relies on genotypic marker information to predict the agronomic performance of future hybrid breeds based on trial records. Because the effect of markers may vary substantially under the influence of different environmental conditions, marker-by-environment interaction effects have to be taken into account. However, this may lead to a dramatic increase in the computational resources needed for analyzing large-scale trial data. A high-performance computing solution, called Needles, is presented for handling such data sets. Needles is tailored to the particular properties of the underlying algebraic framework by exploiting a sparse matrix formalism where suited and by utilizing distributed computing techniques to enable the use of a dedicated computing cluster. It is demonstrated that large-scale analyses can be performed within reasonable time frames with this framework. Moreover, by analyzing simulated trial data, it is shown that the effects of markers with a high environmental interaction can be predicted more accurately when more records per environment are available in the training data. The availability of such data and their analysis with Needles also may lead to the discovery of highly contributing QTL in specific environmental conditions. Such a framework thus opens the path for plant breeders to select crops based on these QTL, resulting in hybrid lines with optimized agronomic performance in specific environmental conditions.
titleNeedles: Toward Large-Scale Genomic Prediction with Marker-by-Environment Interaction
spech/Google Drive/Zotero/De_Coninck_et_al_2016_Needles.pdf

Springer International Publishing

booktitleModeling in Engineering Using Innovative Numerical Methods for Solids and Fluids
seriesCISM International Centre for Mechanical Sciences
titleNumerical Implementation of Phase-Field Models of Brittle Fracture
spech/Google Drive/Zotero/De_Lorenzis_Gerasimov_2020_Numerical_implementation_of_phase-field_models_of_brittle_fracture.pdf

abstractIn phase-field theories of fracture, decompositions of the strain energy density into tensile and compressive parts are often necessary to avoid interpenetration of cracked surfaces and to select physically trustworthy crack paths. General formulations accounting for orthotropy of the well-known spectral and hydrostatic-deviatoric decompositions of the strain tensor (often referred to as Miehe and Amor decompositions) are presented in this study. Additionally, a new principal energy decomposition based on spectral decomposition of the stiffness tensor is proposed for general anisotropic materials. The decompositions are evaluated numerically in a quadratic specimen with an initial stationary edge crack subject to both tensile and shear global remote loading. It is shown that when an isotropic case is considered, solutions agree well with results reported elsewhere for both the spectral and the hydrostatic-deviatoric approaches. The principal energy decomposition results in similar crack paths as the other approaches, with only subtle differences. When orthotropy is considered, however, significant differences in the resulting crack paths as well as global force-displacement behavior are obtained, especially when the crack is subjected to shear loading. For global tensile loading, the decompositions result in similar crack paths and force-displacement relations. The results provide a step forward when developing phase-field fracture theories for brittle materials with an orthotropic nature and highlight the importance of a proper decomposition of the strain energy density.
journaltitleInt J Solids Struct
titleStrain Energy Density Decompositions in Phase-Field Fracture Theories for Orthotropy and Anisotropy
spech/Google Drive/Zotero/Dijk_et_al_2020_Strain_energy_density_decompositions_in_phase-field_fracture_theories_for.pdf;/Users/spech/Google Drive/Zotero/Dijk_et_al_2020_Strain_energy_density_decompositions_in_phase-field_fracture_theories_for2.pdf;/Users/spech/Zotero/storage/YP2UTEMI/S0020768320301396.html
s,Finite element,Fracture,Phase field,phase_field_implementation,phase_field_implementation_paper,Strain energy decomposition

Elsevier BV

journaltitleComposites Part A: Applied Science and Manufacturing
titleComparison of Fracture Properties of Two Wood Species through Cohesive Crack Simulations
spech/Google Drive/Zotero/Dourado_et_al_2008_Comparison_of_fracture_properties_of_two_wood_species_through_cohesive_crack.pdf

Elsevier BV

journaltitleJournal of the Mechanics and Physics of Solids
titleRevisiting Brittle Fracture as an Energy Minimization Problem
spech/Google Drive/Zotero/Francfort_Marigo_1998_Revisiting_brittle_fracture_as_an_energy_minimization_problem.pdf
phase field,phase_field_implementation,phase_field_implementation_paper,reading list

abstractPhase-field modeling of fracture phenomena in solids is a very promising approach which has gained popularity within the last decade. However, within the finite element framework, already a two-dimensional quasi-static phase-field formulation is computationally quite demanding, mainly for the following reasons: (i) the need to resolve the small length scale inherent to the diffusive crack approximation calls for extremely fine meshes, at least locally in the crack phase-field transition zone, (ii) due to non-convexity of the related free-energy functional, a robust, but slowly converging staggered solution scheme based on algorithmic decoupling is typically used. In this contribution we tackle problem (ii) and propose a faster and equally accurate approach for quasi-static phase-field computing of (brittle) fracture using a monolithic solution scheme which is accompanied by a novel line search procedure to overcome the iterative convergence issues of non-convex minimization. We present a detailed critical evaluation of the approach and its comparison with the staggered scheme in terms of computational cost, accuracy and robustness.
journaltitleComputer Methods in Applied Mechanics and Engineering
seriesPhase Field Approaches to Fracture
titleA Line Search Assisted Monolithic Approach for Phase-Field Computing of Brittle Fracture
spech/Google Drive/Zotero/Gerasimov_De_Lorenzis_2016_A_line_search_assisted_monolithic_approach_for_phase-field_computing_of_brittle.pdf;/Users/spech/Zotero/storage/WIFUENKB/S0045782515004235.html

Elsevier BV

journaltitleComputer Methods in Applied Mechanics and Engineering
titleOn Penalization in Variational Phase-Field Models of Brittle Fracture
Mathematics - Numerical Analysis,phase field,phase_field_implementation,phase_field_implementation_paper,reading list

abstractIn the course of an investigation of the effect of surface scratches on the mechanical strength of solids, some general conclusions were reached which appear to have a direct bearing on the problem of rupture, from an engineering standpoint, and also on the larger question of the nature of intermolecular cohesion. The original object of the work, which was carried out at the Royal Aircraft Establishment, was the discovery of the effect of surface treatment—such as, for instance, filing, grinding or polishing—on the strength of metallic machine parts subjected to alternating or repeated loads. In the case of steel, and some other metals in common use, the results of fatigue tests indicated that the range of alternating stress which could be permanently sustained by the material was smaller than the range within which it was sensibly elastic, after being subjected to a great number of reversals. Hence it was inferred that the safe range of loading of a part, having a scratched or grooved surface of a given type, should be capable of estimation with the help of one of the two hypotheses of rupture commonly used for solids which are elastic to fracture. According to these hypotheses rupture may be expected if (a) the maximum tensile stress, (b) the maximum extension, exceeds a certain critical value. Moreover, as the behaviour of the materials under consideration, within the safe range of alternating stress, shows very little departure from Hooke’s law, it was thought that the necessary stress and strain calculations could be performed by means of the mathematical theory of elasticity.
journaltitlePhilosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character
titleVI. The Phenomena of Rupture and Flow in Solids

journaltitleComput. Methods Appl. Mech. Eng.
titleNumerical Aspects of Anisotropic Failure in Soft Biological Tissues Favor Energy-Based Criteria: A Rate-Dependent Anisotropic Crack Phase-Field Model
spech/Google Drive/Zotero/Gültekin_et_al_2018_Numerical_aspects_of_anisotropic_failure_in_soft_biological_tissues_favor.pdf
ase-field,Failure criteria,Fracture,phase_field_implementation,phase_field_implementation_paper,Soft biological tissues

Elsevier BV

journaltitleComputer Methods in Applied Mechanics and Engineering
titleA Primal-Dual Active Set Method and Predictor-Corrector Mesh Adaptivity for Computing Fracture Propagation Using a Phase-Field Approach
spech/Google Drive/Zotero/Heister_et_al_2015_A_primal-dual_active_set_method_and_predictor-corrector_mesh_adaptivity_for.pdf
phase field,phase_field_implementation,phase_field_implementation_paper,reading list,solver

Society for Industrial & Applied Mathematics (SIAM)

journaltitleSIAM Journal on Optimization
titleThe Primal-Dual Active Set Strategy as a Semismooth Newton Method
spech/Google Drive/Zotero/Hintermüller_et_al_2002_The_primal-dual_active_set_strategy_as_a_semismooth_newton_method.pdf
phase_field_implementation,phase_field_implementation_paper,reading list

abstractWe present a new derivation for a phase-field model of cohesive fracture that allows for fully-damaged surfaces to properly transmit tractions under frictionless contact conditions. The model is derived from an energy minimization standpoint, and the governing equations are presented in a general Allen–Cahn form, unifying both brittle and cohesive fracture models. A novel elastic energy split is proposed to enforce frictionless contact conditions along the regularized crack set. A fixed-point iterative algorithm is used to solve the system of equations, and an active-set strategy is incorporated to enforce the crack set irreversibility. The model is then applied to study the mechanical-fracture coupling in the fracture of thin films and soil desiccation problems in which intricate crack networks form. We present results for one-dimensional, two-dimensional, and fully three-dimensional configurations, and examine the relationship between the generalized driving energies for fracture and the resulting fracture patterns. We also investigate the impact of stochastic spatial variations in material parameters on fragment statistics and fracture network morphology. To this end, we develop a random field model for fracture properties and address a forward propagation problem. The versatility offered by the probabilistic framework is finally highlighted by considering a statistical inverse problem based on physical experiments. It is shown, in particular, that the formulation allows key characteristics of observed crack patterns to be reproduced with reasonable accuracy.
journaltitleComputer Methods in Applied Mechanics and Engineering
shorttitleA Phase-Field Model of Fracture with Frictionless Contact and Random Fracture Properties
titleA Phase-Field Model of Fracture with Frictionless Contact and Random Fracture Properties: Application to Thin-Film Fracture and Soil Desiccation
spech/Google Drive/Zotero/Hu_et_al_2020_A_phase-field_model_of_fracture_with_frictionless_contact_and_random_fracture.pdf;/Users/spech/Google Drive/Zotero/Hu_et_al_2020_A_phase-field_model_of_fracture_with_frictionless_contact_and_random_fracture2.pdf;/Users/spech/Zotero/storage/G4NKKE7Z/S0045782520302905.html
ase_field_implementation,phase_field_implementation_paper,Phase-field models,Random field,reading,Soil desiccation,Traction-free surfaces,Uncertainty quantification

Springer Berlin Heidelberg

booktitleElasticity and Plasticity / Elastizität Und Plastizität
spech/Google Drive/Zotero/Irwin_1958_Fracture.pdf

abstractThe simulation of crack initiation and propagation in an elastic material is difficult, as crack paths with complex topologies have to be resolved. Phase-field approaches allow to simulate crack behavior without the need to explicitly model crack paths. However, the underlying mathematical model gives rise to a non-convex constrained minimization problem. In this work, we propose a recursive multilevel trust region (RMTR) method to efficiently solve such a minimization problem. The RMTR method combines global convergence properties of trust region methods with the optimality of multilevel methods. The solution process is accelerated by employing level dependent objective functions, minimization of which provides correction to the original/fine-level problem. In the context of the phase-field fracture approach, it is challenging to design efficient level dependent objective functions as a certain mesh resolution is required to resolve fracture bands. We introduce level dependent objective functions that combine fine level description of the crack path with the coarse level discretization. The overall performance and the convergence properties of the proposed RMTR method are investigated by means of several numerical examples in three dimensions.
journaltitleComputer Methods in Applied Mechanics and Engineering
titleA Recursive Multilevel Trust Region Method with Application to Fully Monolithic Phase-Field Models of Brittle Fracture
spech/Google Drive/Zotero/Kopaničáková_Krause_2020_A_recursive_multilevel_trust_region_method_with_application_to_fully_monolithic.pdf;/Users/spech/Zotero/storage/5523IAMT/S0045782519306085.html
ic scheme,Multilevel,Non-convex functional,phase_field_implementation,phase_field_implementation_paper,Phase-field fracture,reading-list,Trust-region

abstractDistributed energy storage devices are commonly employed as an effective approach for addressing operational challenges introduced by the large scale integration of renewables. However, the modeling of storage devices results in intertemporal coupling of the individual optimal power flow (OPF) problems defined at each subdivision of the time period of interest. The resulting multiperiod optimal power flow (MPOPF) problem becomes intractable prohibiting, forecasting, and planning over long time periods. Interior point (IP) methods have been extensively employed for the solution of OPF and MPOPF problems. This work proposes an efficient IP algorithm, BELTISTOS, particularly designed for MPOPF problems. The structure of the linear system associated with the Karush-Kuhn-Tucker conditions is revisited, and a Schur-complement-based approach tailored to its structure is proposed. Through benchmark cases involving power-grid models of increasing complexity, the BELTISTOS algorithm is demonstrated to provide several orders of magnitude faster solution times than standard optimization methods, such as IPOPT, MIPS, and KNITRO, using significantly less memory.
eventtitleIEEE Transactions on Power Systems
journaltitleIEEE Transactions on Power Systems
titleToward the Next Generation of Multiperiod Optimal Power Flow Solvers
Artificial neural networks,Couplings,Generators,interior point methods,IP networks,Multiperiod optimal power flow,Niobium,Optimization,phase_field_implementation,Planning,power system planning

Elsevier BV

journaltitleComputational Materials Science
titleOn Degradation Functions in Phase Field Fracture Models
spech/Google Drive/Zotero/Kuhn_et_al_2015_On_degradation_functions_in_phase_field_fracture_models.pdf

abstractIncreased use of wood has led to complex timber constructions and new types of engineered wood products. In simulations, however, mainly simplified models are used to describe this material with its strongly varying properties. Therefore, reliable prediction tools for mechanical properties of wooden boards are needed. Those varying properties mainly originate from knots and fiber deviations. Thus, we use fiber directions on board surfaces to reconstruct knots within boards. Combined with a fiber deviation model we assess our model with experiments on different levels: fiber directions on surfaces, strain fields and bending stiffness profiles. This model now better describes fiber patterns near knots and knot clusters. Also, we showed that accurate modeling of the pith is important to avoid large regions of incorrect fiber deviations. Furthermore, modified knot stiffness properties were successfully used to consider pre-cracked knots. Finally, we obtained multiple bending stiffness profiles, where we showed that even local effects can be simulated accurately. We anticipate our tool to be a starting point for improving strength grading models, where effects of knot configurations can be studied more easily than with experiments alone. Furthermore, the presented improvements will render the simulation of realistic failure mechanisms in wooden boards more likely.
journaltitleMaterials & Design
titleA 3D Model for Knots and Related Fiber Deviations in Sawn Timber for Prediction of Mechanical Properties of Boards
spech/Google Drive/Zotero/Lukacevic_et_al_2019_A_3D_model_for_knots_and_related_fiber_deviations_in_sawn_timber_for_prediction.pdf;/Users/spech/Zotero/storage/TY8DBNDP/S0264127519300541.html
anical properties,Fiber course,Knots,phase_field_implementation,phase_field_implementation_paper,Wood

abstractAbstract The computational modeling of failure mechanisms in solids due to fracture based on sharp crack discontinuities suffers in situations with complex crack topologies. This can be overcome by a diffusive crack modeling based on the introduction of a crack phase-field. In this paper, we outline a thermodynamically consistent framework for phase-field models of crack propagation in elastic solids, develop incremental variational principles and consider their numerical implementations by multi-field finite element methods. We start our investigation with an intuitive and descriptive derivation of a regularized crack surface functional that Γ-converges for vanishing length-scale parameter to a sharp crack topology functional. This functional provides the basis for the definition of suitable convex dissipation functions that govern the evolution of the crack phase-field. Here, we propose alternative rate-independent and viscous over-force models that ensure the local growth of the phase-field. Next, we define an energy storage function whose positive tensile part degrades with increasing phase-field. With these constitutive functionals at hand, we derive the coupled balances of quasi-static stress equilibrium and gradient-type phase-field evolution in the solid from the argument of virtual power. Here, we consider a canonical two-field setting for rate-independent response and a time-regularized three-field formulation with viscous over-force response. It is then shown that these balances follow as the Euler equations of incremental variational principles that govern the multi-field problems. These principles make the proposed formulation extremely compact and provide a perfect base for the finite element implementation, including features such as the symmetry of the monolithic tangent matrices. We demonstrate the performance of the proposed phase-field formulations of fracture by means of representative numerical examples. Copyright © 2010 John Wiley & Sons, Ltd.
journaltitleInt. J. Numer. Methods Eng.
titleThermodynamically Consistent Phase-Field Models of Fracture: Variational Principles and Multi-Field FE Implementations
/Google Drive/Zotero/Miehe_et_al_2010_Thermodynamically_consistent_phase-field_models_of_fracture.pdf
coupled multi-field problems,crack propagation,finite elements,fracture,gradient-type damage,incremental variational principles,phase_field_implementation,phase_field_implementation_paper,phase-fields

journaltitleComput. Methods Appl. Mech. Eng.
titleA Phase Field Model for Rate-Independent Crack Propagation: Robust Algorithmic Implementation Based on Operator Splits
spech/Google Drive/Zotero/Miehe_et_al_2010_A_phase_field_model_for_rate-independent_crack_propagation.pdf
ase fields,phase_field_implementation,phase_field_implementation_paper

journaltitleComput. Methods Appl. Mech. Eng.
titlePhase Field Modeling of Fracture in Multi-Physics Problems. Part I. Balance of Crack Surface and Failure Criteria for Brittle Crack Propagation in Thermo-Elastic Solids
spech/Google Drive/Zotero/Miehe_et_al_2015_Phase_field_modeling_of_fracture_in_multi-physics_problems.pdf
ase field modeling,phase_field_implementation,phase_field_implementation_paper,Thermo-elasticity

journaltitleEngineering Fracture Mechanics
titleExtended Finite Element Method for Cohesive Crack Growth
spech/Google Drive/Zotero/Moës_Belytschko_2002_Extended_finite_element_method_for_cohesive_crack_growth.pdf

journaltitleEur. J. of Mech. - A/Solids
titlePhase Field Modelling of Anisotropic Crack Propagation
spech/Google Drive/Zotero/Nguyen_et_al_2017_Phase_field_modelling_of_anisotropic_crack_propagation.pdf
ase field,phase_field_implementation,phase_field_implementation_paper

abstractA node split method based on Partition-of-unity method and mesh modified has been developed for numerical simulation of crack growth. The crack tip could be represented by the enrichment function based on Partition-of-unity method. The crack besides the tip always locates on the edge of element by moving the nodes around crack, and the continuity of crack could be kept on. For an example of 2D quadrilateral mesh, mesh processing has been presented on different condition. The simulation results indicate that applying the node split method to simulate crack growth problem can achieve relatively good results even for sparse grid.
titleA Node Split Method for Crack Growth Problem
spech/Google Drive/Zotero/Peng_Wang_2012_A_Node_Split_Method_for_Crack_Growth_Problem.pdf

SAGE Publications

journaltitleInternational Journal of Damage Mechanics
titleGradient Damage Models and Their Use to Approximate Brittle Fracture
spech/Google Drive/Zotero/Pham_et_al_2010_Gradient_damage_models_and_their_use_to_approximate_brittle_fracture.pdf

journaltitleEngineering Fracture Mechanics
titleCohesive-Zone Modelling of Crack Nucleation and Propagation in Particulate Composites
spech/Google Drive/Zotero/Ponnusami_et_al_2015_Cohesive-zone_modelling_of_crack_nucleation_and_propagation_in_particulate.pdf

abstractWe present ForwardDiff, a Julia package for forward-mode automatic differentiation (AD) featuring performance competitive with low-level languages like C++. Unlike recently developed AD tools in other popular high-level languages such as Python and MATLAB, ForwardDiff takes advantage of just-in-time (JIT) compilation to transparently recompile AD-unaware user code, enabling efficient support for higher-order differentiation and differentiation using custom number types (including complex numbers). For gradient and Jacobian calculations, ForwardDiff provides a variant of vector-forward mode that avoids expensive heap allocation and makes better use of memory bandwidth than traditional vector mode. In our numerical experiments, we demonstrate that for nontrivially large dimensions, ForwardDiff's gradient computations can be faster than a reverse-mode implementation from the Python-based autograd package. We also illustrate how ForwardDiff is used effectively within JuMP, a modeling language for optimization. According to our usage statistics, 41 unique repositories on GitHub depend on ForwardDiff, with users from diverse fields such as astronomy, optimization, finite element analysis, and statistics. This document is an extended abstract that has been accepted for presentation at the AD2016 7th International Conference on Algorithmic Differentiation.
titleForward-Mode Automatic Differentiation in Julia
spech/Google Drive/Zotero/Revels_et_al_2016_Forward-Mode_Automatic_Differentiation_in_Julia.pdf;/Users/spech/Zotero/storage/KVPR7R6C/1607.html
ematical Software,phase_field_implementation

abstractThis paper presents a finite element analysis based on the remeshing technique to predict the dynamic crack propagation and crack arrest in a brittle material, namely Araldite-B. The dynamic fracture toughness criterion is used to form the crack tip equation of motion. Plane stress condition is invoked in the present two-dimensional fracture analysis of the RDCB (rectangular double cantilever beam) specimens and the unknown crack tip position, and velocity is computed during the analysis by the dynamic fracture criterion. According to the new crack tip position, a remeshing algorithm has been used to simulate the dynamic crack growth and arrest. The obtained results including kinetic energy and strain energy, crack tip velocity, dynamic stress intensity factor during crack growth and also crack arrest length have been presented. Comparison of the results with those cited in the literature has shown a good agreement.
journaltitleMaterials & Design
titleFinite Element Analysis of Dynamic Crack Propagation Using Remeshing Technique
spech/Google Drive/Zotero/Shahani_Amini_Fasakhodi_2009_Finite_element_analysis_of_dynamic_crack_propagation_using_remeshing_technique.pdf;/Users/spech/Zotero/storage/ZU7PXBYN/S0261306908003208.html

Chichester, West Sussex, England ; Hoboken, NJ

J. Wiley

titleFracture and Fatigue in Wood

journaltitleComput. Mech.
titleA Phase-Field Crack Model Based on Directional Stress Decomposition

abstractIn the phase field approach for fracture an additional scalar field is introduced in order to describe the state of the material between intact and fully broken. So far, for the loading dependent degradation of stiffness (damage) either the volumetric-deviatoric split of strain [1, 2] or the spectral decomposition [3, 4] is used. In contrast to such an isotropic degradation of stiffness, the fully broken state represents a crack with a particular orientation. Both aforementioned approaches do not take the crack orientation into account. This may lead to the violation of the crack boundary conditions. In order to satisfy these conditions the phase field approach is modified here by taking the orientation of the crack into account. (© 2015 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim)
annotation_eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/pamm.201510068
titleA Novel Treatment of Crack Boundary Conditions in Phase Field Models of Fracture
spech/Google Drive/Zotero/Strobl_Seelig_2015_A_novel_treatment_of_crack_boundary_conditions_in_phase_field_models_of_fracture.pdf;/Users/spech/Google Drive/Zotero/Strobl_Seelig_2015_A_novel_treatment_of_crack_boundary_conditions_in_phase_field_models_of_fracture2.pdf;/Users/spech/Zotero/storage/RP6SAJ5Z/pamm.html

journaltitleInt. J. Non Linear Mech.
titlePhase Field Modeling of Fracture in Anisotropic Brittle Solids
spech/Google Drive/Zotero/Teichtmeister_et_al_2017_Phase_field_modeling_of_fracture_in_anisotropic_brittle_solids.pdf
ase field,Phase field modeling,phase_field_implementation,phase_field_implementation_paper,reading list

journaltitleWood and Fiber Science
titleFracture Zone Characterization - Micro-mechanical Study

abstractA parallel distributed-memory approach for the exact calculation of selected entries of the inverse of a matrix arising in a Best Linear Unbiased Estimation (BLUE) problem in genomic prediction is presented. The particular structure of the matrices involved in this stochastic process, consisting of sparse and dense blocks, requires a framework coupling sparse and dense linear algebra algorithms. Our approach exploits direct sparse techniques based on the Takahashi equations, coupled with distributed LU dense factorizations and Schur-complement computations. The algorithm is validated on several matrices on a Cray XC40 supercomputer.
journaltitleJournal of Computational Science
titleEnhancing the Scalability of Selected Inversion Factorization Algorithms in Genomic Prediction
spech/Google Drive/Zotero/Verbosio_et_al_2017_Enhancing_the_scalability_of_selected_inversion_factorization_algorithms_in.pdf;/Users/spech/Zotero/storage/ZNAPM4KU/S1877750317301473.html
ase_field_implementation,Schur-complement decomposition,Selected inversion

Elsevier BV

journaltitleJ. Mech. Phys. Solids
titleA Unified Phase-Field Theory for the Mechanics of Damage and Quasi-Brittle Failure
spech/Google Drive/Zotero/Wu_2017_A_unified_phase-field_theory_for_the_mechanics_of_damage_and_quasi-brittle.pdf

journaltitleComput. Methods Appl. Mech. Eng.
titleRobust Numerical Implementation of Non-Standard Phase-Field Damage Models for Failure in Solids

Elsevier BV

journaltitleJ. Mech. Phys. Solids
titleA Length Scale Insensitive Phase-Field Damage Model for Brittle Fracture
spech/Google Drive/Zotero/Wu_Nguyen_2018_A_length_scale_insensitive_phase-field_damage_model_for_brittle_fracture.pdf

journaltitleSIAM J. Sci. Comput.
titleActive-Set Reduced-Space Methods with Nonlinear Elimination for Two-Phase Flow Problems in Porous Media

Elsevier BV

journaltitleEngineering Fracture Mechanics
titleX-Ray Computed Tomography Images Based Phase-Field Modeling of Mesoscopic Failure in Concrete
spech/Google Drive/Zotero/Yang_et_al_2019_X-ray_computed_tomography_images_based_phase-field_modeling_of_mesoscopic.pdf