# Fracture Toughness of Periodic Beam Lattices

Gergely Molnár, Julien Réthoré

► To cite this version:

Gergely Molnár, Julien Réthoré. Fracture Toughness of Periodic Beam Lattices. 2024. hal-04793587

HAL Id: hal-04793587

https://hal.science/hal-04793587v1

Preprint submitted on 20 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Fracture Toughness of Periodic Beam Lattices

© Gergely MOLNÁR<fnref n="1" /> and Julien RÉTHORÉ<fnref n="2" />

The study tackles the challenge of accurately modeling fracture behavior in beam lattices, which is essential for designing robust architected materials. Our research focuses on evaluating how the lattice's microstructure and material properties affect fracture toughness. We employed finite element simulations based on the Euler-Bernoulli beam theory to investigate crack propagation, using a failure criterion that initiates beam fracture when maximum axial stress exceeds critical strength. Building on observations from these simulations, we developed a multi-phase-field fracture model with Cosserat elasticity to integrate consistent toughness characteristics into a comprehensive framework for lattice design. This model was validated through experimental tests, ensuring a close match between theoretical predictions and physical reality.

Our findings reveal that the energy release rate remains relatively stable during crack propagation, underscoring its reliability as a measure of the toughness of periodic lattice structures. We discovered that toughness is predominantly influenced by beam height and material properties such as tensile strength and Young's modulus, while slenderness has minimal impact. Additionally, cracks were observed to preferentially propagate along the lattice's structural directions due to stress localization effects, highlighting the importance of the microstructure in fracture behavior.

The implications of this research are significant, suggesting that improved modeling of fracture in lattice structures can enhance material design reliability and optimization. This study bridges the gap between theoretical models and real-world applications, providing valuable insights for the development of advanced materials with tailored fracture properties.

Keywords Cosserat elasticity, beam lattices, fracture toughness, homogenization, metamaterial

## 1 Introduction

Recent advancements in additive manufacturing have enabled the creation of metamaterials with tailored microstructures (Askari et al. 2020). This small-scale fabrication process not only enhances the stiffness of these materials but also strengthens their load-bearing capacity. Additionally, as the size of the structural elements approach the microscopic length scale of the material, the impact of manufacturing defects diminish. Consequently, materials that are macroscopically brittle can exhibit ductile behavior (Ritchie et al. 2009; Zheng et al. 2016; Bauer et al. 2016) at this small scale, leading to the development of materials that are not only lightweight and stiff but also incredibly resistant.

This phenomenon mirrors what is observed in nature (Nepal et al. 2023) where structural hierarchies are evident in the microstructure of various load-bearing components. Examples include cork (Chen et al. 2010), several diatom species (Jang et al. 2013), honeycombs (Mousanezhad et al. 2015), and trabecular bone (Lakes 1993; Ritchie 2011). In each case, smaller beam-like elements form a complex network designed to withstand specific loads.

Given how rare it is to come up with new bulk materials, both scientific and the industrial communities are turning to architected or composite materials to solve future problems. In this regard, architected materials emerge as a promising avenue to address this pressing issue. They not only serve as structural components but also offer additional functionalities such as sensing, adaptation, selfrepair, morphing, and restoration. This versatile nature (Dimas et al. 2013; Libonati et al. 2016; Nepal et al. 2023) makes them applicable across various industrial sectors, including energy storage, pharmaceuticals, wearable electronics, and human-computer

<footnote n="1">¹ Univ Lyon, CNRS, INSA Lyon, LaMCoS, UMR5259, 69621 Villeurbanne, France</footnote>

<footnote n="2">² Nantes Université, Ecole Centrale de Nantes, CNRS, GeM UMR 6183, Nantes, France</footnote>

<footnote>* gergely.molnar@insa-lyon.fr</footnote>
interfaces. Moreover, they are crucial elements in key sectors like automotive, defense and aerospace industries.

Mechanical metamaterials represent a transformative class of materials structured to exhibit properties not found in naturally occurring bulk materials (Jiao et al. 2023). These synthetic materials gain their unique mechanical properties from their designed geometry rather than their composition. While these new materials offer significant advantages, the complexity of their mechanical behavior, particularly under failure, must be thoroughly understood. Without this understanding, current design processes cannot be effectively applied to these advanced materials.

Modeling the elastic properties of beam structures using classical continuum mechanics often fails to capture the size-dependent behavior observed in lattice structures. The Cosserat medium (Cosserat et al. 1909), a type of generalized continuum model, incorporates additional rotational degrees of freedom and couple stresses, which are essential for accurately describing the mechanics of beam lattices and other architected materials. By employing the Cosserat medium, we can better predict the behavior of beam structures, particularly when dealing with fine microstructures where classical theories may fall short(Molnár and Blal 2023). Unfortunately, there are only few lattice types, where analytic formulae are available to calculated the homogenized Cosserat stiffness (Pradel et al. 1998; Sab et al. 2009). Recently, a study showed how to calculate local Cosserat properties on non periodic 2D lattices (Liebenstein et al. 2018a; Liebenstein et al. 2018b). However, until now no unified algorithm is proposed to obtain the Cosserat stiffness of a given type of architecture.

Even if works in the past few decades have concentrated on predicting the strength of homogenized microstructured materials (Pichler et al. 2011; Kolpakov et al. 2020), they struggle to predict the materials resistance in presents of a crack. Due to the theoretical stress singularity in the homogenized model at the crack tip, scale separation cannot be assumed, as the gradient of the stress peak becomes comparable to the scale of the microstructure. While the theoretical concept of homogenized fracture toughness was initially introduced by Roux et al. (2008); Vasoya et al. (2016); Lebihain et al. (2021), it has thus far been demonstrated only for Gaussian random microstructure model materials. We asked ourselves, if the unit cell of a periodic beam lattice contains all the information, how can the materials toughness be calculated? As brittle fracture is a unique problem characterized by its dissipative nature, stress localization at the crack tip, and the coupling with the microstructure, it is still a topical debate if intrinsic fracture properties can be defined independent of the macroscopic structural context.

This theory lends itself well to implementation into the phase-field technique (Bourdin et al. 2000; Molnár and Gravouil 2017), which has demonstrated capability in simulating crack propagation in inhomogenous materials (Bharali et al. 2021). Furthermore, the phase-field approach can be enhanced to account for anisotropy through degradation of elastic stiffness (soft) and fracture energy using penalty functions (hard). The hard representation (Nguyen et al. 2017) requires a penalty function which acts on the gradient of the damage in the phase field formulations. While in the soft model (Bleyer et al. 2018; Scherer et al. 2022) introduces a directional dependent $g_c^M$. The soft model is preferred in our case as it accurately describes the behavior of the unitary cell structure independent of crack length and geometry. While promising, this soft representation has not yet been calibrated to a real cell structures.

In this study, we employed the Euler-Bernoulli beam theory to model crack propagation within lattice structures, focusing on identifying the parameters that influence the critical energy release rate during fracture. By testing the impact of various factors on crack propagation, we were able to determine how the lattice's microstructural and material properties affect its toughness. To extend our findings, we proposed a unified algorithm for calculating Cosserat constants applicable to arbitrary periodic beam lattices, which enhances the accuracy of macroscopic elasticity predictions. Additionally, we utilized a multi-phase-field approach, calibrated using real microstructures, to model the fracture behavior of these lattices. Our theoretical and computational work was complemented by experimental verification, specifically on rectangular beam lattices, to ensure the validity and applicability of our models in practical scenarios.

The paper is structured as follows. First, Section 2 states the elementary question, then Section 3 presents the modeling techniques used to tackle the problem. This is followed by
a comprehensive analysis of the energy release rate in the beam model in Section 4. After presenting the phase-field results in Section 5 and the experimental validation in Section 6, the paper concludes with a final discussion in Section 8.

## 2 Problem statement

This section introduces the mechanical problem under investigation and outlines the objectives of the analysis. Specifically, it compares two approaches to modeling fracture in architected materials: the Euler-Bernoulli beam theory and the phase-field fracture technique.

Given the well-established nature of these methods, this section provides a concise overview of their fundamental principles, defines the necessary input variables, and outlines potential outputs. The primary focus is on comparing how these methods model fracture in beam lattices, a common type of architected material.

Fracture resistance is characterized differently by each approach. In the beam model, the microstructure is explicitly represented, with individual beams modeled according to their geometry and material properties. When a beam fractures, it can be discretely identified. In contrast, the phase-field technique treats the material as a continuum, using homogenized material parameters to describe its behavior. Here, cracks emerge gradually as a damage variable diffuses through the material.

For the beam model, the microstructure is defined by parameters such as the height ($h^m$) and lengths ($L^m$) of the elementary beams. These structures are characterized by their material properties, which describe stiffness (Young's modulus, $E^m$) and strength (the maximum tensile stress they can withstand, $\sigma_c^m$).

On the other hand, the phase-field approach is a diffuse damage model that represents the elastic behavior of the material using homogenized continuum stiffness, described here by the Cosserat theory. This theory is essential for capturing the complex behavior of beam lattices, as demonstrated previously (Molnár and Blal 2023). Fracture initiation and propagation in this model are governed by the critical energy release rate ($g_c^M$, referred to as toughness) and a regularization length scale ($l_c^M$), which controls the extent of damage diffusion.

This analysis seeks to answer three key questions: (i) Can we define a unique fracture toughness for a given beam lattice? If so, what parameters influence this toughness? (ii) Can we develop a phase-field continuum model that accurately describes cracks in beam lattices? (iii) How do the predicted fracture patterns compare with experimental observations?

For consistency, this study focuses on a single lattice type: the square beam lattice with equal beam heights in both directions, analyzed in 2D. This type of structure is unique in that the undamaged elastic Cosserat stiffness constants are available analytically, allowing validation of the numerical technique proposed here. However, the fracture behavior is distinct enough to necessitate a fairly complex homogenized model.

Figure 1 schematically illustrates the problem statement. Our goal is to test whether a unique toughness exists and to identify a homogenized phase-field model (characterized by stiffness, $g_c^M$, and $l_c^M$) that behaves similarly to the beam lattice in terms of fracture (characterized by the parameters $h^m$, $L^m$, $E^m$, and $\sigma_c^m$).

The superscript $\circ^m$ denotes a microscopic (or beam-level) quantity, while $\circ^M$ represents a macroscopic (or Cosserat) homogenized quantity.

## 3 Methods

The following sections are dedicated to detail the three elementary approaches used in this paper. First, the Euler-Bernoulli beam model is explained, then the basic theory of the phase-field model is presented. Finally, a computational homogenization technique is proposed to calculate the elastic and damaged Cosserat stiffness constants of periodic beam lattices.

### 3.1 Beam lattice model

To analyze slender beam structures, the Euler-Bernoulli beam theory offers the most appropriate framework. This classical theory is adapted to characterize the behavior of extended load-bearing elements with one side significantly larger than the others. It relies on three key assumptions:
<figure><img src="image_2.png" /><figcaption>Figure 1: Schematic illustration of the two models used to study crack propagation in beam lattices.</figcaption></figure>

(i) plane cross sections remain planar, (ii) the cross sections are perpendicular to the neutral axis, and (iii) the beam's deformation angles are small. Consequently, the equilibrium can be mathematically described by the following equations:

$$ \begin{aligned} \frac{d^2 M^m}{dt^2} + p_n^m &= 0, \\ \frac{d N^m}{dt} + p_t^m &= 0. \end{aligned} \qquad (1) $$

In this context, $M^m$ and $N^m$ denote the bending moment and normal force, while $p_n^m$ and $p_t^m$ represent distributed loads perpendicular and parallel to the neutral axis (t). These equations are complemented by Neumann and Dirichlet boundary conditions. The Euler-Bernoulli beam theory, in its simplification neglecting shear deformation effects, is particularly suited for characterizing the behavior of slender beams. The linear elastic response can be derived as follows:

$$ \begin{aligned} M^m &= EI\kappa^m, \\ N^m &= EA\epsilon^m, \end{aligned} \qquad (2) $$

Here, $\kappa^m$ and $\epsilon^m$ correspond to the two types of deformations that slender beams can undergo: (i) curvature/bending and (ii) elongation/compression. The parameters E (Young's modulus), I (bending moment of inertia around the z axis), and A (cross-sectional area) are essential in this context.

To establish the connection between degrees of freedom and deformations, we make use of the following assumptions:

$$ \kappa^m = \frac{\partial \phi^m}{\partial t} = \frac{\partial^2 u_n^m}{\partial t^2}, \quad \epsilon^m = \frac{\partial u_t^m}{\partial t}. \qquad (3) $$

Here, $\phi^m$ represents rotation, $u_n^m$ signifies perpendicular displacement, and $u_t^m$ corresponds to parallel displacement concerning the neutral axis.

To determine whether the beams have failed, we compare the maximum axial stress in the beam cross section ($\sigma_n^m$) to the material's tensile strength ($\sigma_c^m$) using the following equation:

$$ \sigma_n^m = \frac{N^m}{A} + \frac{|M^m| h^m}{I^2} \le \sigma_c^m. \qquad (4) $$

Here, $h^m$ is the height of the beam.

Given the linear nature of the problem, the load is increased until the point at which the first beam fails. After registering the force maximum, the beam's stiffness is deactivated, and the load is increased again. In the event of unstable propagation, the load remains constant, and all beam elements in which the stress is larger than $\sigma_c^m$ are iteratively deactivated to track the development of the crack. Once there are no remaining beams experiencing stress levels beyond the material's strength, the loading is resumed.
## 3.2 Phase-field fracture model

Bourdin et al. (2000) introduced the widely used theory for modeling fracture with variational methods, based on damage mechanics (Kachanov 1958) and the regularization of discontinuities with a continuous field (Ginzburg et al. 1950; Cahn et al. 1958). This theory replaced the discrete fracture surface concept from (A. A. Griffith 1921; A. Griffith 1924) with a continuous damage density function within the variational framework of Francfort et al. (1998). It utilized the Mumford and Shah functional (Mumford et al. 1989), part of the broader Ambrosio and Tortorelli elliptic regularization framework (Ambrosio et al. 1990).

In this section we give a general description of the phase-field model used in this study. Then later in Section 5 we will detail how it was used to model a square beam lattice.

### 3.2.1 Basic equation

The phase-field fracture model represents cracks with a damage variable (d) ranging from 0 (undamaged) to 1 (fully fractured). This method allows for the simulation of crack initiation and propagation without explicit crack tracking. The evolution of damage is governed by partial differential equations coupled with linear elasticity. Fundamentally, the energy of the solid body, as shown in eq. (5), is minimized:

$$ \Pi(\mathbf{u}^M, d) = \Psi_{el}(\mathbf{u}^M, \phi^M, d) + \Psi_d(d, \nabla d). \quad (5) $$

In this equation $\Psi_d(d, \nabla d)$ represents the energy consumed by the crack:

$$ \Psi_d(d, \nabla d) = \frac{g_c^M}{c_\omega l_c^M} \int_{\Omega} \left[ \omega(d) + \left(l_c^M\right)^2 |\nabla d|^2 \right] d\Omega, \quad (6) $$

with $g_c^M$ being the macroscopic fracture toughness, $l_c^M$, the internal length scale, $c_\omega$ a normalization constant. In this study we used an AT1 formulation (Pham et al. 2011), thus $\omega(d) = d$ and $c_\omega = 8/3$.

This damage is induced by the crack driving force generated by the elastic strain energy:

$$ \Psi_{el}(\mathbf{u}^M, \phi^M, d) = \int_{\Omega} \psi_{el}^M(\mathbf{u}^M, \phi^M, d) d\Omega = \int_{\Omega} \left[ g(d) \psi_0^M(\mathbf{u}^M, \phi^M) \right] d\Omega, \quad (7) $$

where $\psi_{el}^M$ is the damaged and $\psi_0^M$ is the undamaged strain energy densities and $g(d)$ is the degradation function of form which will be given later.

The mechanical and damage problems are solved in a staggered manner using a single iteration Molnár, Doitrand, Jaccon, et al. 2022; therefore, the time step was controlled automatically Molnár, Doitrand, Estevez, et al. 2020 to capture possible unstable propagation.

### 3.2.2 Multi-phase-field model

Anisotropy in the failure behavior can be induced in a hard and soft manner. The hard includes a penalty function in the damage energy, the soft introduces the degradation of the stiffness in an anisotropic manner. The advantage of the first is that the crack can be enforced to propagate in a direction relatively easily. The inconvenience is that it is not possible to introduce different fracture toughness values in different directions, which may reflect the physical reality. Moreover, if we neglect the gradient of the damage and calculate the homogeneous solution, the maximum stress remains direction-independent, and the penalty function is removed.

Therefore, in this paper, we introduce anisotropy in a soft manner and use multiple damage variables ($d_i$) for different lattice orientations. As a result, in each direction, we can choose the fracture toughness and length scale associated with the beam structure:

$$ \Psi_d(\mathbf{d}, \nabla \mathbf{d}) = \sum_{i=1}^{n} \frac{g_{c,i}^M}{c_\omega l_{c,i}^M} \int_{\Omega} \left[ d_i + \left(l_{c,i}^M\right)^2 |\nabla d_i|^2 \right] d\Omega. \quad (8) $$
In our phase-field calculation, compressive and tensile strain energies are distinguished only along the two directions parallel to the beam structure, which, in our case, align with the global axes. Specifically, if the normal deformations are less than zero (indicating compressive strains), their energies do not contribute to the crack-driving force, and, accordingly, the stiffness in these cases remains undegraded.

Damage irreversibility is enforced using Lagrange multipliers (Lu et al. 2020). More about the implementation can be found in reference (Molnár, Doitrand, Jaccon, et al. 2022).

### 3.2.3 Cosserat elasticity

The classical description of continuum mechanics is ill-suited to characterize the response of materials with an inhomogeneous microstructure, a characteristic microscopic length scale. However, the mechanical behavior of architected materials (lattice structures) is often determined by their specific micro-scale configurations. Therefore, the Cosserat theory (or micropolar elasticity) incorporates rotational degrees of freedom ($\phi$) into the mechanical description.

The Cauchy model is completed with an additional set of equations describing momentum equilibrium:

$$ 
\begin{array}{ll}
\nabla \cdot \sigma = 0 & \text{in } \Omega, \\
\nabla \cdot \mu - \hat{\epsilon} \cdot \sigma = 0 & \text{in } \Omega, \\
\sigma \cdot n = \bar{t} & \text{on } \Gamma_N, \\
\mu \cdot n = \bar{M} & \text{on } \Gamma_N, \\
\mathbf{u} = \bar{\mathbf{u}} & \text{on } \Gamma_D, \\
\phi = \bar{\phi} & \text{on } \Gamma_D.
\end{array}
\qquad (9) $$

In this equation, $\sigma$ is now a non-symmetric ($\sigma_{xy} \neq \sigma_{yx}$) force-stress tensor, $\mu$ is the moment or couple-stress tensor, and $\hat{\epsilon}$ is the three dimensional Levi-Civita symbol. The bar symbol represents external forces ($\bar{t}$), couples ($\bar{M}$), prescribed displacements ($\bar{\mathbf{u}}$), and rotations ($\bar{\phi}$).

The literature (Forest, Pradel, et al. 2001) recounts various ways to define linear elastic behavior using the Cosserat continuum. In this paper, we chose to correlate the complete stress tensor to the deformation components using the following model:

$$ \begin{bmatrix} \boldsymbol{\sigma} \\ \boldsymbol{\mu} \end{bmatrix} = \begin{bmatrix} \sigma_x \\ \sigma_y \\ \sigma_{xy} \\ \sigma_{yx} \\ \mu_x \\ \mu_y \end{bmatrix} = \mathbb{C}(\mathbf{d}) \begin{bmatrix} \boldsymbol{\varepsilon} \\ \boldsymbol{\kappa} \end{bmatrix} = \mathbb{C}(\mathbf{d}) \begin{bmatrix} \varepsilon_x \\ \varepsilon_y \\ \varepsilon_{xy} \\ \varepsilon_{yx} \\ \kappa_x \\ \kappa_y \end{bmatrix}. \qquad (10) $$

As a result, the first four elements of the stress vector correspond to the Cauchy stress components, however paying attention that in the Cosserat case $\sigma_{xy}$ and $\sigma_{yx}$ are not necessarily equal. Finally, the last components are the couple stresses.

Similarly the deformation components can be written as a function of the macroscopic degrees of freedoms, the displacements ($\mathbf{u}$) in the $x$ and $y$ directions and the rotation around the $z$ axis ($\phi$):

$$ \varepsilon_x = \frac{\partial u_x}{\partial x}, \varepsilon_y = \frac{\partial u_y}{\partial y}, \varepsilon_{xy} = \frac{\partial u_x}{\partial y} + \phi, \varepsilon_{yx} = \frac{\partial u_y}{\partial x} - \phi, \kappa_x = \frac{\partial \phi}{\partial x}, \kappa_y = \frac{\partial \phi}{\partial y}. \qquad (11) $$

The 6 × 6 (in 2D) Cosserat stiffness tensor $\mathbb{C}(\mathbf{d})$ as a function of the damage variables is going to be discussed later.
### 3.2.4 Degradation function

To better control soft anisotropy in stiffness degradation (Scherer et al. 2022), the traditional energy degradation function $g(d) = (1 - d)^2$ is replaced by a more complex one, which reverts to the traditional formulation when necessary (Lorentz et al. 2011):

$$g(d) = \left( \frac{1-d}{1+d \cdot \gamma} \right)^2, \qquad (12)$$

where $\gamma$ is a parameter, which describes stress localization in the lattice structure from a macroscopic perspective.

The way how the degradation function acts on the elastic strain energy and the exact values of $\gamma$ are dependent on the micro-structure. The identification procedure is detailed in Section 5 for the rectangular lattice.

## 3.3 Computational homogenization

One of the main disadvantage of the Cosserat theory, that its rigidity is available only for a few types of lattices. Fortunately, for a rectangular beam lattice it can be calculated analytically. To identify the constants in the degradation function described in Section 3.2, we analyze the damaged structure, for which the stiffness constants can no longer be determined analytically. Consequently, we employ a computational homogenization algorithm to calculate the residual stiffness.

We employed a two-step homogenization process to determine the macroscopic Cosserat constants for a specific beam lattice and its representative volume element (RVE).

First, we conducted discrete beam calculations on an RVE with periodic boundary conditions, using the resulting beam forces and moments to calculate the local Cauchy and coupled stress tensors at the junction nodes (Liebenstein et al. 2018a). Next, employing a second-order computational homogenization framework (Geers et al. 2001), we extracted the average stresses and strains. Finally, by correlating the deformation measures with the stress responses, we derived the elastic Cosserat constants.

### 3.3.1 Local stresses

<figure><img src="image_3.png" /><figcaption>Figure 2: Computational Cosserat homogenisation of elastic properties of Euler-Bernoulli beam lattices.</figcaption></figure>

Local stress tensors are derived from the discrete beam forces and moments using the principle of virtual work to establish energy-consistent stress fields within each control volume (Liebenstein et al. 2018a, pp.5-7). According to this principle, the internal work done by forces and moments within the discrete beam network of a control volume $V_i$ (as shown in Figure 2) must be equal to the internal work done by an equivalent continuum stress field over the same volume. To achieve this, we calculate the virtual work associated with nodal displacements and rotations in each control volume, ensuring energy equivalence between the discrete beam model and the continuum representation.
The Cauchy stresses at the i-th junction node can then be obtained using the following equation:

$$ \sigma_i^{\text{meso}} = \frac{1}{V_i} \sum_{j=1}^{n_c} \mathbf{F}_j^m \otimes \mathbf{l}_j, \qquad (13) $$

where $F_j^m$ represents the forces acting on a beam element (assumed constant along the beam) in a Cartesian coordinate system, and $l_j$ is the middle point vector, defined as the difference between the midpoint of two junctions and the junction node itself. This sum runs up to $n_c$, covering all beam segments connected to the junction node. The symbol $\otimes$ denotes the outer product, which produces a matrix by multiplying two vectors, preserving directional relationships between them.

Since the forces are constant along each beam and their directional sum is zero at each junction, it can easily be shown that the divergence of Equation (13) is zero in both directions. Thus, the equilibrium in Equation (9) is respected locally as well.

The Cauchy stresses account only for the contribution of the normal forces. The contribution of the moments, however, is represented by the couple stress tensor, denoted as $\boldsymbol{\mu}_i^{\text{meso}}$, which can be calculated at junction i using the following expression:

$$ \boldsymbol{\mu}_i^{\text{meso}} = \frac{1}{V_i} \sum_{j=1}^{n_c} [\mathbf{M}_j^m - (\mathbf{F}_i^m \times \mathbf{l}_j)] \otimes \mathbf{l}_j \stackrel{2D}{=} \frac{1}{V_i} \sum_{j=1}^{n_c} (\mathbf{M}_j^m - \|\mathbf{F}_i^m \times \mathbf{l}_j\|) \otimes \mathbf{l}_j, \qquad (14) $$

where $\mathbf{M}_j^m$ is the moment vector of the beam. In 2D, this reduces to a scalar moment $\mathbf{M}_j^m$ around the z-axis (denoted $M_z$). The magnitude is used in the second term because the cross product $\mathbf{F}_i^m \times \mathbf{l}_j$ is perpendicular to the xy-plane.

In these calculations, we assume that, (i) within a given control volume (and thus at a given node), the stress fields are homogeneous; (ii) each beam contributes equally to both of its connecting nodes; and (iii) nodes have no spatial extension—the volume is used solely for the final step in the weighting process.

We note that by calculating local, mesoscopic strain values (deformations and curvatures) at each junction node, it would be possible to obtain local elastic constants associated with these nodes. However, this is beyond the scope of this paper, and interested readers are referred to the relevant section of the article (Liebenstein et al. 2018a, p.8).

### 3.3.2 Macroscopic boundary conditions and deformations

The limitations of first-order deformation approximations become evident when modeling Cosserat media, as these approximations assume uniform deformation fields across the RVE. Cosserat media, with their intrinsic rotational and bending effects, require a second-order homogenization approach. Geers et al. (2001) and Kouznetsova et al. (2002) proposed a second-order computational homogenization scheme that incorporates the second term in the Taylor series expansion to approximate the deformed state of an infinitesimal material line element:

$$ \mathbf{x}_j = (\nabla \mathbf{u}^M + \mathbf{I}) \cdot \mathbf{X}_j + \frac{1}{2} {}^3\Phi^M \cdot (\mathbf{X}_j \otimes \mathbf{X}_j), \qquad (15) $$

where $\nabla \mathbf{u}^M$ is the applied displacement gradient, and ${}^3\Phi^M$ is a third order tensor representing second-order deformations (Geers et al. 2001; Kouznetsova et al. 2002), which arise from the finite size of the RVE. It is calculated as ${}^3\Phi^M = \nabla (\nabla \mathbf{u}^M)$, and contains only the applied curvatures in our case (Forest and Sab 1998). The symbol $\otimes$ denotes the outer product.

To simulate periodic, uniform deformations, the RVE is deformed using periodic boundary conditions, enforced by constraining the degrees of freedom between opposite boundary nodes through Lagrange multipliers. As illustrated in Figure 2, for boundary nodes j, the displacement difference can be calculated as follows:

$$ (\mathbf{u}_j^2 - \mathbf{u}_j^1) = (\nabla \mathbf{u}^M + \mathbf{I}) \cdot d\mathbf{X}_j + \frac{1}{2} {}^3\Phi^M \cdot (d\mathbf{X}_j \otimes d\mathbf{X}_j) - d\mathbf{X}_j \qquad (16) $$
To apply bending, rotations are prescribed as follows:

$$ (\phi_j^2 - \phi_j^1) = d\mathbf{X}_j^T \cdot \mathbf{\kappa}^M, \qquad (17) $$

with $d\mathbf{X}_j = \mathbf{X}_j^2 - \mathbf{X}_j^1$ and $\mathbf{\kappa}^M$ representing the prescribed curvature tensor.

To apply a homogeneous rotation field ($\phi^M$), the junction nodes ($n_j$) were rotated clockwise. Thus, following the work of Liebenstein et al. (2018), the macroscopic Cosserat strain tensor $\boldsymbol{\epsilon}^M$ (see details in Section 3.2.3), can be obtained the traditional manner using the following equations:

$$ \boldsymbol{\epsilon}^M = \nabla \mathbf{u}^M + \hat{\boldsymbol{\epsilon}} \frac{1}{V} \sum_{i=1}^{n_{\text{jint}}} \phi_i V_i^m = \nabla \mathbf{u}^M + \hat{\boldsymbol{\epsilon}} \boldsymbol{\phi}^M, \qquad (18) $$

where $V$ is the total volume of the RVE, $n_{\text{jint}}$ is the number of junctions in the RVE, $\phi_i$ are the rotations of these junction nodes, $V_i$ are the volumes defined by the Voronoi cell around each junction (shown in Figure 2a with transparent blue) and $\hat{\boldsymbol{\epsilon}}$ is the Levi-Civita tensor.

### 3.3.3 Macroscopic stresses

In the second step of our homogenization, we define energy-equivalent stresses at the macroscopic scale (the scale of the RVE) using the Hill-Mandel energy condition. This condition states that the mean energy at the microscopic scale must equal the energy density of the RVE at the macroscopic scale. For the current problem, the Hill-Mandel condition is expressed as:

$$ \frac{1}{V} \int_V [\boldsymbol{\sigma}^{\text{meso}} : \boldsymbol{\epsilon}^{\text{meso}} + \boldsymbol{\mu}^{\text{meso}} : \boldsymbol{\kappa}^{\text{meso}}] dV = \boldsymbol{\sigma}^M : \boldsymbol{\epsilon}^M + \boldsymbol{\mu}^M : \boldsymbol{\kappa}^M + {}^3\boldsymbol{Q}^{M:3} \boldsymbol{\Phi}^M. \qquad (19) $$

Applying static equilibrium and conservation of moments in the microstructure (Equation (9)), and assuming no body forces, we can use the divergence theorem on the left side, yielding the following expression for incremental work:

$$ \delta W^{\text{meso}} = \frac{1}{V} \int_{\Gamma} \left[ \mathbf{n}^{\text{T}} \cdot (\boldsymbol{\sigma}^{\text{meso}})^{\text{T}} \cdot \delta \mathbf{x} + \mathbf{n}^{\text{T}} \cdot (\boldsymbol{\mu}^{\text{meso}})^{\text{T}} \cdot \delta \boldsymbol{\phi} \right] d\Gamma + = \frac{1}{V} \int_{\Gamma} \left[ \mathbf{t}^{\text{T}} \cdot \delta \mathbf{x} + \mathbf{M}^{\text{T}} \cdot \delta \boldsymbol{\phi} \right] d\Gamma, \qquad (20) $$

where $\mathbf{n}$ denotes the normal vector on the boundaries, and $\mathbf{t}$ and $\mathbf{M}$ represent surface traction and surface couple, respectively. Substituting from Equation (15) and using the curvature definition from Equation (11), we obtain:

$$ \delta W^{\text{meso}} = \frac{1}{V} \int_{\Gamma} [\mathbf{t} \otimes \mathbf{X}] d\Gamma : \delta \nabla \mathbf{u}^{\text{M}} + \frac{1}{2V} \int_{\Gamma} [\mathbf{t} \otimes \mathbf{X} \otimes \mathbf{X}] d\Gamma : \delta^3 \boldsymbol{\Phi}^{\text{M}} + \frac{1}{V} \int_{\Gamma} [\mathbf{M} \otimes \mathbf{X}] d\Gamma : \delta \nabla \boldsymbol{\phi}^{\text{M}}. \qquad (21) $$

By converting these contour integrals back into volumetric form, we derive the energy-equivalent stresses as follows:

$$ \boldsymbol{\sigma}^M = \frac{1}{V} \int_V \boldsymbol{\sigma}^{\text{meso}} dV, \qquad (22) $$

$$ {}^3\boldsymbol{Q}^M = \frac{1}{2V} \int_V \left[ (\boldsymbol{\sigma}^{\text{meso}})^T \otimes \mathbf{X} + \mathbf{X} \otimes \boldsymbol{\sigma}^{\text{meso}} \right] dV, \qquad (23) $$
$$ \mu^M = \frac{1}{V} \int_V \mu^{\text{meso}} dV. \qquad (24) $$

The detailed simplifications are given in Kouznetsova et al. (2002) and in Appendix A. Substituting Equation (13) and Equation (14) into Equations (22)-(24) allows us to express macroscopic stress measures from the discrete microscopic beam elements as follows:

$$ \sigma^M = \frac{1}{V} \sum_{i=1}^{n_{\text{beam}}} \mathbf{F}_i^m \otimes \mathbf{l}_i^e, \qquad (25) $$

$$ {}^3Q^M = \frac{1}{2V} \sum_{i=1}^{n_{\text{beam}}} \left[ (\mathbf{F}_i^m \otimes \mathbf{l}_i^e)^T \otimes \mathbf{X}_{C,i} + \mathbf{X}_{C,i} \otimes (\mathbf{F}_i^m \otimes \mathbf{l}_i^e) \right], \qquad (26) $$

$$ \mu^M = \frac{1}{V} \sum_{i=1}^{n_{\text{beam}}} \left[ \mathbf{M}_i^m - \left( \mathbf{F}_i^m \times \frac{1}{2} \mathbf{l}_i^e \right) \right] \otimes \mathbf{l}_i^e \stackrel{2D}{=} \frac{1}{V} \sum_{i=1}^{n_{\text{beam}}} \left( M_i^m - \left\| \mathbf{F}_i^m \times \frac{1}{2} \mathbf{l}_i^e \right\| \right) \otimes \mathbf{l}_i^e, \qquad (27) $$

where $\mathbf{l}_i^e$ is the beam vector, and $\mathbf{X}_{C,i}$ denotes the coordinates of the beam element midpoint. With no body forces on the beams, a single beam element between two junction nodes is sufficient, making $\mathbf{l}^e = 2\mathbf{l}$.

We assume that the second-order term does not contribute to the constitutive response in homogenization. Were we to replace the beam lattice with a homogeneous Cosserat model under simple bending, second-order deformations and stresses would arise independently, and the energy would scale with model size (e.g., multiple unit cells). Thus, second-order deformation effects persist irrespective of the model choice.

### 3.3.4 Constitutive response

Following the stiffness calculation based on a Cauchy continuum for atomic-scale samples (Molnár, Ganster, et al. 2016), we present an analogous description for the Cosserat theory applied to beam lattices.

To calculate the Cosserat rigidity of the lattice structure the RVE is deformed in six different ways similar to the deformations of the Cosserat model: 2 axial extensions ($\epsilon_x^M, \epsilon_y^M$), 2 asymmetric shears ($\epsilon_{xy}^M, \epsilon_{yx}^M$), and 2 curvatures ($\kappa_x^M, \kappa_y^M$).

Six different quasi-static deformation cases result in 36 equations. The stiffness matrix is symmetric. Therefore, in 2D, an anisotropic material can be described by 21 unknowns with the Cosserat theory. The six individual equations were rewritten as an overdetermined equation system:

$$ \|\mathbf{M}\mathbf{c} - \mathbf{s}\| = \min \qquad (28) $$

which, relates the Cosserat moduli to the stresses (Equation (25) and Equation (27)). The coefficient matrix ($\mathbf{M}$, size: 36×21) contains six blocks with the applied strain values for each deformation case. The stiffness components are the unknowns ($\mathbf{c}$, size: 21×1):

$\mathbf{c} = [ C_{11} \quad .. \quad C_{16} \quad C_{22} \quad .. \quad C_{26} \quad C_{33} \quad .. \quad C_{36} \quad C_{44} \quad C_{45} \quad C_{46} \quad C_{55} \quad C_{56} \quad C_{66} ]^T$, and the stress values are the constant terms in six blocks for each deformation case ($\mathbf{s}$, size: 36×1):

$\mathbf{s} = [ .. \quad \sigma_x^{M,(j)} \quad \sigma_y^{M,(j)} \quad \sigma_{xy}^{M,(j)} \quad \sigma_{yx}^{M,(j)} \quad \mu_x^{M,(j)} \quad \mu_y^{M,(j)} \quad .. ]^T$. To solve the overdetermined system QR decomposition was used.
<figure><img src="image_4.png" /><figcaption>Figure 3: (a) Plate with a central crack of length $2a_0$, composed of a rectangular beam lattice with beam dimensions $L^m$ and $h^m$. (b) One-quarter model of the symmetric geometry with size $L^M$, with broken beams highlighted in red. (c) Energy release rate $G^M$ as a function of normalized crack advancement.</figcaption></figure>

## 4 Lattice toughness

Fracture in beam lattices is defined by the failure of individual beams. We consider beams failed, when the maximum axial stress defined in Equation (4) equal or overcomes their elementary strength $\sigma_c^m$. In this case, their stiffness, thus the stress which they support goes to zero. Basically, we deactivate the element, only leaving a very small residual stiffness to avoid the singularity of the global stiffness matrix of the model.

The length of the crack (a) is measured by the incremental position of the furthest beam broken. The advantage of this consideration is that it remains consistent with a macroscopic (continuum) approach, where we follow a crack tip. The disadvantage, is that we have trouble following branched cracks and potential crack widening.

Nevertheless, if the crack length is defined, the differential energy release rate ($G^M$) is expressed by the following energy balance equation:

$$G^M (a + \Delta a/2) = - \frac{\Pi_{\text{int}} (a + \Delta a, P) - \Pi_{\text{int}} (a, P) - \Delta \Pi_{\text{ext}} (P)}{\Delta a}, \quad (29)$$

where $\Pi_{\text{int}}$ is the elastic strain energy, $\Delta\Pi_{\text{ext}}$ is the external work, and P represents the applied displacements or external forces on the boundaries. Here, a is the initial crack length, and $\Delta a$ is the crack increment.

Let us consider a simple problem depicted in Figure 3(a). A plate with a $2a_0$ length crack in its middle is subjected to tensile stress on its upper and lower boundaries. While the perpendicular sides are left free to displace. The plate is constructed from a rectangular beam lattice with an elementary beam length of $L^m$ and a beam height of $h^m$. The beams have a Young's modulus of E and a failure strength of $\sigma^m$. As the problem is symmetric, we only model one forth of the geometry depicted in Figure 3(b). The size of the model is taken as $L^M$ with a crack length in the bottom left corner of $a_0$. On the left and bottom sides (except along the crack) symmetric boundary conditions are defined in both the displacement and rotation degrees of freedoms. The load is applied through displacement Dirichlet boundary conditions on the top. The broken beams are highlighted in red.

Finally, Figure 3(c) illustrates $G^M$ as a function of normalized crack advancement for a model with $L^m = 5$ mm, $E = 3$ GPa, $h^m = 1$ mm, and $\sigma_c^m = 100$ MPa. The geometry used has $a_0 = 25L^m$ and $L^M = 150L^m$. The plot shows that, after crack initiation, the energy release rate $G^M$ varies by less than 10%, indicating that it can be considered approximately constant.

Figure 3(c) shows that the energy release rate $G^M$ remains relatively constant as a function of crack advancement, even when the crack shape deviates from its original form. This observation suggests that $G^M$ may be a definable property for a given lattice structure. In the following sections, we will examine how various model and material parameters influence $G^M$ and identify the trends associated with these effects. We will begin by analyzing the impact of model size, followed by an investigation of microscopic dimensions, and finally, we will explore the influence of the elementary material properties.
## 4.1 Macroscopic geometry

A. A. Griffith (1921) originally postulated that the critical energy release rate, $g_c^M$, is a material parameter. Thus, for a given material strength, $g_c^M$ should be unique and independent of the size of the sample. The structural and material properties used in this test are $L^m = 5$ mm, $h^m = 1$ mm, $E = 3$ GPa and $\sigma_c^m = 100$ MPa.

Figure 5 shows $G^M$ as a function of normalized crack advancement for different model sizes with a fixed initial crack length of $a_0 = 25L^m$.

<figure><img src="image_5.png" /><figcaption>Figure 4: Effect of sample with $a_0 = 25L^m$ on the critical energy release rate.</figcaption></figure>

The propagation was followed until the same absolute crack length.

It is clearly visible that after an initial minor fluctuation, the values of $G^M$ are relatively independent of the macroscopic geometry. However, it is important to note that if $L^M$ or $a_0$ approaches $L^m$, scale separation cannot be guaranteed. As a result, the definition of $G^M$ becomes more difficult.

## 4.2 Microstructure

The microstructure significantly influences the homogenized elastic stiffness of a beam lattice (Pradel et al. 1998; Sab et al. 2009). In this section, with a fixed model ratio at $L^M = 150L^m$, an initial crack length of $a_0 = 25L^m$, and material properties $E = 3$ GPa and $\sigma_c^m = 100$ MPa, we will examine the effects of beam height ($h^m$) and slenderness ($L^m/h^m$) on the energy release rate $G^M$.

Figure 5(a) shows the variation of $G^M$ along the propagating crack as a function of the elementary beam height for constant slenderness. It is evident that the energy release rate depends linearly on $h^m$ (with $R^2 = 1.00$). During propagation, $G^M$ remains independent of $a$, as shown previously.

While Figure 5(b) illustrates the effect of slenderness with constant $h^m = 1$ mm. Interestingly, the influence of slenderness is much smaller compared to $h^m$. As slenderness increases, the influence diminishes and $G^M$ converges to a constant value. The mean values of $G^M$ were fitted with a hyperbolic function (with $R^2 = 0.996$).

Thus, we conclude that the beam height $h^m$ is the primary factor affecting $G^M$, while slenderness has only a minor, negligible impact on the energy release rate. This contrasts with its effect on the elastic homogeneous behavior, where $L^m/h^m$ had a significant influence (Pradel et al. 1998; Sab et al. 2009).

Although not shown here, we have found that also the type of structure significantly affects $G^M$. However, due to length constraints, this analysis will be presented in a separate paper.

## 4.3 Material properties

Finally, we turn our attention to the microscopic material properties, which are also known to influence the homogenized properties of the lattice. We use a model size of $L^M = 150L^m$ and an initial crack length of $a_0 = 25L^m$, with a microscopic beam length of $L^m = 5$ mm and beam height of $h^m = 1$ mm.

Firstly, Figure 6(a) demonstrates the effect of failure strength. As expected, $\sigma_c^m$ is the primary variable influencing the fracture toughness of the material. The energy release rate depends quadratically on the tensile strength (with $R^2 = 1$).
<figure><img src="image_6.png" /><figcaption>**Figure 5:** Effect of micro-structure: (a) Beam height with a slenderness ratio of $L^m/h^m = 5$; (b) Slenderness with $h^m = 1$ mm.</figcaption></figure>

Figure 6(b) shows the effect of Young's modulus. A hyperbolic dependence was found (with $R^2 = 1$).

<figure><img src="image_7.png" /><figcaption>**Figure 6:** Effect of material properties: (a) Tensile strength; (b) Young's modulus.</figcaption></figure>

## 4.4 Beam summary

Fracture behavior in beam lattices is defined by the failure of individual beams when the axial stress exceeds their critical strength, $\sigma_c^m$. Once a beam fails, it no longer contributes to the structural stiffness. For a given set of structural and material parameters the energy release rate, $G^M$, remains nearly constant as the crack advances. This invariance suggests that $G^M$ is a characteristic property of the lattice structure itself.

The impact of macroscopic geometry on $G^M$ is minimal, with significant changes only observed when the geometry approaches the scale of the individual beams. This stability of $G^M$ with respect to macroscopic dimensions highlights its role as an intrinsic characteristic of the lattice.

Microstructural and material properties have a more pronounced effect on $G^M$. The beam height $h^m$ significantly influences $G^M$, showing a strong linear relationship. In contrast, the slenderness has a smaller, vanishing effect. Material properties also affect $G^M$: it scales quadratically with the failure strength $\sigma_c^m$ and exhibits a hyperbolic dependence on Young's modulus $E$. These results are consistent with theoretical models Irwin (1958) and phase-field studies (Molnár, Doitrand, Estevez, et al. 2020; Molnár, Doitrand, and Lazarus 2024) on solid, homogeneous materials, suggesting a correlation between discrete microscopic parameters and homogenized phase-field characteristics.
<figure><img src="image_8.png" /><figcaption>Figure 7: Normalized Cosserat stiffness coefficients as a function of the reduction in height of the horizontal beams.</figcaption></figure>

# 5 Phase-field homogenization

The previous section demonstrated that the energy released in rectangular beam lattice is quite constant respect to the length of the crack opened. This result encouraged us to define a continuum model capable of reproducing the crack patterns observed in the discrete beam lattice.

With the anisotropic multi-phase-field model (presented in Section 3.2), the key is to correctly identify how many damage variables do we need, and which variable acts on which stiffness. This way it is going to be clear, how to calculate the crack driving force to induce damage in a given direction.

The square lattice is one of the rare examples, where the undamaged Cosserat stiffness is available analytically. The stiffness matrix results in a diagonal matrix with constants (Pradel et al. 1998; Sab et al. 2009):

$$ 
\begin{aligned}
C_{11} &= C_{22} = \frac{h^m}{L^m} E, \\
C_{33} &= C_{44} = \left(\frac{h^m}{L^m}\right)^3 E, \\
C_{55} &= C_{66} = \frac{(h^m)^3}{12L^m} E.
\end{aligned} 
\qquad (30) $$

We note that these constants are available in this analytic form only for the undamaged state, with equal lattice spacing and equal beam heights in the principal two directions. For cases, when the beams start to break or are already broken, we will use the computational homogenization technique described in Section 3.3. Nevertheless, these constants provide a valuable tool for verifying the numerical technique.

From the structural configuration of a rectangular grid it is easy to see, that there are two principal directions. Thus, we will use two independent damage variables ($n = 2$ in eq. (8)). $d_1$ will correspond to the horizontal grid, while $d_2$, the vertical grid direction.

## 5.1 Damage coupling

Until this point we have not yet discussed how damage variables act on the stiffness matrix. It is critical to do so, as this will decide which deformation will induce which damage.

To analyze the effect of damage on the stiffness properties, we considered an RVE of a rectangular lattice. We sequentially removed the horizontal and vertical beams and recalculated the Cosserat stiffness using the computational homogenization method described in Section 3.3.

Figure 7 illustrates the variation in stiffness as the height of the horizontal beams is reduced. The parameter $h_x^m$ was progressively decreased to simulate damage ($d_1$). As expected, completely removing the horizontal beams causes the stiffness components $C_{11}$ and $C_{55}$ to approach zero. However, it is notable that $C_{44}$ also decreases significantly, indicating that its contribution is lost as well. The remaining stiffness coefficients remain unaffected.

When the vertical beams are removed, a similar pattern emerges in the orthogonal direction: the components $C_{22}$, $C_{33}$, and $C_{66}$ are reduced to zero. These observations suggest that the degradation functions influence the stiffness components in the following manner:
$$ \mathbb{C} = \begin{bmatrix} g(d_1, \gamma_{11}) C_{11} & & & & & \\ & g(d_2, \gamma_{22}) C_{22} & & & & \\ & & g(d_2, \gamma_{23}) C_{33} & & & \\ & & & g(d_1, \gamma_{14}) C_{44} & & \\ & \text{sym.} & & g(d_1, \gamma_{55}) C_{55} & & \\ & & & & g(d_2, \gamma_{66}) C_{66} & \end{bmatrix}, \quad (31) $$

with $\gamma_{ij}$ values representing stress localization. The identification process is presented in the next section.

## 5.2 Homogeneous solution

The traditional phase-field approach needs two parameters to be identified. One is the fracture toughness, the other one is the internal length scale. As a consequence two types of tests are recommended. One represents the "flawless" behavior, when there are no cracks present. the second measures the response when there is a crack.

In the first test, the homogeneous phase-field solution is solved. This step involves neglecting the gradient of the damage in eq. (8), then taking the derivative respect to $d_i$ and solving the resulting equation. In the following we will demonstrate it for $d_1$, but the equivalent method was used for $d_2$. The aim is to define an acceptable $g_{c,1}^M/l_{c,1}^M$ ratio and the $\gamma_{1i}$ localization constants.

$$ \frac{\partial \Pi}{\partial d_1} = 0 \rightarrow \frac{\partial \psi_{el}^M}{\partial d_1} + \frac{3g_{c,1}^M}{8l_{c,1}^M} = 0. \quad (32) $$

In the AT1 model, the maximum stress is reached when $d_1 = 0$. Thus with: $\frac{\partial g(d_1=0,\gamma_{1j})}{\partial d_1} = -2(1+\gamma_{1j})$, we get:

$$ \sum_j 2(1 + \gamma_{1j}) \psi_{0,j}^{cr} = \frac{3g_{c,1}^M}{8l_{c,1}^M}, \quad (33) $$

where $\psi_{0,j}^{cr}$ are critical deformation energies which were identified using the RVE of the beam model. To identify $\gamma_{1j}$ we assumed that the phase-field model is supposed to break at the same state as the RVE in each independent deformation. Thus for example for pure unidirectional extension in the x direction we can write:

$$ (1 + \gamma_{11}) C_{11} (\varepsilon_{x,cr}^M)^2 = \frac{3g_{c,1}^M}{8l_{c,1}^M}, \quad (34) $$

with $\varepsilon_{x,cr}^M$ being the critical deformation applied on the RVE, when the local maximum tensile stress reaches the microscopic tensile strength $\sigma_c$:

$$ \varepsilon_{x,cr}^M = \frac{\sigma_c}{E}. \quad (35) $$

For the rectangular grid, this value can be identified analytically for $\varepsilon_x^M$, however could be and was calculated numerically for arbitrary structures and deformations.

By applying the same procedure on each deformation case, assuming that in each step the $g_{c,1}^M/l_{c,1}^M$ are equal, and normalizing each row, that $\gamma_{11} = 0$ and $\gamma_{22} = 0$, we can identify $\gamma$ to be:

$$ \gamma = \begin{bmatrix} 0 & 0 & 0 & 8 & 2 & 0 \\ 0 & 0 & 8 & 0 & 0 & 2 \end{bmatrix}, \quad (36) $$
<figure><img src="image_9.png" /><figcaption>Figure 8: Homogeneous solution of the phase-field model compared with the beam results.</figcaption></figure>

where each row corresponds to the $i^{th}$ damage variable and each column corresponds to the given member of the stiffness tensor (now diagonal).

From these results it is clearly visible, that shear for example generates a stress peak 3 times larger, than unidirectional extension in the rectangular lattice. This value corresponds exactly to the square root of the ratio $1 + \gamma_{14}$, which associated with the energy difference in the phase-field model between extension and shear.

Of course these parameters would have to be recalculated for a different lattice type, and $\gamma$ would potentially become longer as matrix $C$ might have off-diagonal elements as well.

The homogeneous solution calculated using the beam model and the phase-field formulation is compared in Figure 8 for various load angles. In the beam model, the RVE is subjected to an extension along the direction $\omega$ and an equal magnitude of compression in the perpendicular direction. This deformation represents a pure shear mode in homogeneous materials. Our aim is to evaluate the structural response using both methods to understand how the resistance changes between pure and combined deformation modes, depending on the loading direction relative to the structural orientation.

It is evident that the two methods show good agreement at pure deformation angles, such as $\theta = 0^\circ, 45^\circ, 90^\circ$, which can be attributed to the calibrated constants of $\gamma$. For mixed deformation modes, however, a slight difference is observed. This discrepancy arises because, in the beam model, stresses are summed directly, while in the phase-field approach, the energies are summed.

## 5.3 Effect of initial crack length

The homogeneous solution presented in Section 5.2 helps to identify appropriate ratios for $g_c^M / l_c^M$. However, to fully calibrate these parameters, an additional test case is required. Typically, in phase-field simulations, one uses a scenario where size effects are minimal (no defects) and a case involving a sufficiently large crack, as in classical fracture mechanics (Griffith-type cases).

To investigate the impact of initial defect size, the maximum force was recorded at two key points: when the crack initiated and when the sample reached its maximum load-bearing capacity. Figure 9 shows the moment of the first beam fracture with red circles and the final load-bearing capacity with red crosses. The critical loading is normalized using the homogeneous solution:

$$ \sigma_{c,y}^{M} = \sigma_{c} \frac{h^{m}}{L^{m}}. \qquad (37) $$

When $a_0 < L^m$, the critical loading matches the homogeneous solution since the crack is smaller than the beam spacing. As $a_0$ increases, a power-law behavior is observed in both crack initiation and maximum loading, aligning with linear elastic fracture mechanics. We note that a similar dependence on the crack length was observed regarding the load at failure in three-point bending experiments for cordierite square lattices (Quintana-Alonso et al. 2010).

Thus, by setting $l_c^M$ equal to $L^m$, the phase-field results closely replicate those observed in the beam model. It is important to note that, due to the gradual appearance of damage in the model, pinpointing the exact loading state corresponding to the first beam fracture is challenging. For the remainder of the analysis, we therefore set $l_{c,i}^M = L^m$.
<figure><img src="image_10.png" /><figcaption>Figure 9: Critical loading normalized by the homogenized macroscopic strength as a function of the initial crack length.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 10: (a) Experimental geometry. (b) Principal loading directions and structural orientations in a square beam lattice.</figcaption></figure>

When eq. (35), (30) and $L^m$ is substituted into eq. (34), we get:

$$g_{c,1}^{M} = \frac{8h^{m}\sigma_{c}^{2}}{3E}. \qquad (38)$$

From this equation we recover the correlations found in the beam simulations. The critical fracture toughness has indeed a quadratic dependence on the tensile strength, a hyperbolic dependence on Young's modulus, and a linear dependence on the beam height. Interestingly, the correlation identified through the homogeneous phase-field solution fails to capture the vanishing effect of slenderness observed in the beam model.

## 6 Experimental validation

The phase-field model was initially calibrated to the beam model for tensile opening, but beam-architected materials are inherently anisotropic and may exhibit different behaviors under varying loading directions. To investigate the primary mechanical couplings present in a rectangular beam lattice, we conducted a series of tests inspired by the work of Ayatollahi et al. (2009). The basic concept involves extending a rectangular plate at two of its opposite corners, as shown in Figure 10(a), which allows for the application of both tensile and shear loading in a tensile testing machine with relative ease.

In the original study by Ayatollahi et al. (2009), the load was applied through pinholes in homogeneous materials. However, due to the weakened nature of the material in our lattice structure, a concentrated load could cause the sample to fracture around the point of load application. To address this, we replaced the pinhole with a solid section, depicted in dark gray in Figure 10(a), and applied the load using clamps.

In a rectangular grid, the orthotropic microstructure introduces a third notable direction, in addition to the orientation of the loading and the crack: the orientation of the microstructure
itself. To validate our numerical models, we selected four configurations designed to test the
relative loading between these directions:

- **Case 1:** Tensile opening with microstructure parallel to the crack,

- Case 2: Shear opening with microstructure parallel to the crack,

- Case 3: Shear opening with microstructure oriented at 45° to the crack,

- Case 4: Tensile opening with microstructure oriented at 45° to the crack.

These elementary cases are illustrated in Figure 10(b).

The experimental samples were cut from 5 mm thick polymethyl methacrylate (PMMA) sheets using a laser cutter. The microscopic unit spacing was chosen to be *L*<sup>*m*</sup> = 5 mm, with a beam height of *h*<sup>*m*</sup> = 0.66 mm. The macroscopic size of the samples was *L*<sup>*M*</sup> = 141 mm. Prior to the experiments, a series of numerical tests were conducted to determine the appropriate length of the clamping zone. The clamping zone length and the initial crack length were carefully set to ensure that fracture initiation occurred at the initial crack. If the clamped zone was either too small or too large, the crack tended to initiate at the edges of the clamped zone rather than at the intended location. Consequently, we selected a clamped width of *L*<sub>c</sub> = 10*L*<sup>*m*</sup> = 50 mm and an initial crack length of *a*<sub>0</sub> = 9*L*<sup>*m*</sup> = 45 mm.

The samples were loaded gradually at a rate of 0.02 mm/s to minimize the influence of viscoelastic effects. The results from these tests are displayed in the first column of Figure 11, with the cracks highlighted in red.

To replicate the asymmetric crack initiation observed in the experiments, the tensile strength and fracture toughness in the numerical models were varied by 5 %. The results from the Euler-Bernoulli beam model are shown in the second column of Figure 11, while the corresponding phase-field simulations are presented in the third column.

The results from all models exhibit coherence and agreement, with the crack consistently favoring the direction of the original microstructure. This behavior aligns with expectations, as the shear contribution (third and fourth columns of eq. (36)) is the most significant factor influencing crack propagation.

Our analysis demonstrates that both simulation techniques accurately represent physical
reality and can be used as predictive tools for modeling fracture in lattice structures.

To further explore the limitations of our model, we conducted additional tests on two geometries with varying slenderness, focusing on the tensile case with a parallel material structure to the crack (Case 1). In the first scenario, the unit length was reduced to 2.5 mm, while keeping all other dimensions constant. In the second scenario, the overall geometry was doubled in size, with the beam height remaining unchanged.

The results of these tests are presented in Figure 12, where the crack is highlighted in red. It is evident that in the case of thick beams with a slenderness ratio of 3, the crack propagated horizontally, consistent with the behavior of solid (homogeneous) materials. This result diverged from the predictions made by both numerical methods. However, in the more slender sample, the crack followed the patterns observed in the numerical simulations, reinforcing the applicability of our models. This outcome underscores the limitations of the Euler-Bernoulli beam theory, which remains accurate up to a slenderness ratio of 5 (Molnár and Blal 2023).

7 Discussion

The paper presents a complex phase-field formulation to model fracture in beam lattices, following a step-by-step modeling strategy. We advocate starting with the simplest approach that captures the phenomenon of interest. Initially, we employed a Cosserat continuum to model the elastic response of the beam lattice, but this proved insufficient as cracks in uniaxial tension propagated horizontally, similar to bulk materials. The next step was to introduce anisotropy (specifically, orthotropy for square lattices) in the damage phase-field formulation, aiming to replicate the homogeneous beam response using a soft anisotropic approach with multiple damage variables to control toughness in different directions. However, this still did not produce the correct crack paths. Ultimately, we modified the traditional degradation function by enriching it with coefficients to capture stress localization, leading to accurate results for both homogeneous solutions and crack paths under various loading conditions.
<figure><img src="image_12.png" /><figcaption>**Figure 11:** Comparison between experiments (left) and simulations using beam theory (center) and phase-field models (right). Red highlights the crack.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 12:** Experiments on other samples with different slenderness. Cracks are highlighted in red.</figcaption></figure>
# 8 Conclusion

In this study, we examined the fracture behavior of rectangular beam lattices using finite element simulations based on Euler-Bernoulli beam theory. We employed a straightforward failure criterion where beams fracture when the maximum axial stress exceeds their critical strength, leading to progressive crack propagation as beams fail. To quantify the fracture toughness, we calculated the energy release rate $G^M$ using an energy balance approach. Our results indicate that $G^M$ remains relatively stable with crack advancement, underscoring its robustness as a material property intrinsic to the lattice structure. This stability suggests that $G^M$ is a reliable indicator of the lattice's overall toughness.

Building on this observation, we developed a multi-phase-field fracture model to more accurately represent the homogenized toughness of these architected materials. This model integrates the consistent toughness characteristics observed in our simulations into a comprehensive framework for lattice design and analysis.

To validate our theoretical and computational findings, we conducted an experimental campaign to compare the results from our simulations with physical tests. This comparison provided valuable insights into the accuracy and applicability of our models, bridging the gap between theoretical predictions and practical outcomes.

Our findings reveal that $G^M$ is largely insensitive to changes in macroscopic geometry, such as the size of the model or the initial crack length, provided these dimensions are sufficiently larger than the beam dimensions. This stability emphasizes that $G^M$ is primarily influenced by the lattice's microstructure and material properties. Specifically, we found that the beam height $h^m$ and material properties such as tensile strength and Young's modulus significantly affect $G^M$, whereas the beam length $L^m$ has a less pronounced impact. The linear, quadratic, and hyperbolic dependencies of $G^M$ on beam height, tensile strength, and Young's modulus, respectively, align with theoretical expectations and phase-field models, thereby confirming the lattice's homogenized fracture characteristics.

In subsequent phases of the study, we employed phase-field homogenization techniques to model the lattice behavior, with the goal of replicating the observed crack patterns and toughness metrics. The phase-field parameters were derived from the homogeneous phase-field solution and size-effect tests conducted in tensile opening scenarios. Our phase-field model, calibrated against the discrete beam model, demonstrated good correspondence to experimental results across various loading directions. This validation underscores the efficacy of both the beam and phase-field approaches in capturing the anisotropic fracture behavior of beam lattices, providing a reliable framework for predicting lattice performance in fracture.

We observed that, in all cases, the crack preferentially propagated along the original directions of the rectangular beam lattice. This phenomenon can be explained by stress localization, which was most pronounced in shear. This localization effect caused higher stress concentrations along specific directions, leading to preferential crack growth in those orientations. Such behavior underscores the significance of lattice orientation and loading conditions in determining fracture paths within beam lattices.

# A Macroscopic-microscopic stress relationship

In this section, we derive the relationship between Equation (21) and Equations (22)-(24). The symbol $\otimes$ denotes the outer product.

Applying the divergence theorem, the first contour integral in Equation (21) is converted to a volumetric form:

$$ \frac{1}{V} \int_{\Gamma} [t \otimes X] d\Gamma = \frac{1}{V} \int_{\Gamma} [n^T \cdot (\sigma^{\text{meso}})^T \otimes X] d\Gamma = \frac{1}{V} \int_{V} \nabla \cdot [(\sigma^{\text{meso}})^T \otimes X] dV. \quad (39) $$

Using the microscopic equilibrium, defined by the local equilibrium at each junction, $\nabla \cdot \sigma^{\text{meso}} = 0$ (see Equation (9)$_1$), and the identity $\nabla X = I$, we establish the following:
$$ \nabla \cdot [(\sigma^{\text{meso}})^T \otimes \mathbf{X}] = \nabla \cdot (\sigma^{\text{meso}})^T \otimes \mathbf{X} + \sigma^{\text{meso}} \cdot \nabla \mathbf{X} = \sigma^{\text{meso}}. \quad (40) $$

As a result, the energy contribution of the first term in Equation (21) linked to the displacement gradient simplifies to the average microscopic stress:

$$ \frac{1}{V} \int_V \sigma^{\text{meso}} dV : \delta \nabla \mathbf{u}^{\text{M}}. \quad (41) $$

The higher-order term derivation follows similarly:

$$ \begin{aligned} \frac{1}{2V} \int_{\Gamma} [\mathbf{t} \otimes \mathbf{X} \otimes \mathbf{X}] d\Gamma &= \frac{1}{2V} \int_{\Gamma} \{[\mathbf{n}^T (\sigma^{\text{meso}})^T] \otimes \mathbf{X} \otimes \mathbf{X}\}^{T_{1,2}} d\Gamma \\ &= \frac{1}{2V} \int_V \{\nabla \cdot [(\sigma^{\text{meso}})^T \otimes \mathbf{X} \otimes \mathbf{X}]\}^{T_{1,2}} dV, \end{aligned} \quad (42) $$

where the subscript $^{T_{1,2}}$ represents left conjugation, thus the permutation of indexes I and J, $a_{Ijk}^{T_{1,2}} = a_{Jik}$. Furthermore, with the equality:

$$ \begin{aligned} \nabla \cdot [(\sigma^{\text{meso}})^T \otimes \mathbf{X} \otimes \mathbf{X}] &= [\nabla \cdot (\sigma^{\text{meso}})^T] \otimes \mathbf{X} \otimes \mathbf{X} + \sigma^{\text{meso}} \cdot \nabla \mathbf{X} \otimes \mathbf{X} + (\mathbf{X} \otimes \sigma^{\text{meso}} \cdot \nabla \mathbf{X})^{T_{1,2}} \\ &= \sigma^{\text{meso}} \otimes \mathbf{X} + (\mathbf{X} \otimes \sigma^{\text{meso}})^{T_{1,2}}, \end{aligned} \quad (43) $$

we reach the expression in Equation (23). The third and final part of Equation (21) includes the diagonal couple stress tensor:

$$ \begin{aligned} \frac{1}{V} \int_{\Gamma} (\mathbf{M} \otimes \mathbf{X}) \, d\Gamma : \delta \nabla \phi^M &= \frac{1}{V} \int_{\Gamma} [\mathbf{n}^T \cdot \boldsymbol{\mu}^{\text{meso}} \otimes \mathbf{X}] \, d\Gamma : \delta \nabla \phi^M \\ &= \frac{1}{V} \int_V \nabla \cdot [\boldsymbol{\mu}^{\text{meso}} \otimes \mathbf{X}] \, dV : \delta \nabla \phi^M \\ &= \frac{1}{V} \int_V [\nabla \cdot \boldsymbol{\mu}^{\text{meso}} \otimes \mathbf{X}] \, dV : \delta \nabla \phi^M \\ &\quad + \frac{1}{V} \int_V \nabla \mathbf{X} \cdot \boldsymbol{\mu}^{\text{meso}} \, dV : \delta \nabla \phi^M. \end{aligned} \quad (44) $$

To simplify the first term, we use the fact that the geometric center is located at the origin of the Cartesian coordinate system, making the integral $\frac{1}{V} \int_V \mathbf{X} \, dV = 0$. Thus, it follows that:

$$ \frac{1}{V} \int_V [\nabla \cdot \boldsymbol{\mu}^{\text{meso}} \otimes \mathbf{X}] \, dV = 0. \quad (45) $$

While, the second term simplifies to yield the macroscopic couple stress tensor, as shown in Equation (24).

We have not yet accounted for the applied uniform micro-rotation. This case is treated with an alternate boundary condition, where the uniform rotation is applied to all nodes in the RVE. When applying a uniform rotation ($\boldsymbol{\phi} = \boldsymbol{\phi}^M$), bending cannot be applied simultaneously due to the intersection with the side boundaries, resulting in $\nabla \boldsymbol{\phi}^M = 0$. Writing the work of the external moments under these conditions and applying the divergence theorem, we obtain:
$$ 
\begin{aligned}
\frac{1}{V} \int_{\Gamma_{\phi}} [\mathbf{M}^T \cdot \delta\boldsymbol{\phi}] \, d\Gamma &= \frac{1}{V} \int_{\Gamma_{\phi}} [\mathbf{n} \cdot \boldsymbol{\mu}^{\text{meso}}] \, d\Gamma : \delta\boldsymbol{\phi}^{\text{M}} \\
&= \frac{1}{V} \int_{V} \nabla \cdot \boldsymbol{\mu}^{\text{meso}} \, dV : \delta\boldsymbol{\phi}^{\text{M}} \qquad (46) \\
&= -\hat{\varepsilon} \int_{V} \boldsymbol{\sigma}^{\text{meso}} \, dV : \delta\boldsymbol{\phi}^{\text{M}}.
\end{aligned}
$$

Considering the strain definition in Equation (18), we can show that, by including both the traditional deformations and rotations, the work-conjugate stress of the potentially asymmetric strain tensor is obtained from the volume average of the mesoscopic values.

It is important to note that if the asymmetric part of $\nabla \mathbf{u}^M$ equals $\boldsymbol{\phi}^M$, the RVE undergoes a rigid body rotation, resulting in zero associated strain and, consequently, zero stress.

## References

Ambrosio, L. and V. M. Tortorelli (1990). "Approximation of functional depending on jumps by elliptic functional via $\Gamma$-convergence". *Communications on Pure and Applied Mathematics* 43.8, pp. 999-1036. DOI: 10.1002/cpa.3160430805

Askari, M., D. A. Hutchins, P. J. Thomas, L. Astolfi, R. L. Watson, M. Abdi, M. Ricci, S. Laureti, L. Nie, S. Freear, R. Wildman, C. Tuck, M. Clarke, E. Woods, and A. T. Clare (2020). "Additive manufacturing of metamaterials: A review". *Additive Manufacturing* 36, p. 101562. DOI: 10.1016/j.addma.2020.101562

Ayatollahi, M. R. and M. R. M. Aliha (2009). "Analysis of a new specimen for mixed mode fracture tests on brittle materials". *Engineering Fracture Mechanics* 76.11, pp. 1563-1573

Bauer, J., A. Schroer, R. Schwaiger, and O. Kraft (2016). "Approaching theoretical strength in glassy carbon nanolattices". *Nature materials* 15.4, pp. 438-443. DOI: 10.1038/nmat4561

Bharali, R., F. Larsson, and R. Jänicke (2021). "Computational homogenisation of phase-field fracture". *European Journal of Mechanics-A/Solids* 88, p. 104247

Bleyer, J. and R. Alessi (2018). "Phase-field modeling of anisotropic brittle fracture including several damage mechanisms". *Computer Methods in Applied Mechanics and Engineering* 336, pp. 213-236

Bourdin, B., G. Francfort, and J.-J. Marigo (2000). "Numerical experiments in revisited brittle fracture". *Journal of the Mechanics and Physics of Solids* 48.4, pp. 797-826. DOI: 10.1016/S0022-5096(99)00028-9

Cahn, J. W. and J. E. Hilliard (1958). "Free energy of a nonuniform system. I. Interfacial free energy". *The Journal of chemical physics* 28.2, pp. 258-267

Chen, B., Q. Yuan, J. Luo, and J. H. Fan (2010). "Fibre reinforced cellular microstructure of cork wood". *Plastics, Rubber and Composites* 39.2, pp. 86-90. DOI: 10.1179/174328910X126088518326

Cosserat, E. and F. Cosserat (1909). *Théorie des corps déformables*. Librairie Scientifique A. Hermann et Fils

Dimas, L. S., G. H. Bratzel, I. Eylon, and M. J. Buehler (2013). "Tough Composites Inspired by Mineralized Natural Materials: Computation, 3D printing, and Testing". *Advanced Functional Materials* 23.36, pp. 4629-4638. DOI: https://doi.org/10.1002/adfm.201300215

Forest, S., F. Pradel, and K. Sab (2001). "Asymptotic analysis of heterogeneous Cosserat media". *International Journal of Solids and Structures* 38.26-27, pp. 4585-4608. DOI: 10.1016/S0020-7683(00)00295-X

Forest, S. and K. Sab (1998). "Cosserat overall modeling of heterogeneous materials". *Mechanics research communications* 25.4, pp. 449-454

Francfort, G. and J.-J. Marigo (1998). "Revisiting brittle fracture as an energy minimization problem". *Journal of the Mechanics and Physics of Solids* 46.8, pp. 1319-1342. DOI: 10.1016/S002-2-5096(98)00034-9
Geers, M. G. D., V. Kouznetsova, and W. A. M. Brekelmans (2001). "Gradient-enhanced computational homogenization for the micro-macro scale transition". *Le Journal de Physique IV* 11, pp. 145-152

Ginzburg, V. L. and L. D. Landau (1950). "Theory of superconductivity". *Zh. Eksp. Teor. Fiz.; (USSR)* 20.12

Griffith, A. (1924). "The theory of rupture". *First Int. Cong. Appl. Mech*, pp. 55-63

Griffith, A. A. (1921). "The Phenomena of Rupture and Flow in Solids". *Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences* 221.582-593, pp. 163-198. DOI: 10.1098/rsta.1921.0006

Irwin, R. G. (1958). "Fracture". *Elasticity and Plasticity / Elastizität und Plastizität*. Berlin, Heidelberg: Springer Berlin Heidelberg, pp. 551-590. DOI: 10.1007/978-3-662-43081-1_5

Jang, D., L. R. Meza, F. Greer, and J. R. Greer (2013). "Fabrication and deformation of three-dimensional hollow ceramic nanostructures". *Nature materials* 12.10, pp. 893-898. DOI: 10.1038/nmat3738

Jiao, P., J. Mueller, J. R. Raney, X. Zheng, and A. H. Alavi (2023). "Mechanical metamaterials and beyond". *Nature communications* 14.1, p. 6004. DOI: 10.1038/s41467-023-41679-8

Kachanov, L. M. (1958). "Rupture time under creep conditions"

Kolpakov, A. and S. Rakin (2020). "Homogenized strength criterion for composite reinforced with orthogonal systems of fibers". *Mechanics of Materials* 148, p. 103489

Kouznetsova, V., M. G. D. Geers, and W. A. M. Brekelmans (2002). "Multi-scale constitutive modelling of heterogeneous materials with a gradient-enhanced computational homogenization scheme". *International journal for numerical methods in engineering* 54.8, pp. 1235-1260

Lakes, R. (1993). "Materials with structural hierarchy". *Nature* 361.6412, pp. 511-515. DOI: 10.1038/361511a0

Lebihain, M., L. Ponson, D. Kondo, and J.-B. Leblond (2021). "Effective toughness of disordered brittle solids: A homogenization framework". *Journal of the Mechanics and Physics of Solids* 153, p. 1044-1063

Libonati, F., G. X. Gu, Z. Qin, L. Vergani, and M. J. Buehler (2016). "Bone-inspired materials by design: toughness amplification observed using 3D printing and testing". *Advanced Engineering Materials* 18.8, pp. 1354-1363

Liebenstein, S. and M. Zaiser (2018a). "Determining Cosserat constants of 2D cellular solids from beam models". *Materials Theory* 2.1, pp. 1-20. DOI: 10.1186/s41313-017-0009-x

Liebenstein, S. and M. Zaiser (2018b). "Determining Cosserat constants of 2D cellular solids from beam models". *Materials Theory* 2, pp. 1-20

Lorentz, E. and V. Godard (2011). "Gradient damage models: Toward full-scale computations". *Computer Methods in Applied Mechanics and Engineering* 200.21, pp. 1927-1944. DOI: 10.1016/j.cma.2010.06.025

Lu, Y., T. Helfer, B. Bary, and O. Fandeur (2020). "An efficient and robust staggered algorithm applied to the quasi-static description of brittle fracture by a phase-field approach". *Computer Methods in Applied Mechanics and Engineering* 370, p. 1132-1148

Molnár, G., A. Doitrand, and V. Lazarus (2024). "Phase-field simulation and coupled criterion link echelon cracks to internal length in antiplane shear". *Journal of the Mechanics and Physics of Solids* 188, p. 1056-1075. DOI: 10.1016/j.jmps.2024.105675

Molnár, G. and N. Blal (2023). "Topology optimization of periodic beam lattices using Cosserat elasticity". *Computers & Structures* 281, p. 1070-1080. DOI: 10.1016/j.compstruc.2023.107037

Molnár, G., A. Doitrand, R. Estevez, and A. Gravouil (2020). "Toughness or strength? Regularization in phase-field fracture explained by the coupled criterion". *Theoretical and Applied Fracture Mechanics* 109, p. 1027-1042. DOI: 10.1016/j.tafmec.2020.102736

Molnár, G., A. Doitrand, A. Jaccon, B. Prabel, and A. Gravouil (2022). "Thermodynamically consistent linear-gradient damage model in Abaqus". *Engineering Fracture Mechanics* 266, p. 108390. DOI: 10.1016/j.engfracmech.2022.108390

Molnár, G., P. Ganster, J. Török, and A. Tanguy (2016). "Sodium effect on static mechanical behavior of MD-modeled sodium silicate glasses". *Journal of Non-Crystalline Solids* 440, pp. 12-25
Molnár, G. and A. Gravouil (2017). "2D and 3D Abaqus implementation of a robust staggered phase-field solution for modeling brittle fracture". *Finite Elements in Analysis and Design* 130, pp. 27–38. doi: 10.1016/j.finel.2017.03.002

Mousanezhad, D., S. Babaee, H. Ebrahimi, R. Ghosh, A. S. Hamouda, K. Bertoldi, and A. Vaziri (2015). "Hierarchical honeycomb auxetic metamaterials". *Scientific reports* 5.1, pp. 1–8. doi: 10.1038/srep18306

Mumford, D. and J. Shah (1989). "Optimal approximations by piecewise smooth functions and associated variational problems". *Communications on pure and applied mathematics* 42.5, pp. 577–685. doi: 10.1002/cpa.3160420503

Nepal, D., S. Kang, K. M. Adstedt, K. Kanhaiya, M. R. Bockstaller, L. C. Brinson, M. J. Buehler, P. V. Coveney, K. Dayal, J. A. El-Awady, et al. (2023). "Hierarchically structured bioinspired nanocomposites". *Nature materials* 22.1, pp. 18–35. doi: 10.1038/s41563-022-01384-1

Nguyen, T. T., J. Réthoré, and M.-C. Baietto (2017). "Phase field modelling of anisotropic crack propagation". *European Journal of Mechanics-A/Solids* 65, pp. 279–288

Pham, K., H. Amor, J.-J. Marigo, and C. Maurini (2011). "Gradient damage models and their use to approximate brittle fracture". *International Journal of Damage Mechanics* 20.4, pp. 618–652. doi: 10.1177/1056789510386852

Pichler, B. and C. Hellmich (2011). "Upscaling quasi-brittle strength of cement paste and mortar: A multi-scale engineering mechanics model". *Cement and Concrete Research* 41.5, pp. 467–476

Pradel, F. and K. Sab (1998). "Cosserat modelling of elastic periodic lattice structures". *Comptes Rendus de l'Académie des Sciences-Series IIB-Mechanics-Physics-Astronomy* 326.11, pp. 699–704. doi: 10.1016/S1251-8069(98)80002-X

Quintana-Alonso, I., S. Mai, N. Fleck, D. Oakes, and M. Twigg (2010). "The fracture toughness of a cordierite square lattice". *Acta materialia* 58.1, pp. 201–207

Ritchie, R. O. (2011). "The conflicts between strength and toughness". *Nature Materials* 10.11, pp. 817–822. doi: 10.1038/nmat3115

Ritchie, R. O., M. J. Buehler, and P. Hansma (2009). "Plasticity and toughness in bone". *Physics today* 62.6, pp. 41–47. doi: 10.1063/1.3156332

Roux, S. and F. Hild (2008). "Self-consistent scheme for toughness homogenization". *International journal of fracture* 154.1, pp. 159–166

Sab, K. and F. Pradel (2009). "Homogenisation of periodic Cosserat media". *International journal of computer applications in technology* 34.1, p. 60. doi: 10.1504/IJCAT.2009.022703

Scherer, J.-M., S. Brach, and J. Bleyer (2022). "An assessment of anisotropic phase-field models of brittle fracture". *Computer Methods in Applied Mechanics and Engineering* 395, p. 115036

Vasoya, M., V. Lazarus, and L. Ponson (2016). "Bridging micro to macroscale fracture properties in highly heterogeneous brittle solids: weak pinning versus fingering". *Journal of the Mechanics and Physics of Solids* 95, pp. 755–773

Zheng, X., W. Smith, J. Jackson, B. Moran, H. Cui, D. Chen, J. Ye, N. Fang, N. Rodriguez, T. Weisgraber, and C. M. Spadaccini (2016). "Multiscale metallic metamaterials". *Nature materials* 15.10, pp. 1100–1106. doi: 10.1038/nmat4694
# B Backmatter information

The following subsections will be added to the accepted version of the paper.

**Funding** This work was partly supported by the French Research National Agency program through the grant ANR-16-CE30-0007-01.

**Authors' contributions** G. Molnár carried out most of the study, performed numerical simulations, and drafted the manuscript. J. Réthoré conceived of the study and helped with the experimental testing. Both authors participated in its design, coordination, and critical review of the manuscript. Both authors read and approved the final manuscript.

**Competing interests** The authors declare that they have no competing interests.