Computational modeling of progressive damage and rupture in fibrous biological tissues: application to aortic dissection

This study analyzes the lethal clinical condition of aortic dissections from a numerical point of view. On the basis of previous contributions by Gültekin et al. (Comput Methods Appl Mech Eng 312:542–566, 2016 and 331:23–52, 2018), we apply a holistic geometrical approach to fracture, namely the crack phase-field, which inherits the intrinsic features of gradient damage and variational fracture mechanics. The continuum framework captures anisotropy, is thermodynamically consistent and is based on finite strains. The balance of linear momentum and the crack evolution equation govern the coupled mechanical and phase-field problem. The solution scheme features the robust one-pass operator-splitting algorithm upon temporal and spatial discretizations. Based on experimental data of diseased human thoracic aortic samples, the elastic material parameters are identified followed by a sensitivity analysis of the anisotropic phase-field model. Finally, we simulate an incipient propagation of an aortic dissection within a multi-layered segment of a thoracic aorta that involves a prescribed initial tear. The finite element results demonstrate a severe damage zone around the initial tear and exhibit a rather helical crack pattern, which aligns with the fiber orientation. It is hoped that the current contribution can provide some directions for further investigations of this disease.


Introduction
Aortic dissection is a lethal condition characterized by delamination and separation of adjacent lamellae within the media of the aorta. The annual incidence of aortic dissection ranges from 2 to 9 per 100,000 patients according to Clouse et al. (2004) and Howard et al. (2013). It is of utmost importance to better understand the underlying mechanisms contributing to the development of the disease. Toward this understanding, mathematically robust and physically relevant computational models can identify certain aspects of this intricate phenomenon.

Histology of the aortic wall
The aorta, the main conduit for blood delivery, is an elastic artery which consists of three distinct layers: the tunica intima, the tunica media and the tunica adventitia. The intima is the main layer involved in metabolic processes and is comprised of mono-layered endothelial cells supported by loose connective tissue. The intima provides a negligible mechanical contribution to the wall resistance in young and healthy individuals. Separated from the intima by the internal elastic lamina is the media, which supports the aortic wall against the physiological blood pressure. The media is composed of as many as 70 fenestrated medial lamellar units, hosting two adjacent elastic lamellae (involving elastin), which are interspersed with collagen fibers, smooth muscle cells, glycosaminoglycans (GAGs) and proteoglycans (PGs), as illustrated in Fig. 1a. As reported by Schriefl et al. (2012) the media mainly involves two families of collagen fibers organized successively in the stacked laminar units, each of which contains a single fiber family. The outermost layer of the artery is the adventitia, which prevents the excessive expansion of the wall during systole in a cardiac cycle. In contrast to the media, the adventitia is formed of a network of collagen fibers arranged in a helical structure and is in general stiffer than the media Ross and Pawlina 2011;Humphrey 2013).
A closer look into the structural features of the media and the adventitia shows that the elastin and the collagen fibers are the most significant extra-cellular matrix proteins, as illustrated in Fig. 1a. Out of several damage and fracture mechanisms that occur in a variety of materials, collagen fiber pull-out, collagen fiber bridging, collagen fiber/matrix debonding and matrix cracking seem likely to occur in the aortic wall in mode I, mode II and mixed-mode fracture, as illustrated in Fig. 1b. Such a hypothesis is justified by some studies summarized in Sects. 1.3 and 1.4. That the ultimate rupture stress values for single human collagen fibrils generally reach hundreds of MPa s (see, e.g., Svensson et al. 2013), which is further evidence substantiating our theory that during the rupture of an aortic wall, the fiber-matrix and/or matrix-matrix interactions keeping the wall in register are disrupted, i.e., fibers themselves do not break.

General aspects of aortic dissection
An aortic dissection may appear as an acute or a chronic type, depending on the duration of clinical symptoms (Criado 2011). Of all factors that can be associated with an aortic dissection, hypertension appears to be the most common predisposing factor, diagnosed in 77% of patients (Mussa et al. 2016). Among other factors are the cystic medial necrosis, connective tissue disorders, e.g., Marfan syndrome and Ehlers-Danlos syndrome, aneurysm, trauma, e.g., car accidents, cardiac catheterization, male sex, and individuals aged 60 to 70 years (Dunning et al. 2000;Criado 2011;Tsamis et al. 2013).
In most cases the initial phase of an aortic dissection manifests itself as an initial tear in the intima due to some structural weakening in a localized part of the endothelium, see Fig. 2a, the cause of which is rather elusive. An aortic dissection, however, may also commence within the wall due to an intramural hemorrhage or a hematoma formation in the media (Thubrikar et al. 1999;Khan and Nair 2002). Two-thirds of all initial tears happen within the ascending aorta at about 2 cm above the aortic root near the sinotubular junction (Pasta et al. 2012). The second common location is the isthmus of the descending thoracic aorta located beyond the origin of the left subclavian artery, as depicted in Fig. 2a  In between the elastic lamellae is the embedded collagenous lamella (denoted by the letter C) that involves smooth muscle cells, adhesion molecules as well as GAGs and PGs not shown here due to illustrative reasons. This unit corresponds to the micro-scale outlook of the aortic wall with the typical thickness of the elastic lamella E, about 1.5 μ m, and the collagenous lamella C, about 12 μm , reconstructed from Ross and Pawlina (2011); b possible damage and fracture mechanisms that are likely to occur in the aorta during mode I, mode II and mixed-mode fracture. 1: collagen fiber pull-out; 2: collagen fiber bridging; 3: collagen fiber/matrix debonding; 4: matrix cracking; reconstructed from Anderson (2005) (Cherry and Dake 2009). From the hemodynamics point of view, the most likely reason that makes these areas prone to aortic dissection is the turbulent pulsatile blood flow imposed by the geometry, i.e., the blood pumped out of the left ventricle has to navigate the turn of the aortic arch (Rajagopal et al. 2007;Cherry and Dake 2009). Upon the initial tear formation, which may be seen as a radial propagation across the sub-layers of the media, the dissection changes its course and follows a rather helical path. The dissection separates the lamellar units and creates a false lumen, see Fig. 2a. In any case, a crack propagation across the entire lamellae is not energetically favored (Ottani et al. 2001). A significant amount of the blood now enters the dissected part of the wall, causing even more tearing as blood jets through the tear in each cardiac cycle. Thereafter, the dissection may continue to propagate down toward the abdominal aorta or create an exit (reentry) tear so that the blood can flow back into the true lumen (Cherry and Dake 2009). In the absence of an exit tear, more and more blood enters the false lumen, which decreases the blood volume through the true lumen, leaving the remaining wall susceptible to a dilatation. Over time this may cause other severe pathologies such as the rupture of the remaining aortic wall. Aortic dissections are classified according to the Stanford system and the DeBakey system, as depicted in Fig. 2b. In the Stanford classification, type A involves the ascending segment of the aorta, whereas in type B the ascending aorta is not affected. The DeBakey system, on the other hand, characterizes three different types of aortic dissection: type 1 involves the ascending and descending part of the aorta, and commonly extends beyond the arch distally; in type 2 the only part affected is the ascending aorta; type 3 describes an aortic dissection excluding the ascending part (Tsamis et al. 2013).

Experimental studies on aortic dissection
Early experimental investigations on aortic dissections have mostly followed methods such as pressurization, direct tension and peel tests. Carson and Roach (1990) and Roach and Song (1994) infused fluid via a needle through the media extracted from porcine aortas, causing a bleb formation starting at a non-physiological pressure value of 579 mmHg. Notably, the dissection generated continued to propagate at physiological pressure values. In another study by Tam et al. (1998), the pressure was found to be inversely correlated with the initial tear depth under static conditions. All studies stated above reported on the energy release rates (mJ/cm 2 ), i.e., the energy expended to create a unit area of ruptured material, a significant concept in fracture mechanics. Later the experimental studies focused more on the direct tension and peel tests, which were presented by Sommer et al. (2008), Tong et al. (2011), Pasta et al. (2012 and Wang et al. (2014), among others. All the investigations provided energy release rates, which were found to vary over the different aortic layers and along the different orientations for which the tests were conducted, i.e., circumferential and longitudinal directions.
Microscopical images of the dissected surfaces in the study of Sommer et al. (2008) provided visual evidences for the fiber pull-out and fiber/matrix debonding mechanisms. Nonetheless, direct tension, peel and trouser tests led to a delamination of the adjacent medial lamellae, which is undoubtedly driven by normal stresses ( rr , , zz ), i.e.,  Fig. 2 a Schematic view of an aortic dissection which starts with an initial tear near the left subclavian artery and propagates downwards, while following a path that results in the formation of a false lumen next to the true one; b classifications of the aortic dissection: Stanford type A involves the ascending aorta and Stanford type B excludes the ascending part. In the DeBakey system, type 1 includes the entire aorta, type 2 contains only the ascending aorta, and type 3 excludes the ascending aorta, reconstructed from Tsamis et al. (2013) 1 3 mode I fracture. Later, Sommer et al. (2016) found that ultimate in-plane shear stresses ( u r , u rz ) are about an order of magnitude lower than ultimate out-of-plane shear stresses ( u z ), a result which is based on simple shear tests performed on medial specimens cut from the human thoracic aortas. This result can be attributed to the collagen fibers embedded in the ground matrix that are barely deformed in their mean orientation in the case of in-plane shear tests. Hence, the aortic dissection is less likely to propagate in the radial direction running across several lamellae, a conclusion which was also proposed by Haslach et al. (2018). The authors of Haslach et al. (2018) concluded that the dissection propagation, which can be explained by the relative slip between the adjacent medial lamellae in the circumferential-longitudinal plane, is an in-plane shear driven process evoking mode II fracture.

Numerical studies on aortic dissection
To date, aortic dissection has been the subject matter of several numerical models in solid mechanics. Ferrara and Pandolfi (2010) applied a cohesive zone model (CZM) together with an anisotropic traction-separation law and investigated one of the peel tests conducted by Sommer et al. (2008).
In an experimental and computational study, Leng et al. (2018) employed the CZM to fit the model parameters to load-displacement curves obtained from a shear dominated mixed-mode and a mode I peel test along both circumferential and longitudinal directions. Apart from that, Noble et al. (2017) carried out experimental and computational analyses of a catheter-induced delamination. By using the partition of unity finite element method,  studied peel tests. Other studies, adopting the extended finite element method (XFEM), were presented by Wang et al. (2017Wang et al. ( , 2018 where peel tests similar to the aforementioned contributions were examined numerically. Moreover, the authors therein presented an inflation test of a residually stressed plane strain solid model of a hollow circle representing the cross section of an aortic wall, with varying opening angles. An identical blood pressure was applied on both the inner layer and the tear edge. The tear edge describes a prescribed circumferentially dissected zone where the nodal enrichments are introduced. As a result, the critical pressure value for initiation of aortic dissection propagation was found to increase with the opening angle; therefore, residual stresses seem to protect the artery from tear propagation. Even though several simplified models (mostly 2-D) exist using CZM and XFEM, which are relatively easy to handle in terms of tracking the discontinuities by means of remeshing and nodal/element enrichment functions, more complex and histologically representative 3-D geometries lead to an arduous task when CZM and XFEM are used.
In contrast, the crack phase-field approach (CPFA) by Francfort and Marigo (1998) circumvents the modeling of discontinuities (Miehe et al. 2015a(Miehe et al. , 2016Ambati et al. 2015;Borden et al. 2016). Resembling the gradient damage models, CPFA contains the critical fracture energy g c (critical energy release rate), an essential ingredient of fracture mechanics. Within CPFA, Gültekin et al. (2016) and  studied rupture in human aortic tissues induced by uniaxial extension, simple shear and peeling. They assumed that distinct fracture mechanisms were involved in the overall macroscopic fracture process and they introduced separate critical fracture energies, g iso c for the ground matrix, and g ani c for the fibrous content. These represent the up-scaled, homogenized, macro-resistance of the interactions within the ground matrix and between the matrix and collagen fibers against damage and rupture. A more elaborate account of the extant computational solid models can be found in .
Aside from the benchmark solid models presented in the literature, there are a few patient-specific studies dealing with computational fluid dynamics (CFD) for exploring the underlying hemodynamics of the aorta. Tse et al. (2011) found a blood pressure difference as high as 210 Pa between the true and the false lumen which pinpointed the proximal ascending aorta and the distal aortic arch as the areas subject to the vortex flow, these coinciding with the prevalent initial tear locations, see Sect. 1.2. Besides, on the basis of a helical blood pattern observed solely in the ascending aorta, the authors inferred that the clinically observed helical dissection propagation is the result of the vortex flow. Cheng et al. (2013) showed that the aortic morphology, the initial tear size and the position influence the flow and other hemodynamic parameters. Although most of the CFD studies consider the wall shear stress as a direct contributor to the development of the aortic dissection, such a hypothesis does not seem probable as the wall shear stresses range between 3 and 10 Pa, which are negligibly small compared with the stresses experienced within the wall. It is also worth mentioning some fluid-solid interaction models within the context of aortic dissection (Qiao et al. 2015;Malvindi et al. 2017).

Aim of the present study
In view of our previous findings and contributions (Gültekin et al. 2016;, the present study delivers a computational protocol with novel features for investigating the nascent aortic dissection, and addresses certain mechanical aspects of the protocol, based on the phase-field modeling of progressive damage and rupture.
The article is organized as follows. Section 2 offers a tour of the continuum mechanical and algorithmic framework in terms of geometry, kinematics, governing balance equations and the solution strategy. Section 3 outlines a parameter identification procedure for experimental data obtained from a diseased human aorta and provides a sensitivity analysis of an anisotropy parameter for the extension of a squared single-edge notched domain. In addition, certain physical aspects of an nascent aortic dissection are studied in an idealized cylindrical tube with a prescribed initial tear. Sections 4 and 5 provide a critical discussion and overview of some open problems and possible improvements regarding the modeling concept discussed herein.

Multi-field framework for rupture
This section is devoted to anisotropic phase-field modeling of fracture. The primary field variables, namely the crack phase-field d and the deformation map , are introduced along with their governing equations, i.e., the evolution of the crack phase-field and the balance of linear momentum. Subsequently, an account on the numerical edifice is briefly given which features the operator-splitting algorithm. For the related continuum mechanics see, e.g., the books by Ogden (1997) and .

Primary field variables of the multi-field problem
Let  ⊂ ℝ 3 be a continuum body at time t 0 ∈  ⊂ ℝ + and  ⊂ ℝ 3 at current time t ∈  ⊂ ℝ + in the Euclidean space. The coupled problem of rupture is expressed by the bijective deformation map t and the auxiliary crack phase-field d, i.e., where t maps a material point ∈  onto a spatial point ∈  , see Fig. 3, while d, a thermodynamic measure of damage in the solid, interpolates between the intact ( d = 0 ) and the ruptured ( d = 1 ) state of the material corresponding to a domain regulated by the length scale parameter l, as illustrated in Fig. 4. Note that the crack phase-field d is formulated in the reference configuration .

Kinematics of the mechanical problem
Let the gradient operators ∇(•) and ∇ x (•) denote the gradients with respect to the reference and the spatial coordinates and , respectively. The continuous manifolds are locally equipped with the covariant reference and spatial metric tensors = IJ I ⊗ J and = ij i ⊗ j , respectively, where IJ and ij denote the Kronecker deltas. The bodies  and  admit the deformation gradient and the left Cauchy-Green tensor such that for non-penetrable deformations, i.e., J > 0 , where J = det , see Fig. 3. The energy stored in a hyperelastic isotropic continuum is characterized by the three independent invariants Nonlinear deformation of an anisotropic solid with the reference configuration  ⊂ ℝ 3 and the spatial configuration  ⊂ ℝ 3 . The surface boundary associated with  is defined by  ⊂ ℝ 2 =  ∪  t and  ∩  t = � , and the respective surface boundary in  is given by  ⊂ ℝ 2 =  ∪  t and  ∩  t = � . The surface tractions ̃ and ̃ are applied on  t and  t with unit normals and pointing outward, respectively. The map ∶  ×  →  transforms a material point ∈  onto a spatial point = ( , t) ∈  at time t. The anisotropic micro-structure of is rendered by two families of fibers with unit vectors and ′ . Likewise, the anisotropic micro-structure of is described by and ′ as the spatial counterparts of and ′ , respectively The anisotropic structure of the aortic wall is described by two reference unit vectors and ′ representing the mean fiber orientations, see Fig. 3, with their spatial counterparts, i.e., = and � = � . This idealization of the tissue micro-structure leads to the respective Eulerian form of the structure tensors, i.e., with their Lagrangean counterparts, i.e., = ⊗ and � = � ⊗ � . Next we introduce the (physically meaningful) additional invariants which measure the squares of stretches along each mean fiber direction.

Kinematics of the phase-field problem
Let Γ be a discontinuous boundary such that Γ ∈  ⊂ ℝ 2 at time t 0 , which characterizes a sharp crack surface, i.e., Γ = ∫ Γ dA , as indicated by a thick solid curve in Fig. 4. Instead of tracking such an interface, the phase-field approach approximates the surface integral by a volume integral engendering a regularized crack surface Γ l (d) , i.e., The tensor variable denotes the rotations in the orthogonal group (3) , which contains rotations and reflections. This approximation can be extended to a class of anisotropic materials such that describes the anisotropic crack surface density function subject to the condition (d, ∇d) = (d, ∇d) , ∀ ∈  ⊂ (3) , in which  designates a symmetry group as a subset of (3) . In (7), the second-order anisotropic structure tensor  reads which aligns the crack with the orientation of fibers in the continuum as the phase-field evolves. Therein, the anisotropy parameters M and M ′ regulate the transition from weak to strong anisotropy for two families of fibers. Hence, in the limit case, i.e., i → ∞ , i ∈ {M, M � } , the crack path perfectly lies parallel to the fiber directions. For isotropic solids, the parameters M = M � are zero, whereas for a general anisotropic continuum with several families of fibers, each of them lying in an open range, i.e., −1 < i < ∞ , i ∈ {M, M � , …} , dictated by the ellipticity condition for Γ l (d) (Teichtmeister et al. 2017;).
Schematic view of the diffusive crack topology in  and  with the crack phase-field d ∶  ×  → [0, 1] , while the sharp crack surface Γ smears out in the respective solid domain, denoted by Γ l (d) , which is regularized by the length scale parameter l. The material anisotropy is imparted by the two families of fibers with the mate-rial unit vectors and ′ along with their spatial counterparts and ′ , respectively. The anisotropic crack phase-field problem is further characterized by a traction-free Neumann-type boundary condition ∇d ⋅ = 0 on 

Euler-Lagrange equations of the phase-field problem
In view of (7), we can state the minimization principle for the regularized crack surface Γ l (d) as subject to the Dirichlet-type boundary constraint Upon the minimization of the regularized crack surface functional, we derive the Euler-Lagrange equations, i.e., where the divergence term interpolates d between the intact and the ruptured state of the material, while denotes the unit surface normal pointing outwards in the reference configuration.

A particular form of the degradation function
The macroscopic damage accumulated in the anisotropic solid manifests itself in the mechanical response in the sense of a degradation, which may assume generic distinct functional forms in accordance with the additively split isotropic and anisotropic mechanical contributions, i.e., satisfying the following growth conditions where i ∈ {iso, ani} . The first condition ensures monotonic degradation, while the second and third conditions set the limits for the intact and the ruptured states, and the final condition ensures the saturation of g i (d) as d → 1 provided that a iso and a ani , which control the rate of the mechanical degradation with respect to the evolution of d, lie in the open range (a iso , a ani ) ∈ (1, ∞) . We further emphasize that a iso and a ani can also be functions of stretch or stress and identified via ex vivo biomechanical experiments, see Sect. 4 for a discussion. For the sake of simplicity, we restrict ourselves to the quadratic degradation, i.e., a iso = 2 and a ani = 2 , as given in Miehe et al. (2010b), which also retrieves the multi-field formulation of the fracture problem presented in Gültekin et al. (2016) and .

A particular form of the anisotropic constitutive model
We now briefly describe the specific form of the effective Helmholtz free-energy function for the hyperelastic anisotropic mechanical response of the aortic wall, which can be split into an isotropic and an anisotropic part, i.e., for which the effective isotropic and the anisotropic parts can be expressed as functions of invariants, i.e., The effective isotropic part follows from the neo-Hookean model accounting for the mechanical response of the ground matrix, whereas the effective anisotropic response features the hyperelasticity of the two families of collagen fibers. The explicit forms of Ψ iso 0 and Ψ ani 0 are proposed in  (it is straightforward to adopt here any anisotropic model, see, e.g., Holzapfel and Ogden 2017a, b). In what follows, the mechanical response of the degrading wall due to macroscopic damage in the isotropic and the anisotropic parts is stated via the degradation function in (12) such that which modifies the undamaged energy storage function given in (14). The respective expressions for the Kirchhoff stress and the elasticity tensors are presented in .

Governing equations of the anisotropic fracture
This section is devoted to the governing equations of the coupled multi-field problem of fracture where the classical balance of linear momentum is accompanied by the evolution equation for the crack phase-field, for the strong form of the boundary-value problem. More details can be found in . For the canonically compact gradient-damage formulations of the boundaryvalue problems in regard to the standard dissipative solids, the interested reader is referred to Mielke and Roubíček (2006) and Miehe (2011). ,

Rate-dependent variational formulation based on power balance
As a point of departure, we introduce the viscous rate-type potential Π as The first term  on the right-hand side of (17) represents the rate of energy storage functional, i.e., where the work conjugate variables to and d are the Kirchhoff stress tensor and the scalar energetic force f, respectively, i.e., The second term  on the right-hand side of (17) is a viscous regularized dissipation functional due to fracture, i.e., where the artificial viscosity ≥ 0 regulates the scalar viscous over-stress that reads for which the Macaulay brackets ( < x >= (x + |x|)∕2 ) in (20) filter out the positive values. Note that in (21), g c stands for the critical fracture energy. Finally, the last term  on the right-hand side of (17) denotes the (classical) external power functional acting on the body according to where 0 , ̃ and ̃ represent the material density, the prescribed spatial body force and the spatial surface traction, respectively. Now, with the rate-type potential Π at hand, we propose a mixed variational principle of the evolution problem as with the admissible domains for the primary variables Afterwards, the variation of the potential Π with respect to the fields {̇ ,ḋ, } along with simple algebraic manipulations via elimination and substitution of the respective terms (see  for more details) yields the strong form of the field equations, i.e., The first equation in (25) simply describes the balance of linear momentum, whereas the latter states the evolution equation for the crack phase-field in which  indicates the crack driving source term such that

Energy-based anisotropic failure criterion
Following Gültekin et al. (2016) and , an anisotropic failure criterion is now used. We begin with the assumption that two distinct failure processes govern the cracking of the ground matrix and the fibers whereby the anisotropic structure tensor  in (8) is additively split into distinct forms as We introduce now g iso c and g ani c corresponding to the critical fracture energies attributed to the ground matrix (isotropic) and the fibrous content (anisotropic) of the aortic wall, respectively, which homogenize the mechanical resistance of the respective interactions against rupture. Such a model consideration is suitable for the description of the mechanical response of fibrous biological tissues. The crack driving source term in (26) can therefore be decomposed as For a rate-independent case where → 0 , the above expressions (27) and (28) in conjunction with (25) 2 lead to distinct evolution equations of the crack phase-field in relation to the ground matrix and the fibrous content, i.e., What remains now is to superpose the two distinct failure processes emanating from (29), which leads to the rate-independent evolution equation of the phase-field, i.e., along with the specific form of the dimensionless crack driving source term Relation (31) yields an irreversible and positive crack driving source term such that the maximum positive value of (s) − 1 is taken into account in the entire deformation history s ∈ [0, t] , and the Macaulay brackets filter out the positive values for (s) − 1 keeping the solid intact until the failure surface is reached. Next, we specify the rate-dependent case in view of (30), i.e., where the evolution of the crack is characterized by the balance between the crack driving force and the geometric resistance to the crack (Miehe et al. 2015b).

A note on the weak formulation and numerical implementation
We now give a brief account to the staggered solution procedure of the multi-field problem associated with the primary field variables ( , t) and d( , t) . An identical temporal as well as spatial discretization scheme is employed for the mechanical and the phase-field problem so as to transform the continuous integral equations into sets of discrete algebraic equations. This set of algebraic equations is solved by a one-pass operator-splitting algorithm in a Newton-type iterative solver for the nodal degrees of freedom.
Focusing on a typical time increment = t n+1 − t n in a solution process, where t n+1 and t n stand for the current and previous time steps, respectively, all field variables at time t n are assumed to be known, e.g., n and d n , together with the crack driving source term  n , which is stored as a history variable. In the sequel, the unknown field variables at time t n+1 are sought. Note that for the sake of simplicity, all field variables without a subscript such as and d are hereinafter evaluated at time t n+1 . (31)

One-pass operator-splitting algorithm
The mechanical and crack phase-field sub-problems can be decoupled by means of a one-pass operator-splitting algorithm, i.e., for the time increment . Such an algorithm is extremely robust, although it slightly underestimates the speed of the crack evolution, see Miehe et al. (2010a) and  for details.

Spatial discretization of the mechanical problem
The following notation follows that of Miehe et al. (2010bMiehe et al. ( , 2010a in which a finite element discretization h of the solid  is considered, where h denotes the mesh size composed of E h finite element domains  h e ∈ h and N h global nodal points. In accordance with the discretization h , the finite element interpolations of the deformation map and the deformation gradient can be expressed by the state vector in relation to the nodal position vector ∈ ℝ , where ∈ {1, 2, 3} indicates the space dimension, while serves as a symbolic representation of the global interpolation matrix comprising the shape functions and its derivatives. 1 For a known phase-field d at t n+1 , the algorithmic potential energy functional related to (25) 1 is for which the algorithmic form of the variational principle is described as  The respective Euler equation features nonlinear elasticity at finite strains, which is solved by a Newton-type iteration based on a sequence of updates

Spatial discretization of the phase-field problem
In an analogous way to that in Sect. 2.8.3, we write the finite element interpolations of the phase-field and its gradient by in relation to the nodal phase-field vector d ∈ ℝ , where ∈ {1, 2, 3} , while d serves as a symbolic representation of the global interpolation matrix comprising the shape functions and their derivatives. For a known at time t n+1 , in view of (25) 2 the algorithmic potential energy density functional states Then, the discretized form of the variational principle is written as The respective Euler equation is linear and can therefore be solved in closed form:

Representative numerical examples
Beginning with the material parameter identification for the elastic response, the anisotropic features of the phase-field model are first examined by means of sensitivity analyses, which is followed by the simulation of a nascent aortic dissection within a multi-layered thoracic aortic wall.

Identification of material parameters
We first fit the elastic constitutive response to experimental data obtained using uniaxial and in-plane simple shear tests performed on medial strips, which are cut from aneurysmatic human thoracic aortas (Sommer et al. 2016). In particular, specimens subjected to uniaxial extension are tested in the circumferential -and longitudinal z-directions, referred to as ( ) and (zz) modes, while in-plane simple shear tests are carried out on the radial r plane along the -and z-directions, indicated by ( r ) and (rz) modes, respectively. The elastic parameters are estimated through a nonlinear least-squares analysis, which is based on a single-objective function 2 ( ) characterizing the sum of squares of the analytical model predictions of the Cauchy stresses n (ij) minus the experimentally determined values ̄n (ij) , i.e., The objective function 2 ( ) is minimized with respect to the set of the fitting parameters = { , k 1 , k 2 , } of the constitutive model as stated by . Note that represents the angle between the mean fiber direction and the circumferential -direction. The angle is here conside r e d a s a f i t t i n g p a r a m et e r. F u r t h e r m o r e , = {( ), (zz), (r ), (rz)} denotes the set of the aforementioned modes describing the test (ij) along with the associated number of data points N (ij) exp . A MATLAB ® 2016 built-in function referred to as lsqnonlin is implemented in order to compute the minimization problem. The set of elastic parameters identified is summarized in Table 1 together with the correlation coefficients R 2 (ij) and the root-mean-square error according to which is used as a measure for the 'goodness of fit' (Holzapfel et al. 2005). Therein, q specifies the number of fitting parameters , whereas ̄m ean (ij) is the arithmetic mean of the corresponding Cauchy stresses for each mode. The associated hyperelastic constitutive responses are depicted with a reference to experimental data in Fig. 5, which shows favorable agreement.

Sensitivity analysis of the anisotropy parameter ! M
To demonstrate how sensitive the crack path is with respect to the anisotropy parameter M , a plane strain boundary-value problem is analyzed. In particular, a squared single-edge notched domain is considered with , 38,800 quadrilateral elements, which are connected by 39,295 nodes. The material exhibits anisotropy characterized by a single family of fibers with mean orientation , orientated at an angle = 45 • with respect to the x-axis, see Fig. 6a. While the bottom edge of the domain is fixed in the y-direction ( u y = 0 ), the top edge is subjected to a monotonically increasing vertical displacement ( u y =ū).
The material parameters used for the simulations are = 1.0 kPa , k 1 = 1.0 kPa , k 2 = 1.0 together with the bulk modulus = 3.0 kPa . The phase-field parameters are chosen as g iso c ∕l = 10 −2 kPa and g ani c ∕l = 10 −2 kPa , with the length scale l = 0.1 mm satisfying l > 2h in order to resolve the diffusive crack surface, where h refers to the minimum mesh size. Figure 6b-h depicts the influence of the anisotropy parameter M on the crack pattern, starting with M ≈ −1 up to M = 500 , where the crack starts to follow the orientation of fibers ( = 45 • ), as M increases, which complies with the findings of Teichtmeister et al. (2017) in the small strain context. In fact, the crack path becomes almost parallel to the fibers for M = 500 . To elaborate on the results obtained, let us substitute (8) into (7) 2 which reshapes the anisotropic crack surface density to From (44) it is easily discernible that M serves as a penalty parameter, which enforces ∠(∇d; ) = 0 • as M → ∞ , making the crack path aligned with the mean orientation of the fiber family. It is worth emphasizing that a slight crack kinking is observed for M ≈ −1 which may be a result of the limit imposed by the ellipticity condition and possible multiple minima of the energy encountered on the path. The anisotropy parameter M is responsible for the transition of the fracture mechanism from fiber bridging to matrix cracking as its value increases, see Fig. 1b. Stability becomes an issue for cases where M ≥ 10 , which causes a termination of the computation prematurely for the standard displacement-based Q1 finite element formulation. Such a predicament is not observed within the small strain context (Li et al. 2015;Teichtmeister et al. 2017). The reason may be found in both the finite strain framework employed here, and the exponential stiffening attributed to the fibers that might have hampered the model stability for d → 1. A theoretical prediction for the crack angle is suggested by Takei et al. (2013). These authors translate the maximum energy release rate concept, see Erdoğan and Sih (1963), into a graphical representation by analogy with the so-called Wulff's plot for crystal growth (Herring 1951). The graphical construction consists of a polar plot of the inverse anisotropic critical energy release rate G −1 c ( , ) and a line plot of the anisotropic energy release rate G( ) imposed by loading, which is tangent to the polar plot. This tangency marks the angle with which the crack propagates in the anisotropic continuum. Although this method is applicable within the small strain context, its use for finite deformations is debatable, as the experimental analysis in Takei et al. (2013) neglects elastic stretching. Figure 7a, b shows the corresponding force-displacement curves along with the sensitivity of the crack angle with respect to the anisotropy parameter M . The force required for fracture is remarkably elevated by the increase of M , which can be attributed to the increased effective length scale parameter due to M (Gültekin et al. 2016), thereby resulting in a greater geometric resistance against fracture. Aside from that, a notably sensitive character of the crack angle associated with relatively low values of M is followed by a saturation-type behavior for larger values of M (>100). Note that in order to obtain the traceable curve in Fig. 7b, additional computations with varying M were performed.

Aortic dissection propagation
This example marks a proof of concept in regard to the 3-D modeling of an aortic dissection propagation upon its initiation, which delineates a helical pattern within the multi-layered wall structure, specifically inside a medial sublayer in the neighborhood of the prescribed initial tear due to stress concentrations.

Geometry and material
A segment with H = 40 mm length is isolated from the human ascending aorta possessing typical geometrical values, measured at the end-diastolic phase, with inner and outer aortic radii R i = 15 and R o = 17.5 mm , respectively, as reported by Mao et al. (2008). The geometrical setup is tailored for an idealized geometry, which features a cylindrical tube consisting of 6 layers. Starting from the endothelium, the first four layers (with color codes pink, blue, cyan, green, as illustrated in Fig. 8a) belong to the combination of intima and media, while the outermost two layers represent the adventitia (with color codes yellow and orange). Each of the medial and adventitial sub-layers has reference thicknesses of T med = 0.375 and T adv = 0.5 mm , respectively.
The initial tear size and tear-shape are assumed to be a priori known and span three medial sub-layers across the thickness of the wall, i.e., from the endothelium up to media 4 sub-layer. Notches with varying length R i ∕180 • , where ∈ {30 • , 60 • } , with a width of w = 2 mm , are incorporated into the solid model to examine the influence of the initial tear size on the progression of the dissection, as depicted in Fig. 8b.
The parameters identified in Sect. 3.1 stand for a degenerated media whose constitutive response exhibits a mechanical degradation which is expressed by the related material parameters. In particular, the degenerated media corresponds to the sub-layer media 3. Since respective tests on the healthy medial specimens from the same sample are lacking, we have increased the values of the constitutive parameters , k 1 , k 2 by 20% and attributed them to media 1, media 2 and media 4, as summarized in Table 2. Because of the lack of experimental data of the adventitial layer, we assume a relatively stiffer response of the healthy adventitial layer with respect to the healthy medial layer, and use (rather) arbitrary constitutive parameters.
As for the phase-field parameters, a direct measurement is generally impeded by the size effect during a rupture test. Most of the extant studies solely report the ultimate stress and stretch values that evoke a rather rudimentary information on the tissue strength. Therefore, arbitrary values for g i c ∕l are considered for the sake of proving the concept elucidated in the present study. Nonetheless, the anisotropy parameters are specified in the light of the sensitivity analysis, as described in Sect. 3.2.

Mesh and fiber orientation
The corresponding finite element meshes consist of fournode tetrahedral elements, see Table 3, with a constant length scale parameter l = 0.1875 mm . Figure 9a depicts a typical meshed geometry for the problem considered. In addition, for the sake of simplicity, the fitted angle = 44.71 (see Table 1) between the orientation of the fibers { , � } and the circumferential direction for media 3 is applied to each of the sub-layers of media and adventitia in a discrete sense, as visualized in Fig. 9b, c, respectively. Fig. 8 a Idealized geometry of an extracted 3-D segment obtained from a human ascending aorta composed of four medial sub-layers associated with the color codes pink, blue, cyan and green, and two adventitial sub-layers represented by yellow and orange color. Note that media 3 refers to a degenerated layer (with lower strength); b sliced view of the geometry depicting the prescribed initial tear size depending on the varying parameter that regulates the length of the tear. All dimensions are in millimeters  Table 2 Elastic and crack phase-field parameters related to the individual layers used for the extension-inflation-torsion analysis  Table 3 Total number of nodes and elements pertaining to each geometry given in Fig. 8b designed according to the parameters controlling the length ( ) and width (w) of the initial tear

Boundary and loading conditions
The solid domain is subjected to an extension-inflation-torsion test with appropriate boundary conditions, see Fig. 10a, which is performed in two loading cycles. The first loading cycle refers to a physiological state, while the second cycle refers to a (rather assumed) supra-physiological loading state. In particular, the physiological aortic pressure p ranges between 80 and 120 mmHg, while the supra-physiological state reaches an assumed peak pressure value of 600 mmHg reproducing hypertension, the most common predisposing factor of patients suffering aortic dissection (Mussa et al. 2016). Both of them are applied on the inner surface of the wall in a saw-tooth loading manner, as depicted in Fig. 10b. Arteries are significantly pre-stretched, and the axial deformations are close to zero during pressure cycles (Schulze-Bauer et al. 2003). Following the study by Horný et al. (2014) on the age-dependence of the axial pre-stretch values of human abdominal aortas, a representative axial displacement of û z = 8 mm , which is the equivalent of an axial stretch of z = 1.2 , is applied during the physiological state and maintained during the supra-physiological loading cycle, see Fig. 10c. Experimental measurements suggest an end-systolic twisting angle ̂ for a healthy left ventricle which ranges between 8 and 12 • in the physiological state, as reported by Carreras et al. (2011). These values can also be predicted for the ascending part of the aorta. Nonetheless, a higher value of the twisting angle may occur. We assume here a peak twist angle of 10 • and 30 • with regard to the physiological and supra-physiological scenarios, as illustrated in Fig. 10d.

Simulations and numerical results
In view of the loading scenario described in Sect. 3.3.3, all simulations start with a time increment of = 10 −2 , which is decreased up to = 10 −4 when a stability issue occurs. The 3-D evolution of the crack phase-field d in the degenerated medial sub-layer (media 3) at instants A , B , and C is depicted in Fig. 11. In particular the damage zone of the thoracic aortic segment is shown, where d ≥ 0.8 for varying initial tear sizes, as specified in Fig. 8b. As a matter of fact, none of the geometric descriptions result in an acute/excessive damage zone around the tear at instant A referring to the peak physiological loading state indicated in Fig. 10. The crosssectional view obtained from the top plane at z = 40 mm in Fig. 12 also illustrates the situation at instant A. The analysis shows a crack initiation around the initial tear due to stress concentration, which propagates in a specific manner by aligning with the direction of the first fiber family, as seen in Fig. 11, where d ≥ 0.8 ; the crack follows a rather helical path in the 3-D domain. A special focus is now given to the related stress distributions in Fig. 13a-c (for the case of Fig. 11b). These figures indicate a significant loss of the load-bearing capacity of the degenerated medial sub-layer (media 3) at instant B within the damage zone, compare with Fig. 11b. Such a mechanical degradation is undoubtedly accompanied by a loss of intactness in the respective sub-layer. In reality it causes the blood to enter the wall through the initial tear which will peel off media 1 and media 2 from media 4, adventitia 1 and adventitia 2. A cascade of supra-physiological cycles would trigger even more tearing as the blood jets through the medial sub-layer yielding a false lumen next to the true one. It is also worth highlighting that larger initial tears provide larger damage zones, which are associated with higher stress concentrations, as

Discussion
In the light of the mechanical tests documented by Sommer et al. (2016) and Haslach et al. (2018), the focus is placed on the ubiquitous (elastic) mechanical factors involved in aortic dissection, particularly on the normal and in-plane shear stresses. The elastic material properties are identified from experimental data, as depicted in Fig. 5. However, it is most likely the case that a certain amount of damage, e.g.,  Fig. 10 a Displacements are constrained at the bottom plane where z = 0 mm , along the x-, y-, and z-directions, required to twist the specimen at the top plane at z = 40 mm ; loading conditions for the extension-inflation-torsion test are realized by one physiological and one supra-physiological cycle in a saw-tooth manner with regard to b aortic pressure p ; c axial displacement û z (remains constant after the peak in the physiological cycle is reached); d twisting angle ̂ . Snapshots of the results are shown at instants A , B , and C at time t ∈ {0.4, 1.2, 1.6} corresponding to the peak physiological, supraphysiological loading states, and to the end of the simulation stress softening, is induced prior to the ultimate stresses, which motivates further studies on the parameter quantification. Another important observation we identified during the analyses is that displacement driven tests such as uniaxial extension, shear and peel tests seem to overestimate the rupture properties. In particular, for a dissecting aortic tissue, basically no severely damaged zone is achieved when the energy release rates are used from the literature (Sommer et al. 2008(Sommer et al. , 2016Tong et al. 2014;Leng et al. 2018).
Separation of the medial lamellae seems to follow a twofold mechanism. Inhomogeneity in the respective mechanical properties results in in-plane circumferential and longitudinal shear components ( r , rz ) during inflation of the aortic segment. They are probably responsible for the rupture of interfacial bonds between two adjacent lamellae which is in line with the fracture mechanisms described in Sect. 1.1 (fiber pull-out, fiber bridging, fiber debonding, matrix failure). This rupture enables the blood to enter the interface through the initial tear, while the lamellae have mostly lost their mechanical resistance as d ≥ 0.8 . In a sense, the crack propagation seems to follow a mode II type of fracture rather than mode I. In a nutshell, inhomogeneous in-plane shear deformations catalyzed by the heterogeneous material properties evoke mode II fracture in the form of the failure mechanisms shown in Fig. 1b. That forms the basis for the separation of the medial lamellae leading to aortic dissection as observed at the macro-scale.
A systematic characterization of the elastic and rupture properties of the aorta in health and disease is of particular importance to fracture models in order to cope with the elusive phenomenon of aortic rupture. In fact, constituent-specific mechanical tests (uniaxial extension, shear, peel tests etc.) on adjacent tissue strips extracted from the ascending and descending parts of the aorta should be performed after elastase (breaks down elastin) and collagenase (breaks peptide bonds in collagen). Such enzymatic removals (enzymolyses) of elastin and collagen have been studied by, e.g., Roach and Burton (1957) and more recently by Schriefl et al. (2015). To a certain extent this enables a better understanding of the mechanical role of elastin/collagen and would Fig. 11 3-D evolution of the crack phase-field in the degenerated media 3 at instants A , B and C , according to Fig. 10, in regard to varying tear sizes: a = 30 • , w = 2 ; b = 60 • , w = 2 (see Fig. 8b). Isosurfaces are used to visualize the damage zone corresponding to d ≥ 0.8 allow a refined quantitative assessment of the individual rupture behavior via the identification of the degradation parameters a iso , a ani (see Sect. 2.5) and the critical fracture energies g iso c , g ani c (see Sect. 2.7.2). The determination of the layer-specific elastic and rupture properties of the ascending and descending aortas would definitely enhance our understanding of the role of altered mechanical wall properties, and better inform computational models.
It has been speculated that PGs contribute to the mechanics of arterial walls by linking the individual collagen fibers together. In this respect, matrilysins (Ross and Pawlina 2011) can be used before mechanical tests to better decipher the role of PGs on the mechanical wall response. The reader is also referred to the computational study of Roccabianca et al. (2014), which presents finite element simulations that support a hypothesis for the initiation of local delaminations that subsequently propagate as dissections. In particular, according to the hypothesis of Humphrey (2013), the authors of Roccabianca et al. (2014) show that the pooling of GAGs/PGs in the medial layer of a thoracic aorta can lead to significant stress concentrations through intra-lamellar Donnan swelling pressures that disrupt the normal cell-matrix interactions and delaminate the layered micro-structure of the aortic wall. The pooling of GAGs/PGs may be initiated by an increased transforming growth factor-(TGF-). As a consequence, a significant reduction of the radial elastic properties due to elastic fiber breakage may take place between the elastic laminae before the event of an aortic dissection (MacLean et al. 1999). Such a condition can initiate a local delamination and/or altered mechanosensitive cellular response leading to dysregulated wall homeostasis (Humphrey 2013). Consequently, constitutive models need to be improved by including more realistic mechanisms, as described here, and we need to better identify under which load combination the crack is initiated so that clinical events can better be simulated.
It should also be underlined that the simulations we documented in Sect. 3.3 provide limited validity from a quantitative point of view. Nevertheless, to show the capability of the algorithm as a proof of concept to determine the incipient progression of a dissection is the ultimate goal of this example. The phase-field approach can also be combined with XFEM in order to visually capture the delamination phenomenon and the formation of a false lumen. Thereby, the global problem would probably lead to a dynamical problem which might be a computationally troublesome task to undertake via implicit solvers; therefore, an explicit analysis is recommended.
In our previous study we demonstrated the numerical performance of the crack phase-field model by analyzing uniaxial extension and simple shear fracture tests performed on tissue samples obtained from a human aneurysmatic thoracic aorta (Gültekin et al. 2016). We need to collect more experimental data from human tissues harvested from dissected aortas in order to better validate our model via comparison. Future modeling efforts should also be informed by imaging data extracted from cardiac computed tomography or cardiac magnetic resonance angiography (Baliga et al. 2014;Doyle and Norman 2016). The best would be having access to imaging data from the aorta before and after Fig. 12 Circumferential evolution of the crack phase-field in the degenerated media 3 at instants A , B and C , according to Fig. 10, when viewed from the top plane at z = 40 mm in regard to varying tear sizes: a = 30 • , w = 2 ; b = 60 • , w = 2 (see Fig. 8b). Isosurfaces are used for the cross-sectional view to visualize the damage zone corresponding to d ≥ 0.8 Fig. 13 Distributions of a circumferential Cauchy stress , b in-plane shear stress r and c in-plane shear stress rz with the initial tear size = 60 • , w = 2 (see Fig. 8b) obtained at the instants A , B and C . Transparency is used with slices extracted on planes at z = 10 , 15, 20, 25, 30, then at x = 0 , and y = 0 (for the related coordinate system see a) the event of a dissection from the same patient. Such data should also capture the origin of the crack with a resolution allowing detailed 3D reconstruction. Such geometrical data would also better allow to show the potential of the proposed methodology.

Conclusion
The present contribution addresses the incipient aortic dissection that takes place in a degenerated medial sub-layer in the vicinity of a prescribed initial tear. The aortic wall segment comprises several layers of intima/media and adventitia with various mechanical properties. The proposed crack phase-field model when combined with appropriate hyperelastic constitutive relations captures the anisotropic rupture behavior of the aortic wall. Together with improved imaging techniques and computational power, high-end mechanistic models can further be calibrated and optimized to shed more light on the mechanical aspects of aortic dissections.