Numerical analysis of crack path stability in brittle porous materials

.


Introduction
Engineering porous (cellular) materials, manufactured or natural ones, such as paper, foam, tissue, fabric, honeycombs or bone, have often excellent mechanical and physical properties, besides having a relative low weight.They are frequently subject to complex mechanical loading conditions, e.g. when used as cores in sandwich structures [1], bone-repairing materials [2], packaging solutions [3], or filters [4].To design light-weight materials able to resist catastrophic fracture, a deeper understanding of the fracture mechanisms in cellular materials is needed.With such knowledge, it would be possible to architect the microstructure, e.g. by optimizing pore size distributions in microstructural regions subject to different loading conditions.
Fracture in homogeneous brittle materials can, in some cases, be reasonably well described by traditional linear elastic fracture mechanics (LEFM) [5].According to LEFM, a stationary straight crack subject to a pure opening load in its orthogonal direction will grow along its crack plane, while under a mixed opening/shear load, the crack will grow along a path that maximizes the local opening of the crack while reducing shear movements of the crack lips in the tip region, cf.[6,7].For linear elastic crack growth in homogeneous materials, an assumption of an opening kink in a direction giving a pure opening mode is frequently used, which is equivalent to assuming crack growth in a direction that minimizes shear crack growth, cf.[8][9][10].In strong contrast to homogeneous materials, porous (or cellular) materials present a different set of challenges.Their fracture behaviour can sometimes be predicted by LEFM [11], nevertheless, in many other cases, it cannot [12][13][14][15][16][17][18][19][20].A crack in porous materials can, for instance, initiate at locations far away from a present crack-tip [21].The crack paths in cellular materials may differ from those in homogeneous materials [22].Experimental observations of fracturing porous materials reveal crack kinking also under pure opening loadings, Fig. 1, in strong contrast to what is observed in homogeneous materials.
An important parameter of cellular, or porous, structures is the ratio of global density of the structure, ρ, to that of the material in the solid cell edges and wall, ρ s , i.e., ρ = ρ/ρ s .In this study, ρ (0 < ρ ≤ 1) is referred to as relative density.Thus, relative density is related to porosity, another commonly used term, as 1 − ρ.Some experiments on crack paths in porous materials (mostly rigid foams) subject to mixed mode loading have been reported.In [23,24] polyurethane (PU) foams of varying relative densities (0.085, 0.123 and 0.256) were tested using asymmetric semi-circular bend (ASCB) and single edge-crack (SEC) specimens.The observed crack paths were reported to deviate from LEFM, however, for the foam with the highest relative density examined, the crack paths were found to be closer to the theoretically predicted.Similarly, in [25] PU foams with relative densities of 0.085, 0.11, and 0.17 were tested using compact-tension-shear (CTS), edge-notched disc bend (ENDB) and ASCB specimens.The experiments revealed dissimilar fracture trajectories and a conclusion was that the fracture paths could not be predicted by the commonly used criteria [25].Further, mixedmode fracture in rigid polyvinyl chloride (PVC) foams were analysed in [22] using CTS specimens, and the observed fracture paths were reported to be scattered [22,27,28].Similarly, fracture of PU foams (with relative densities of 0.08, 0.123, and 0.256) under mixed-mode loading conditions were studied using asymmetrical four-point bending (AFPB) and ASCB specimens [30].Also, here it was concluded that the fracture trajectories were scattered and could not be predicted by any fracture criteria.Likewise, in [28] were PU foams with a relative density of nearly 0.3 tested under mixed-mode loading conditions.The measured fracture toughness' were scattered and reported to not follow any criteria.The fracture behaviour in porous materials differs from homogenous materials because their microstructure greatly influences crack initiation and propagation [1,[31][32][33].In a homogeneous material, (nearly) singular stresses occur at a crack-tip and drive the crack growth process.In porous materials, however, the stress field in the vicinity of a crack-tip is altered by the microstructure.At intrinsically weak locations, such as defects, large pores, or bad arrangements of pores, local stresses can be elevated to a similar level, or even higher, than the stress field in the crack-tip region [11].Thus, in porous materials, the crack growth mechanism is seemingly governed by a competition between stresses at weak locations in the microstructure and stresses at crack-tips, if such are present, [12,18,19,21].The pores can actually shield crack-tips and make the material seemingly tougher, cf.[17].Similarly, the crack paths are also influenced by inherent variations in the porous microstructure.Struts between pores can break and link adjacent pores into  several micro-cracks, as depicted in Fig. 2. Such micro-cracks can coalesce to form a larger crack close to, or far from, a present crack and result in stress concentrations which in turn affect the crack path [34][35][36].Hence, the crack path in porous materials may primarily be governed by its microstructure rather than the remote loading [34,37,38].However, as the relative density increases, this trend should be altered, since the porous materials become closer to a homogeneous continuum.
Whether or not LEFM can be applied to cellular materials, and under which conditions, is still an open question.Intuitively, the relative density (or porosity) could be a crucial variable since it indicates how close porous materials are to the solid of which it is made.In the limit when the relative density approaches one (or the porosity approaches zero), the material approach a homogeneous continuum and LEFM assumptions should prevail.On the other hand, when the relative density is low, the material is composed of large empty pores rather than solid material.To understand the influence of the relative density (or porosity) on the fracture process in porous materials, and to which extent a porous microstructure can be considered as a continuum on a global scale, fracture in idealized two-dimensional porous materials of different relative densities subject to mixed-mode global loading is investigated here.

Theory
Consider an elastic brittle porous body Ω, subject to boundary tractions T on part of the boundary ∂Ω T ⊂∂Ω and boundary displacements U on part of the boundary ∂Ω U ⊂∂Ω, Fig. 3.It is assumed that a small deformation elasticity theory is applicable.Let the displacement field be denoted by u.The total potential energy, in the absence of any volume forces, is then given by [39,40], where ψ e is the elastic strain energy density and ψ s is the surface energy density, represented as a regularized monotonic damage field d ∈ [0, 1] defined over the entire body, where d = 0 corresponds to undamaged material and d = 1 to completely broken material.A crack density functional is introduced, often referred to as AT1 (Ambrosio-Tortorelli 1), cf.[41][42][43], which is a linear approximation of the Munford-Shah functional [53] according to where l is a regularization parameter, often considered as an internal characteristic length, that governs the crack topology and controls the width of a regularized crack, G c is the critical energy release rate and ∇ is the spatial differential operator.The crack profile, d(ξ), acting along ξ, is obtained from the variational minimization of Eqn.Moreover, following e.g.[39], the strain energy density is split into an active/tensile component ψ + e , which is degraded and an inactive/compressive component ψ − e , which does not degrade.Assuming a degradation function of the form (1 − d) 2 , cf. [44,45], The split of the strain energy density into two parts aims to ensure that the material do not degrade under compression and stiffness is kept in the case of crack closure.However, the choice of the split is not obvious, and many different variants have been suggested to overcome diverse unphysical behaviour in compression and shear [46].Here, a hydrostatic-deviatoric split, cf.[47], is applied in which the hydrostatic and deviatoric deformation modes are the independent modes of deformation which are split into the positive and negative contributions according to where ∊ = 1 2 [∇u + (∇u) T ] is the linearized strain tensor, e the linearized strain deviator, the parameter α = 1 if the volume expands, i.e. tr(∊) > 0, and α = 0 otherwise.The bulk and shear modulus are given by K and μ, respectively.One may notice that this split assumes that all strain energy associated with shape deformation affects damage, also under compressive load.Moreover, to impose irreversibility of the crack evolution in computer simulations, a history field H is used in place of ψ + e , such that H is the maximum positive strain energy density experienced [44,46], i.e., where t denotes time.The history field (5) ensures that it is the largest strain energy density experienced in the material during its loading history that determines the present damage state, i.e. the material cannot heal.Now, using the above relations, the total energy functional of the system is written as, cf.[45], The first variation of ( 6) should be zero according to the principle of least action.The Euler-Lagrange equations of ( 6) with respect to the displacement and damage fields are where the Cauchy stress tensor and the consistent material stiffness tensor are given by σ = ∂ψ e /∂∊ = (1 − d) 2 ∂ψ + e /∂∊ +∂ψ − e /∂∊ and C = ∂ 2 ψ e /∂∊ 2 , respectively.See e.g.[52] for a comprehensive review of the topic.

Model
Consider a cellular brittle elastic body with a quadratic cross-section with edge length L containing a stationary edge crack of length L/2 positioned at the cross section's mid-height, Fig. 5.A Cartesian (x 1 , x 2 , x 3 ) coordinate system is introduced with its origin coinciding with the cross-section's center and the crack-tip.It is assumed that the body has a uniform topology in the x 3 -direction, consequently the model is restricted to open porous structures.The extent of the body in the x 3 -direction is further assumed to be large, thus plane strain conditions prevail.Moreover, the body is loaded so that both shear and tensile stresses may be present.The boundaries x 1 = ±L/2 are held stress free while the boundary conditions acting at the edges x 2 = ±L/2 are given as prescribed displacements, i.e., where Δ > 0 and the load parameter β ∈ [0, 1] determines the ratio of global opening load vs global shearing load.When β = 0 the load is pure opening while when β = 1 the load is pure shear.When 0 < β < 1, the load is a combination of both opening and shearing, i.e. a mixed-mode loading.
The pores (or cells) are randomly distributed in the cross-section domain, i.e., there is no preferred orientation in any direction, and the body is considered transversely isotropic on the global scale.For a certain relative density ρ, the cells in the body have approximately the same size.The in-plane relations between global and microscopic (i.e., walls/edges) material properties follow the approximations [11,15,31], where the global Young's modulus Ê, the global shear modulus μ and the global critical energy rate Ĝc are scaled counterparts to the microscopic (local) properties E, μ and G c .In the limit when the relative density ρ → 1, the body becomes homogeneous (and the cells vanish).

Numerical method
The equations are solved using well-established finite element algorithms implemented in a Matlab [48] code.The cellular bodies are discretized with standard bi-linear 4-node isoparametric elements assuming plane strain conditions.An in-house mesh algorithm generates a two-dimensional mesh of randomly positioned polygon pores using a Voronoi tessellation technique.The tessellation process is controlled to generate uniformly distributed cells with small variance in size, which means the study is valid for materials with nearly evenly distributed cells of similar size and with a uniform topology in the out-of-plane direction.
The initial edge cracks are formed by simply removing all the finite elements crossing the crack line x 1 < 0 and x 2 = 0.The meshes always have a cell located with its centre at the origin, i.e., the initial crack-tip is placed inside a pore/cell to allow comparison between different discretized structures.Generally, cracks can propagate in cellular materials either by separation of the pore/cell walls, like in wood [15,50] due to the special lamella structure of the cells, or by breaking the cell walls [11][12][13][14]21], like in rigid foams, or a combination of the two mechanisms.The transition between cell wall separation and cell wall breaking depends on the architecture of the microstructure.However, in this study, the fracture toughness is assumed equal everywhere in the porous solid, meaning no preferred fracture mechanism is selected and the fracture positions are solely determined by the potential energy minimization problem (6).
The number of pores in the discretized bodies is about 20 in each direction x 1 and x 2 , Fig. 6.Thus, there are sufficient pores between the initial crack-tips and the domain boundaries to reduce boundary effects.The number of cells in each direction was selected to give consistent results while aiming to keep the computational cost at a reasonable level, based on a convergence study provided in the Appendix.All initial cracks in the cellular bodies cut about 10 cells which are sufficient to localize initial damage to the crack-tip region, cf.[12,15,21].
The final meshes, resembling cellular microstructures similar to e.g.rigid open foams, some softwoods or parts of trabecular bone, have relative densities ranging from ρ = 0.1 to ρ = 1.Fig. 6a shows a small part of a structural mesh, while Fig. 6b and 6c display examples of complete meshes with small and large relative densities, respectively.The number of elements in each unique mesh (body) varies depending on its specific relative density and the random generator seed used in the Voronoi tessellation procedure.The mesh in Fig. 6b contains about 56, 000 elements while the mesh in Fig. 6c has roughly 47, 000 elements.The meshes shown in Fig. 6 were made using the same random seed in the Voronoi tessellation procedure.The number of elements along a wall's thickness direction is always at least four to reasonably well capture bending moments of the slender walls, Fig. 6a.The element sizes varies within each mesh due to its irregular topology.However, the largest element edge-length h in each structural mesh is always less than half of the length parameter l, i.e., h < l/2.An element edge-length h < l/2 is crucial to obtain a correct resolution of the regularized crack surface, cf.[42,44,52].A convergence study done with different relative densities and length ratios l/h is provided in the Appendix and it is clear the selected element size is sufficient to obtain consistent results.While the length l is constant in all models to allow for a fair comparison, it is much smaller than the average pore size ̅̅ ̅ c √ in each unique mesh, l < ̅̅ ̅ c √ /4 (c is the average pore area in each mesh).All bodies are loaded with prescribed displacements according to Eqn. (8) in small incremental steps.The coupled displacementdamage field problem (8) can be solved either by a so-called monolithic scheme [37] or a so-called staggered scheme [49].Since the former is known to be rather numerical unstable, cf.[49], the displacement and phase-field equations are solved incrementally in a staggered manner, using Newton-Raphson iterations [42,46].A detailed flowchart of the staggered scheme is given in Fig. 7.
The governing quasi-static equations are, where t K i,j is the finite element stiffness matrix at load step t, corresponding to the present stress and damage states, t R is the load vector representing the driving forces including the body and surface forces at step t, t F i− 1 is the nodal force vector corresponding to the displacement vector t U i− 1 and damage vector t D j− 1 .The superscripts i and j denote the iteration number and the vectors ΔU i and ΔD j represent the increments in the nodal displacements and nodal damages at iterations i and j.
In all models are the following in-plane global material and geometric properties applied: Ê = 210 GPa, μ = 80.77GPa, Ĝc = 2.7 kJ/m, L = 1 mm and l = 12 µm.The properties are selected to represent a typical engineering material, cf.[52].The study aims to compare crack paths at different loading conditions and relative densities.It is underlined that as long as the material properties are kept consistent in the models, the selected values will not influence the results.The bodies are subject to increasing displacement increments of Δ/L = 10 − 4 .When a body is sufficiently loaded, the initial crack is anticipated to propagate along different paths depending on the loading parameter β and the unique microstructure.

Results and discussion
The two special situations when a cellular body is loaded when β = 0 or β = 1 were first investigated since these are two wellknown examples found in numerous studies concerning phase-field fracture modelling, cf.[40,41].Fig. 8 presents typical crack paths obtained in bodies with different relative densities ρ subject to pure global opening load (β = 0).In contrast, Fig. 9 shows typical crack paths in the bodies when they are subject to pure global shearing load (β = 1).For visibility, the crack paths are illustrated by The results shown in Figs. 8 and 9 reveal that computed crack paths/trajectories (i.e., resulting damage fields) in cellular bodies can be different to those obtained in homogeneous bodies having the same global material properties and boundary conditions, depending on their relative densities.The crack paths in the homogeneous models (ρ = 1) under opening (β = 0) and shearing (β = 1) loads are consistent with common criteria in LEFM, e.g.assuming crack opening (kinking) along directions of maximum tensile normal stress, cf.[5,29].However, while the crack paths in the high relative density body (ρ = 0.7) essentially capture the crack paths in the homogenous body, the crack paths in the low relative density body (ρ = 0.1) are very different from the other two cases.Notably, for this unique microstructure and inherent material properties, the computed crack path in the low relative density body subject to global opening load is more similar to the crack path in the homogeneous body subject to a shearing load than to the crack path in the homogeneous body subject to opening load.
To get a deeper insight into the observable differences of the crack paths in porous and homogeneous bodies, bodies with different relative densities subject to varying global mixed-mode loading parameters β are displayed in Fig. 10.One may observe, in Fig. 10, that the crack paths obtained in the bodies with lower relative densities (i.e., at higher porosities) are seemingly stochastic while at higher relative densities the crack paths are quite similar to those obtained in homogeneous bodies.This phenomenon was observed irrespective of the specific random generator seed used when discretizing the meshes in the Voronoi tessellation algorithm.It is emphasized that crack paths in bodies with a specific relative density were not exactly the same because of the random nature (different seeds used in the random generator give different microstructural architecture).For the same reason, the onset of damage growth in bodies having the same relative density are not obtained at exactly the same global load.Nevertheless, the results in Fig. 10 clearly illustrate the instability of the preferred crack paths in porous bodies.
The different crack paths observed in Figs.8-10 are explained by the fact that the deformation fields in cellular bodies are different compared to those in homogeneous bodies because of their open, or sparse, microstructures, which modify the distribution of internal forces.The lower relative density of the microstructure, the larger is the deviation from a classical homogeneous continuum because of the stochastic nature.In a cellular body, the walls/edges alter the strain field in the microstructure and the cells (or pores) limit the strain and stress fields resulting in a mechanical smoothening effect, cf.[17,18].Judging from Figs. 8-10, the selected crack paths in the lower relative density bodies seem to have no preferential direction, the crack trajectories seem independent of loading conditions.This observation suggests that at low relative densities, the crack paths seem to be governed by local microstructural features rather than the remote global loading.
Furthermore, as shown in Fig. 10, at relative densities around 0.5, there seems to be a shift in the fracture behaviour.At relative densities above 0.5, the crack paths are similar to those in a homogeneous body (even though they are diffuse and scattered, because of the cells), while at lower relative densities the crack paths seem more arbitrary (i.e. the random architecture of the microstructure has higher influence).To further illustrate the phenomenon and try to overcome the microstructural randomness, Fig. 11 shows accumulated fracture positions in 10 unique bodies at each relative density and load parameter β.Each dot in Fig. 11 represents the position of a fully damaged node.The more spread the spatial positions of the dots are, the more spread and diffuse are the crack paths for the specific load parameter β and relative density ρ.Thus, the crack paths are more geometrically stable where the dots are less spread, or more focused, less geometrically stable, where the dots are widely spread.One may further observe that the crack paths somewhat converge to those obtained in the homogeneous body (ρ = 1) when the relative densities increase.The crack paths are more stochastic, or geometrically unstable, at relative densities below about 0.5, and the lower the relative density is, the sparser the material is, the less geometrically stable are the crack trajectories.
Thus, the fracture patterns in Fig. 11 demonstrate the influence of the relative density (or porosity) on the crack path stability in transversely isotropic brittle porous materials.At low relative densities, the stochastic microstructures have larger influence on the preferred crack paths than the global remote load.At higher relative densities, in more homogeneous materials, the remote global load governs the preferred crack path while the influence from the microstructure is limited.This observation explains the crack kinking phenomenon visualized in the experiment illustrated in Fig. 1.A low-density rigid foam subject to global opening mode load is expected to have a large variance of preferred directions, as illustrated by the case β = 0 and ρ = 0.1 in Fig. 11, and crack growth along other directions than along the original crack plane straight ahead is not surprising.
Moreover, the findings presented in Figs.8-11 are consistent with earlier experimental and numerical studies, cf.[12,15,21].A sufficient crack length to obtain fracture localization was found to be about four times the average pore size [12,21].A similar phenomenon was also observed in [26] where continuum and high-resolution microstructural phase-field models were compared.
Those results indicate that for microstructures with low relative densities (or large pore sizes compared to the length l), the microstructures will influence the crack behaviour to such an extent that traditional continuum modelling is inadvisable: the microstructures have to be taken into account.Other experimental studies on fracture in brittle cellular materials subject to mixed-mode loading conditions report that both the fracture toughness and the fracture paths deviate from what is expected from commonly used fracture criteria, cf.[25,27,30,55].The relative densities in those studies are reported to be up to 0.3, meaning the experimental observations are well in agreement with what is found in this study.
It is worth noting that the mechanical properties of open cellular materials may not necessarily scale with the relative density as given by the scaling law in Eq. ( 9).For instance, the Young's modulus can possibly scale linearly, quadratically, or cubically with the relative density, cf.[50].The scaling effect was investigated, and some force-displacement curves obtained in the simulations are shown in Fig. 12 for the cases β = 0 and β = 1 when using the same random generator seed in the mesh generation process but different relative densities.As is observed, the global stiffness is reasonable well captured, having in mind the very large difference in stiffness between a body with a relative density of ρ = 0.2 and a homogeneous body (ρ = 1).According to Eqn. (9), the global stiffness of a body with relative density ρ = 0.2 is estimated to be 25 times lower than the global stiffness of a homogeneous body when having equal microscopic (local) Young's modulus.Nonetheless, the small differences observed in Fig. 12 could be explained by mainly-two phenomena: 1) there is a strong biaxiality (Poisson's) effect which makes the homogeneous body, and the bodies with higher relative densities, stiffer, and 2) the initial crack-tip is centred inside a cell meaning that for especially low-density bodies the cell size is quite large, so the load carrying area in front of the initial crack is dramatically decreased (the cell size scale with the relative density as ̅̅̅̅̅̅̅̅̅̅̅ One would be able to compensate for both these effects, however, it would also decrease the transparency in the results.Nevertheless, as is shown in Fig. 12, the fracture strengths (peak-loads) are fairly consistent.Also, it is noted that the force--displacement curves, such as those shown in Fig. 12, might drop vertically at post-peak load.These vertical drops indicate incremental load steps in which significant damage growth take place in load-carrying regions, leading to substantial global stiffness loss.This is a typical feature in phase-field fracture models solved with a fully staggered scheme together with the linear Ambrosio-Tortorelli approximation of the Munford-Shah functional, cf.[41][42][43][51][52][53].One might interpret this phenomenon as unstable crack growth, even though the term often is associated with dynamic phenomena.
The findings presented here are believed to be valid for all transversely isotropic brittle cellular materials which have a uniform topology in the out-of-plane direction (e.g., rigid open foams) and a random distribution of in-plane cell positions.The model does not distinguish between different failure modes in the microstructure (e.g.splitting or breaking of walls) since the fracture toughness is assumed to be equal everywhere in the microstructure, meaning no preferred fracture mechanism is selected, the fracture positions are solely determined by the potential energy minimization problem (6).However, the model can quite easily be modified, following e.g.[50], to account other failure mechanisms such as delamination in cell walls.

Conclusions
A numerical phase-field fracture model is used to analyse crack path stability phenomena in porous solids with randomly distributed cells.Cellular bodies with different relative densities (or, porosities) are subject to varying global mixed-mode loading.The bodies contain initially straight edge-cracks.When the bodies are loaded, cracks grow and the final crack paths are contrasted.
In bodies with low relative density (high porosity), the crack paths appear to be stochastic and governed by microstructural features.An initial straight crack can kink even when the body is subject to pure global opening mode load, while the growing crack would be expected to follow its original crack plane as its homogeneous counterpart.The sparse microstructure influences the distribution of internal forces to a great extent, i.e., the cell walls/edges alter the deformation field and high stressed positions may appear at remote distances from an existing crack.Thus, the crack paths in low relative density bodies are governed by weak links in the stochastic cellular microstructure rather than the remote global load.
In strong contrast, in cellular bodies with higher relative densities -in more homogeneous materials-the deformation fields are more continuum-like and the crack paths are predominantly controlled by the remote global load and the influence from the stochastic microstructure is limited.Cracks growing in high relative density cellular bodies are thus expected to follow similar paths as those obtained in homogeneous bodies (continua).There is a transition of the observed fracture behaviour at relative densities around 50%.At lower relative densities than 50%, the stochastic nature of the microstructure seems to govern the preferred paths for growing cracks.At higher relative densities, the global load seems to govern the fracture paths rather than the microstructure.Hence, when the relative density of porous materials is higher than roughly 50%, the crack paths tend to be geometrically stable and follow those in continua, which can be used as a "rule of thumb" for engineering purposes.
The findings in this study can be used to e.g.optimize the performance of manufactured cellular materials by controlling, or manipulating, their microstructures.For example, porous bone implants can be designed with varying densities, to have higher density in load-carrying regions and lower density in other regions to obtain higher permeability.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.9).The internal length was l = 2h, where h is the maximum element edge length in the discretized structures.

Pore density in the problem domain
Fig. A2 shows crack paths computed in quadratic structures with N c = 10, 20 and 40 pores per unit length L at relative densities ρ = 0.1 and 0.7, subject to global load β = 0 (global opening load) and β = 1 (global shearing load).The porous bodies had initial straight edge cracks with their tips positioned at the center of the problem domains, Fig. 5.The boundary conditions and material parameters are as described in Section 3. The internal length was l = 2h, where h is the maximum element edge length in respective structure.
Figs. A1 and A2 reveal that a mesh pore density of N c = 20 pores per unit length L seems to give consistent results, both in terms of global stiffness and crack paths.When N c = 20, the global stiffness of the bodies is similar to those in bodies having higher pore density.The computed crack paths are similar to those obtained in structures having more pores (N c = 40) as well as less pores (N c = 10).

Internal length l
Fig. A3 shows crack paths computed in structures with N c = 20 pores per unit length L at relative densities ρ = 0.1 and 0.7, subject to global shearing load β = 1.The structures had initial straight edge cracks with their tips positioned at the center of the problem domains, Fig. 5.The boundary conditions and material parameters are as described in Section 3. The internal length was l = h, 2h and 3h, where h is the maximum element edge length in the respective body.One may notice that the width of the cracks increases with increasing l while the crack paths are similar for the three different lengths.
Finally, Fig. A4 shows global stress-strain relations from the simulations resulting in the crack paths in Fig. A3.The results are consistent.The fracture (peak) global stress decreases with increasing l (as they should using the linear Ambrosio-Tortorelli crack functional) while the crack paths are similar for the three different lengths.

Fig. 2 .
Fig. 2. Micro-cracks initiate at weak points in the porous structure.A possible crack path is indicated as a dashed line due to clustered microcracks.

Fig. 3 .
Fig. 3. Representation of a crack.A sharp crack (left) and a diffuse crack (right).

Fig. 4 .
Fig.4.Profile of the diffusive phase-field crack as an approximation to a sharp crack.Here, ξ is a spatial coordinate axis directed orthogonal to the crack plane.

Fig. 6 .
Fig. 6.Examples of cellular meshes with different relative densities ρ.A close-up of a small region in a mesh with ρ = 0.1 (a).Complete meshes with ρ = 0.1 (b) and ρ = 0.5 (c).The same seed was used for the random generator when creating the meshes (b) and (c).

Fig. 10 .
Fig. 10.Typical crack paths obtained in a range of bodies with different relative densities subject to varying loads.

Fig. 11 .
Fig. 11.Accumulated fracture positions in 10 unique bodies at each relative density and load.

Fig. A1 .
Fig. A1.Global stiffness at varying number of pores per unit length L and relative density ρ.

Fig. A2 .
Fig. A2.Computed crack paths in structures having a different number of pores per unit length, N c , and relative density ρ, subject to global opening load (β = 0) and global shearing load (β = 1).

Fig. A3 .
Fig. A3.Computed crack paths in bodies with internal length l = h, 2h and 3h, relative density ρ = 0.1 and 0.7, and subject to global shearing load (β = 1).The number of pores per length L was N c = 20.

Fig. A1 shows
Fig.A1shows computed global Young's modulus E * estimated from uniaxially loaded quadratic domains of side-length L having approximately N c number of pores per unit length L. The porous bodies have relative densities ρ = 0.1, 0.4 and 0.7.The calculated values are normalized with the Young's modulus of a homogeneous continuum Ê.The local in-plane Young's modulus E was obtained using the scaling relation, Eqn.(9).The internal length was l = 2h, where h is the maximum element edge length in the discretized structures.