Strength increase during ceramic biomaterial-induced bone regeneration: a micromechanical study

Bone tissue engineering materials must blend in the targeted physiological environment, in terms of both the materials’ biocompatibility and mechanical properties. As for the latter, a well-adjusted stiffness ensures that the biomaterial’s deformation behavior fits well to the deformation behavior of the surrounding biological tissue, whereas an appropriate strength provides sufficient load-carrying capacity of the biomaterial. Here, a mathematical modeling approach for estimating the macroscopic load that initiates failure of a hierarchically organized, granular, hydroxyapatite-based biomaterial is presented. For this purpose, a micromechanics model is developed for downscaling macroscopically prescribed stress (or strain) states to the level of the needle-shaped hydroxyapatite crystals. Presuming that the biomaterial fails due to the quasi-brittle failure of the most unfavorably stressed hydroxyapatite needle, the downscaled stress tensors are fed into a suitable, Mohr-Coulomb-type failure criterion, based on which the macroscopic failure load is deduced. The change of the biomaterial’s composition in response to placing it in physiological solution, caused by growth of new bone tissue on the granules’s surfaces, on the one hand, and by resorption of the hydroxyapatite crystals, on the other hand, is taken into account by means of suitable evolution laws. Numerical studies show how the macroscopic load-carrying capacity of the biomaterial is influenced by its design parameters. The presented modeling approach could prove beneficial for the design process of the studied biomaterials (as well as similarly composed biomaterials), particularly in terms of optimizing its mechanical performance.

type failure criterion, based on which the macroscopic failure load is deduced. The change of the biomaterial's composition in response to placing it in physiological solution, caused by growth of new bone tissue on the granules's surfaces, on the one hand, and by resorption of the hydroxyapatite crystals, on the other hand, is taken into account by means of suitable evolution laws. Numerical studies show how the macroscopic loadcarrying capacity of the biomaterial is influenced by its design parameters. The presented modeling approach could prove beneficial for the design process of the studied biomaterials (as well as similarly composed biomaterials), particularly in terms of optimizing its mechanical performance. Keywords Continuum micromechanics · Elastic limit · Multiscale modeling · Bone ingrowth · Tissue engineering Angle defining the orientation of the spherical coordinate system (e r , e ϑ , e ϕ ) ψ Angle defining orientation of vector n

Introduction
The field of bone tissue engineering aims at the reinforcing or even replacing diseased (or for other reasons malfunctioning) bone tissue by scaffold structures that are specifically engineered, for blending in the targeted physiological environment, i.e. the immediate vicinity of bone tissue, as well as possible (Burg et al. 2000;Reichert and Hutmacher 2011). From a mechanical point of view, careful tuning of such scaffold structures (and of the materials which they are made of) is called for because contradictory requirements must be brought in line-scaffold structures must be stiff enough to sustain all relevant mechanical load cases, but also soft enough to facilitate, through mechanobiological couplings (Klein-Nulend et al. 2005;Porter et al. 2009;Velasco et al. 2015), the integration into their bony environment. In this regard, two mechanical properties are of particular interest, both on material and structural levels: the stiffness, governing the elastic deformation behavior and therefore the forces attracted by the involved macro-and microstructures; as well as the strength, indicating the stress level that induces material failure.
In the present paper, we study one specific scaffold material that has been developed as bone replacement material with the human mandible as targeted application area (Komlev et al. 2002(Komlev et al. , 2003. This biomaterial is produced in form of porous, pre-cracked granules, composed of hydroxyapatite as main constituent, but also including various kinds of pore spaces of distinctively different characteristic lengths. After exposing this biomaterial to the targeted physiological environment prevailing in the immediate vicinity of mandibular bone tissue, two mechanisms are triggered, causing a progressing change of the material's composition over time. On the one hand, bone tissue grows on the granule surfaces, while, on the other hand, concurrently the hydroxyapatite crystals are resorbed-in the long run, the scaffold material merges with the surrounding bone tissue. In a first approach to analyzing their mechanical behavior, these granules underwent micro-computed tomography (μCT), and the resulting scans served as basis for combined Finite Element/micromechanicsbased simulations (Dejaco et al. 2012(Dejaco et al. , 2016. Here, as a (computationally more efficient) complement, we present a three-step, fully continuum micromechanicsbased macro-to-meso-to-micro (stress and strain) downscaling scheme, linking in the end the quasibrittle failure of single micrometer-or sub-micrometersized hydroxyapatite crystal needles to the overall strength of both millimeter-sized biomaterial scaffolds and composites comprising biomaterial scaffold and bone tissue, respectively. For this purpose, a number of homogenization concepts are adapted, extended, and combined, considering the pioneering contributions of Eshelby (1957), Hill (1963Hill ( , 1965, Laws (1977Laws ( , 1985, Hervé and Zaoui (1993); and also considering more recent contributions of Deudé et al. (2002), Dormieux et al. (2004), Fritsch et al. (2006, Bertrand andHellmich (2009). Following Fritsch et al. (2009a, b), we feed the stress of the most unfavorably loaded hydroxyapatite needle into a suitable, Mohr-Coulombtype failure criterion, and deduce then therefrom the corresponding ultimate macroscopic load bearable by the aforementioned granular, hydroxyapatite-based biomaterial (optionally containing ingrown bone tissue).
After introducing the fundamental modeling concept, together with the chosen model representation of the studied biomaterial, see Sect. 2, a mathematical model for downscaling of the mechanical loading, from the macroscopic to the hydroxyapatite needle scale, is presented, see Sect. 3. Then, a suitable failure criterion is elaborated in Sect. 4.1, and numerical studies show how the macroscopic mechanical loading inducing single hydroxyapatite needle failure changes with varying biomaterial composition. In order to simulate bone regeneration (which occurs after having placed the biomaterial in the targeted physiological environment), involving bone growth and scaffold resorption, suitable evolution laws are introduced, and the effects of different material input parameters on the model-predicted development of the load-carrying capacity over time are studied, see Sect. 5. A brief discussion closes the paper, see Sect. 6.

Characterization of the multi-porous
hydroxyapatite tissue engineering scaffold material The biomaterial investigated in this paper is a granular scaffold material, with the granules composed of carbonate-containing hydroxyapatite, the chemical composition of which reads as Ca 10 (PO 4 ) 6 (OH) 1.9 (CO 3 ) 0.05 . This material is produced based on the effect of immiscible liquids (Komlev et al. 2002(Komlev et al. , 2003, giving access to granules with diameters ranging from 50 to 2000 µm, with the technically relevant granule diameter being approximately 1800 µm (Dejaco et al. 2012). Several morphological features of these granules can be observed, see the column on the left-hand side of Fig. 1. Firstly, the granules contain pores of two different characteristic lengths: small pores, with a characteristic length ranging from less than one to several micrometers (Dejaco et al. 2016)-these pores are termed "micropores" hereafter; and large pores, with a characteristic length of several hundred micrometersthese pores are termed "mesopores" hereafter. A composite of randomly oriented hydroxyapatite crystals and the micropores constitutes the "base material" of the granules. Increasing the observation scale by several orders of magnitude, one can discern, besides the mesopores, cracks pervading the granule body. Finally, the scaffold material is made up of the above described granules, with pore space in-between-due to the characteristic length of these pores, which is approximately equal to the granule diameter, they are termed "macropores" in the remainder of this paper. the depicted images have been acquired by means of scanning electron microscopy (hierarchical level I) and μCT imaging techniques (hierarchical levels II and III) 2.2 Fundamentals of continuum micromechanics: the representative volume element A method particularly well suited for modeling the mechanical behavior of the material described in Sect. 2.1 is continuum micromechanics (Hill 1963;Zaoui 1997Zaoui , 2002, where a material is understood as a micro-heterogeneous body filling a macro-homogeneous representative volume element (RVE) with characteristic length , d, d standing for the characteristic length of inhomogeneities within the RVE, and L, L standing for the characteristic lengths of geometry or loading of a structure built up by the mater-ial defined on the RVE. It should be noted the aforementioned requirements of "much larger" ( ) and "much smaller" ( ), respectively, have been shown to be already satisfied if the respective characteristic lengths are separated by a factor of two to three and five to ten, respectively (Drugan and Willis 1996;Kohlhauser and Hellmich 2013).
In general, the microstructure within an RVE is too complicated to be described in complete detail. Therefore, quasi-homogeneous subdomains with known physical properties (such as volume fractions and mechanical properties) are reasonably chosen. They are called material phases. The homogenized (upscaled) behavior of the material on the observation scale of the RVE, i.e. the relation between homogeneous deformations acting on the boundary of the RVE and resulting macroscopic (average) stresses, can then be estimated from the mechanical behavior of the material phases, their volume fractions within the RVE, their characteristic shapes, and their interactions. If a single material phase is micro-heterogeneous itself, its mechanical behavior can be estimated by introduction of RVEs within this phase, with characteristic lengths 1 ≤ d, comprising again inhomo-geneities with characteristic length d 1 1 , and so on. Such an approach is referred to as multi-step homogenization and provides, eventually, access to "universal" phase properties at sufficiently low observation scales (Fritsch and Hellmich 2007).

Micromechanical modeling
Having in mind the concept of "separation of scales", as introduced in Sect. 2.2, the following three-level micromechanical representation emerges for the biomaterial under investigation: On hierarchical level I, a microporous, overall isotropic, hydroxyapatite polycrystal is composed of spherical micropores (with volume fraction φ polyHA micro ), which interact mutually with randomly oriented cylindrical hydroxyapatite crystals (with volume fraction f polyHA HA = 1 − φ polyHA micro ). Typically, the microporosity amounts to φ polyHA micro = 0.445 (Dejaco et al. 2012). The characteristic length of the polycrystalline RVE I is in the order of 10 µm, see the bottom row of Fig. 1, with a scanning electron micrograph of the granule nanostructure on the left-hand side and the corresponding RVE I on the right-hand side. In terms of stiffness upscaling, the mutual mechanical interaction of all phases within RVE I calls for a self-consistent homogenization scheme, as introduced by (Fritsch et al. 2006), giving access to the stiffness tensor of the microporous hydroxyapatite polycrystal, C polyHA , based on the composition and morphology of RVE I, as well as on the stiffness tensors of the hydroxyapatite crystals, C HA , and of the micropores, C microφ .
On hierarchical level II, penny-shaped cracks (with vanishing volume fraction) and spherical mesopores (with volume fraction φ gran meso ) are embedded in the polycrystal matrix with properties arising from the structure of RVE I, this matrix filling within RVE II the volume fraction f gran polyHA = 1 − φ gran meso . Typically, the mesoporosity comes to φ gran meso = 0.189 (Dejaco et al. 2012). The characteristic length of RVE II is in the order of 1 mm, see the middle row in Fig. 1, with a microcomputed tomography (μCT) image of the microstructure within a granule on the left-hand side and the corresponding RVE II on the the right-hand side. The distinctive matrix-inclusion morphology of RVE II-i.e. cracks and mesopores can be considered as inclusions embedded in the hydroxyapatite polycrystal matrixsuggests the use of a Mori-Tanaka-type homogenization scheme (Mori and Tanaka 1973;Benveniste 1987) for stiffness homogenization; mathematical treatment of the penny-shaped cracks has been dealt with by Deudé et al. (2002), Dormieux et al. (2004). The stiffness tensor of the pre-cracked, mesoporous granule material, C gran is then governed by the composition and morphology of RVE II, as well as by the stiffness tensors of the hydroxyapatite polycrystal matrix, C polyHA , accessible from stiffness homogenization across RVE I, and of the mesopores, C mesoφ , and by the density of cracks, quantified by the so-called crack density parameter e (Budianksy and O'Connell 1976).
On hierarchical level III, a macroporous conglomerate material consisting of mesoporous, cracked hydroxapatite granules and newly grown bone tissue emerges, see the top of Fig. 1: granules with the stiffness of RVE II described above and filling volume fraction f congl gran , are surrounded by layers of newly grown bone tissue, with volume fraction f congl bone and stiffness derived from the ultrasonic tests of Ashman and van Buskirk (1987). These coated spherical elements are assembled, in mutual contact, to a granular conglomerate with macropores, with volume fraction φ congl macro , in-between. At the time of granule implantation, no bone tissue has been formed yet, and this initial configuration is characterized by f congl bone = 0. For estimating the macroscopic stiffness tensor of the bone-scaffold conglomerate, C congl , the homogenization approach for an nlayered spherical inclusion proposed by Hervé and Zaoui (1993) is specialized for n = 1 (relating to bone tissue), adapted for the case that the stiffness of this layer is transversally isotropic, see (Bertrand and Hellmich 2009), and further combined with a selfconsistent homogenization scheme, in order to account for mutually interacting coated spheres with porous space in-between-in absence of any explicit "matrix phase". This homogenization step is thus based on the composition and morphology of RVE III, as well as on the stiffness tensors of the granule material, C gran , accessible from stiffness homogenization across RVE II, of the bone tissue, C bone , and of the macropores, C macroφ ; the underlying mathematical framework is described at length in (Scheiner et al. 2016).

Downscaling of stresses from macro-to microscale
The stiffness homogenization scheme for the herein investigated biomaterial scaffold for mandibular regeneration, described in Sect. 2.3, constitutes the basis for stress downscaling. For conciseness, the derivation of this (mathematically extensive) model is not repeated here in minute detail; instead, we focus on elaborating how stress downscaling is achieved, while assuming that the homogenized stiffness tensors on all observation scales, compare Fig. 1, namely the stiffness tensor of the nanoporous hydroxyapatite matrix, C polyHA , the stiffness tensor of the pre-cracked, mesoporous granule material, C gran , and the stiffness tensor of the bone-scaffold conglomerate, C congl , are known from the respective homogenization steps. In particular, we first clarify the mechanical input parameters the employed downscaling approach is based on (see Sect. 3.1), and present then briefly the three downscaling steps in terms of the underlying mathematical framework: 1. From the macroporous bone-scaffold conglomerate to the pre-cracked and mesoporous granules (see Sect. 3.2); 2. From the pre-cracked and mesoporous granules to the microporous, polycrystalline hydroxyapatite matrix (see Sect. 3.3); and 3. From the microporous, polycrystalline hydroxyapatite matrix to single, arbitrarily oriented hydroxyapatite crystal needles (see Sect. 3.4).

Definition of mechanical input parameters
As for the underlying main elementary constituent, i.e. hydroxyapatite, the respective stiffness tensor, C HA , is defined via the bulk modulus, k HA , and the shear modulus, μ HA , C HA = 3k HA K + 2μ HA J, with K being the volumetric part of the fourth-order unit tensor I, and J the corresponding deviatoric part, K + J = I.
Numerical values for k HA and μ HA are found based on the experiments performed by Katz and co-workers (Katz and Ukraincik 1971;Gilmore and Katz 1982), who revealed the Young's modulus and Poisson's ratio of hydroxyapatite, E HA = 114 GPa and ν HA = 0.27, see also (Hellmich and Ulm 2002;Hellmich et al. 2004). Through standard relations of continuum (Mang and Hofstetter 2000), one finally obtains k HA = 82.61 GPa and μ HA = 44.88 GPa. Furthermore, all pore spaces are assumed to be drained at all times, thus For defining the stiffness tensor of newly formed bone tissue, C bone , we consider the orthotropic stiffness tensor determined for mandibular bone by means of ultrasound measurements by Ashman and van Buskirk (1987), and follow then the strategy described in (Bertrand and Hellmich 2009) for converting the ultrasound-based, anisotropic stiffness tensor into the transversally isotropic stiffness tensor related to the bone tissue growing on granules, see also (Scheiner et al. 2016) for details, yielding with superscript (e r , e ϑ , e ϕ ) indicating that this stiffness tensor is expressed in a spherical coordinate system, owing to the fact that in the particular case, bone tissue is added in form of a concentric shell on top of the spherical granules.
3.2 From the macroporous bone-scaffold conglomerate to the pre-cracked, mesoporous granules ( Fig. 1, hierarchical level III) Based on the pioneering work of Hervé and Zaoui (1993), the aforementioned stiffness homogenization scheme for granular hydroxyapatite-based biomaterials, see Sect. 2.3, gives, on the one hand, access to the volumetric stress tensor of the granule material, ε congl,vol gran , in response to a macroscopically applied volumetric strain tensor E vol congl , and, on the other hand, to the deviatoric stress tensor of the granule material, ε congl,dev gran , in response to a macroscopic strain tensor representing pure shear, E shear congl , where k gran is the bulk modulus, μ gran the shear modulus, and ν gran the Poisson's ratio of the granule material. Furthermore, Eqs. (2) and (3) are governed by a number of material properties, namely k gran,1 , μ gran,1 , μ gran,2 , and D k . These material properties are functions of the conglomerate's composition (quantified by volume fractions f congl gran and f congl bone ), as well as the stiffness tensors of the conglomerate consisting of scaffold material and bone tissue, C congl , of the granule material, C gran , of the added bone matrix, C bone , and of the macropores, C macroφ . The somewhat unusual form of the downscaling relations given by Eqs. (2) and (3) as well as of the lengthy mathematical expressions for k gran,1 , μ gran,1 , μ gran,2 , and D k , given in detail in (Scheiner et al. 2016), result from the nature of the employed coated inclusionproblem of Hervé and Zaoui (1993), where, in contrast to the classical Eshelby-problem, the stresses and strains are not homogeneous throughout the inclusion.
Considering now that any strain (and, of course, also stress) tensor can be decomposed into volumetric and pure shear strain tensors, Eqs. (2) and (3) can be utilized for downscaling a general macroscopic strain tensor, E congl to the corresponding stress tensor experienced by the granule material, where tr is the trace operator, tr E congl = E congl,11 + E congl,22 + E congl,33 , and I is the second-order unit tensor. Given that the underlying constitutive law on the observation scale of the bone-scaffold conglomerate is linear elastic, congl = C congl : E congl , Eq.(4) Fig. 2 Image of the studied bone-scaffold conglomerate inserted into the upper jaw, based on which the macroscopic mechanical boundary conditions of the conglomerate consisting of bone tissue-coated granules are defined: load application occurs in direction e 1 (resembling mastication), deformation in direction e 3 is presumably prevented, the material is stress-free in direction e 2 , the conglomerate material is furthermore free of shear stresses and strains can be straightforwardly applied for also downscaling macroscopic stress tensors to the granule material level, through substituting E congl by (C congl ) −1 : congl . In the following, we want to exemplify stress downscaling, quantitatively, through prescribing macroscopically the predominant loading type expected for the studied biomaterial in physiological conditions. Typically, the biomaterial investigated in this paper is inserted from the buccal (exterior) mandible surface, and the granules-filled bone defect is covered afterwards by means of a bioresorbable membrane, see Fig. 2 (where the membrane has not yet been put in place). Considering that the macroscopic loading is usually prescribed in terms of stress component congl,11 < 0 (relating to mastication), and that the further mechanical boundary conditions for the bonescaffold conglomerate are E congl,11 < 0, E congl,22 > 0, E congl,33 = E congl,12 = E congl,13 = E congl,33 = 0, congl,33 < 0, and congl,22 = congl,12 = congl,13 = 0 gives access, via the linear elastic constitutive law, to the corresponding macroscopic stress tensor: In the following, the macroscopic stress tensor component in direction of base vector e 1 is set to congl,11 = −10 MPa, while Eq. (5) gives access to the second nonzero component of the macroscopic stress tensor. Furthermore, the crack density parameter for quantifying the occurrence of cracks on hierarchical level II is set to e = 10 (representing a low to moderate crack density), see Sect. 3.3 for further details. The downscaling relation given by Eq. (4) provides then the corresponding stress tensor experienced by the granule material, see Fig. 3, where also the effect of increasing bone volume fraction is illustrated. Notably, the macrocopically applied biaxial loading relates to a triaxial loading on one observation scale below.
3.3 From the pre-cracked, mesoporous granules to the microporous hydroxyapatite polycrystal (Fig. 1, hierarchical level II) For this downscaling step we can make use of classical strain and stress downscaling as defined in the framework of continuum micromechanics. As derived elsewhere, see e.g. (Zaoui 2002;Dormieux et al. 2006), continuum micromechanics allows to downscale strain tensors through the so-called concentration (or localization) tensors: for the present case, this implies that the strain tensor of the microporous hydroxyapatite polycrystal, ε gran polyHA , is related to the strain tensor of the granule material, ε where A gran polyHA is the strain concentration tensor of the microporous hydroxyapatite polycrystal. The latter is estimated by means of Eshelby's matrix-inclusion problem (Eshelby 1957); considering the matrixinclusion-type morphology discernible on hierarchical level II, compare Fig. 1, A gran polyHA is defined by (Deudé et al. 2002;Dormieux et al. 2004) where P polyHA sph is the fourth-order Hill tensor related to spherical inclusions embedded in a matrix exhibiting a stiffness tensor C polyHA , see (Eshelby 1957;Zaoui 2002) for how P polyHA sph is computed. Hence, in contrast to the situation in Sect. 3.2, we here employ the Eshelby's classical matrix-inclusion problem (Eshelby 1957). Finally, e is the so-called crack density parameter (Budianksy and O'Connell 1976), e = N (r crack ) 3 , with N as the number of cracks per volume, and r crack as the (average) crack radius, and Q is a tensor defined via the Poisson's ratio of the microporous hydroxyapatite polycrystal, ν polyHA , through ) with the stress concentration tensor of the microporous hydroxyapatite matrix, B gran polyHA , following from the strain concentration tensor A Downscaling of stresses according to Eq. (9) leads to further magnification of stresses; see Fig. 4 for the nonzero stress tensor components of the microporous granule material, obtained through application of Eq. (9) to the stresses depicted in Fig. 3.

From the microporous hydroxyapatite polycrystal
to the single hydroxyapatite crystal (Fig. 1, hierarchical level III) The hydroxyapatite crystals making up the the microporous hydroxyapatite polycrystal are oriented arbitrarily, in all space directions, defined by in a spherical coordinate system by Euler angles ϑ and ϕ, see with σ gran polyHA following from Eq. (9).
Feeding the composition-dependent stresses depicted in Fig. 4 into Eq. (11) reveals substantial (hydroxyapatite needle orientation-dependent) stress magnification effects, see Fig. 6. In particular, Fig. 6a-f show the variations between the different components of the hydroxyapatite stress tensor, for φ congl macro = 0.4 and f congl bone = 0. In order to also highlight the extreme values of the stress tensor components as functions of the scaffold composition, Fig. 6g-i show the minimum and maximum values of all stress tensor components over all hydroxyapatite needle orientations, see also Table 1 for an overview of the observed maximum and minimum stress tensor components. Notably, minimum stresses stresses (i.e. compressive stresses with maximum magnitude) always occur in direction e 1 , while maximum stresses (i.e. tensile stresses with maximum magnitude) always occur in direction e 2 .

Failure criterion suitable for hydroxyapatite needles
Based on the three-step scheme presented in Sect. 3, a macroscopically applied mechanical loading, prescribed in terms of macroscopic strains E congl or macroscopic stresses congl , were downscaled to the corresponding stress state experienced by a single, arbitrarily oriented hydroxyapatite needle, σ polyHA HA (ϑ, ϕ). From (Fritsch et al. 2009a, b), we adopt that hydroxyapatite needle failure is governed by the normal stress in needle direction, and by the shear stress in planes orthogonal to the needle direction, σ HA,rn (ϑ, ϕ, ψ) = e r (ϑ, ϕ) · σ polyHA HA (ϑ, ϕ) · n(ϑ, ϕ, ψ) .
In Eqs. (14) and (15), e r (ϑ, ϕ) denotes the vector defining the direction of a particular needle, i.e. the base vector of the employed spherical coordinate system in radial direction, while vector n(ϑ, ϕ, ψ) denotes the direction orthogonal to e r (ϑ, ϕ), additionally governed by angle ψ, compare Fig. 5. As for definition of vector n(ϑ, ϕ, ψ), we consider that the needle orientationdependent base vectors (e r , e ϑ , e ϕ ) are defined in the Cartesian base system (e 1 , e 2 , e 3 ) as e r = sin ϑ cos ϕ, sin ϑ sin ϕ, cos ϑ T , e ϑ = cos ϑ cos ϕ, cos ϑ sin ϕ, − sin ϑ T , and that the transformation tensor Q from base system (e 1 , e 2 , e 3 ) to base system (e r , e ϑ , e ϕ ) reads as Then, vector n(ϑ, ϕ, ψ), dependent on angles ϑ, ϕ, and ψ, and expressed in a Cartesian base frame, follows as As furthermore proposed in (Fritsch et al. 2009a, b), the failure criterion for a single hydroxyapatite needles takes into account both tensile strength σ ult,t HA and shear strength σ ult,s HA . The two strength values are accessible through analyzing the experiments of Akao et al. (1981) and Shareef et al. (1993), revealing σ ult,t HA = 52.2 MPa and σ ult,s HA = 80.3 MPa, see (Fritsch et al. 2009a). Mathematically, the failure surface related to the hydroxyapatite needles reads as Eq. (19) takes into account that for each needle orientation the angle ψ inducing the maximum tangential stress must be found, as well as the needle orientation inducing the maximum value for the combination of normal and tangential stresses. For this purpose, angles ϑ, ϕ, and ψ are varied between ϑ = 0 . . . π, ϕ = 0 . . . 2π , and ψ = 0 . . . 2π . Substituting into Eq. (19) the relation between σ polyHA HA and congl , according to the downscaling scheme elaborated in Sect. 3 allows Choice of initial value for congl,11 .
-If f HA ( congl ) = 0, then the load iteration is completed, and the current magnitude for congl,11 induces failure of the material.
to alternatively express the failure criterion in terms of macroscopic stresses, 4.2 Computation of composition-dependent macroscopic loading inducing hydroxyapatite needle-failure The failure criterion presented in Sect. 4.1, for estimating the macroscopic loading that leads to quasi-brittle failure of the most unfavorably stressed hydroxyapatite needle, is evaluated in an iterative manner, see Table 2.
Let us e.g. consider, for the sake of demonstration, a biomaterial configuration defined by φ congl macro = 0.4 and f congl bone = 0. Then, the iterative approach sketched in Table 2 reveals that a macroscopic stress tensor with non-zero components congl,11 = −20.31 MPa and congl,33 = −3.09 MPa is related to f HA ( congl ) = 0, thus inducing failure of the most unfavorably stressed hydroxyapatite needle. The stress state of the latter is illustrated in Fig. 7a and b, in terms of the normal and maximum shear stress components, as functions of the needle orientation. The corresponding values of the failure function f HA ( congl ), as obtained through insertion of the orientation-dependent normal and maximum shear stress components into the failure criterion, is again a function of the needle orientation, as depicted in Fig. 7(c). In this figure, the needle orientations that are actually evoking the most unfavorable stress states for the aforementioned macroscopic load-ing inducing f HA ( congl ) = 0, are indicated, namely (ϕ = 1.466, ϑ = 1.676) and (ϕ = 4.608, ϑ = 1.466), respectively; or (ϕ = 84.00 • , ϑ = 96.03 • ) and (ϕ = 264.03 • , ϑ = 84.00 • ), respectively.
Furthermore, the magnitude of the macroscopic, material failure-inducing loading increases with increasing bone volume fraction ( f congl bone ), as well as with decreasing volume fraction of the macropores (φ congl macro ), see Fig. 7d.

Development of failure-inducing macroscopic stresses during bone regeneration
Once immersed into its targeted physiological environment, i.e. the immediate vicinity of (mandibular) bone tissue, the studied biomaterial undergoes compositional changes, due to two distinct mechanisms. On the one hand, the granules become coated with a layer of newly forming bone tissue whose thickness is growing with time. Presuming that the growth of bone tissue occurs linearly over time (Cancedda et al. 2007), the following bone volume fraction evolution law can be deduced: where t is the time after placing the scaffold in a bony environment, k growth is the bone growth rate, and r gran is the radius of the granules. For the presently studied material, the typical bone growth rate is k growth =  (Scheiner et al. 2016). The timedependent bone tissue volume fraction defined by Eq. (21) enters, as given in great detail in (Scheiner et al. 2016), the material functions k gran,1 , μ gran,1 , μ gran,2 , and D k appearing in Eqs. (2)(3)(4). In this way, Eq. (21) induces time-dependency in the strength of the investigated conglomerate material. In the same way, Eq. (21) induces time-dependencies in the homogenized stiffness expressions relating to the RVEs depicted in Fig. 1, as given in (Scheiner et al. 2016).
On the other hand, the contact of the bone-coated granules with physiological solution also leads to resorption of hydroxyapatite. As is known from analysis of X-ray microtomography scans, this resorption process occurs via growth of the micropores (Czenek et al. 2014), with φ polyHA micro,0 as the microporosity before resorption sets in, and k res as scaffold resorption rate. Notably, φ polyHA micro enters stress and strain downscaling at Eq. (13) of this paper, thus influencing strain concentration tensorA polyHA HA (ϑ, ϕ), and therefore again the timedependent strength and stiffness properties of the con-  (Scheiner et al. 2016).
In the following, a series of parametric studies are presented, for elucidating the effects of model parameter variations on the predicted development of the macroscopic failure-inducing stress tensor component in direction e 1 , ult congl,11 . For the sake of clarity, all model input parameters (and the variations considered for the numerical studies presented in this paper) are summarized in Table 3. The accordingly performed simulations highlight strong sensitivites of the macroscopic loading that causes failure of the hydroxyapatite needles, on the considered parameter variations, see Fig. 8. As regards the granule radius, it is striking that decreasing the radius initially implies accelerating the bone ingrowth-related increase of the loadcarrying capacity, see Fig. 8a. However, the smaller the granule radius the sooner the complete pore space is filled with bone tissue, which is indicated by a peak of ult congl,11 (the actual magnitude of which depends on the macroporosity), followed by prolonged decrease of the load-carrying capacity, caused by scaffold resorption. Furthermore, it is interesting to note that the granule radius influences strongly the development of ult congl,11 directly after placing the granules in the targeted physiological environment (initiating scaffold resorption and bone formation), while the actual value of ult congl,11 at a later point of time is only governed by the macroporosity. A high macroporosity supports the strength development at mature bone regeneration states since it provides substantial space for newly added, load-carrying bone tissue; but it compromises strength before onset of bone regeneration as well as at early regeneration states as then a comparably low amount of bone-scaffold conglomerate needs to carry all (or most) of the loading.
The scaffold resorption rate, in turn, governs both the load-carrying capacity of the bone-coated granules when the macroporosity is completely filled with new bone tissue, as well as the long-term development afterwards. Zero resorption implies that the value of ult congl,11 is maintained at a constant (maximum) level (related to complete filling of the pore space by bone matrix), whereas a non-zero resorption rate causes a long-term decrease of ult congl,11 after reaching the aforementioned maximum value, see Fig. 8b. It should be noted that this long-term decrease is caused by assuming that bone growth merely occurs on the outer surface of the (bone-covered) granules. However, in reality, it can be assumed that after substantial resorption of the hydroxyapatite crystals the physiological solution enters the micro-and mesopore spaces, leading eventually also to bone formation within the granule body. The omission of this potential additional bone formation effect can be deemed as limitation of our model.
Similar to the granule radius, the bone formation rate only influences the time span until the whole macropore space is filled with bone tissue; the long-term development of the load-carrying capacity is unaffected, see Fig. 8c.
Remarkable effects are revealed when varying the crack density parameter e, see Fig. 8d. In particular, decreasing the crack density parameter leads to a significantly increased stiffness of the granule material. Considering that in a composite material (such as the material studied in this paper) stiffer constituents attract larger fractions of a macroscopically applied loading than softer constituents, which attract lower fractions of the macroscopic loading, a granule material containing less cracks hence transfers higher stresses to the hydroxyapatite needles than those with more cracks. This eventually implies that increasing the crack density in the granule material leads to an increased load-carrying capacity. This possibly counterintuitive conclusion straightforwardly suggests that future extensions of the here presented model should comprise formulation of a failure criterion related to the bone tissue as well, in order to improve the significance of the model-predicted load-carrying capacity. φ congl macro = 0.3, r gran = 300 μm φ congl macro = 0.4, r gran = 300 μm φ congl macro = 0.5, r gran = 300 μm φ congl macro = 0.3, r gran = 500 μm φ congl macro = 0.4, r gran = 500 μm φ congl macro = 0.5, r gran = 500 μm φ congl macro = 0.3, r gran = 1000 μm φ congl macro = 0.4, r gran = 1000 μm φ congl macro = 0.5, r gran = 1000 μm In this paper, a continuum micromechanics-based model was presented for estimating the macroscopic loading acting onto a hydroxyapatite-based granular biomaterial, developed for application as bone replacement material (in mandibular bone), that leads to quasibrittle failure of the material's main constituent, i.e. hydroxyapatite crystals. The parametric studies presented in Sect. 5 show how the load-carrying capacity of the studied biomaterial develops over time once placed in the targeted physiological environment, i.e. the immediate vicinity of mandibular bone, consider-ing the growth of new bone tissue on the surface of the scaffold material, and resorption of the hydroxyapatite needles. Thereby, main emphasis was on highlighting the influence of specific design parameters of the production process, e.g. the exact chemical composition of the biomaterial might influence the rates of bone ingrowth and scaffold resorption, the crack density may be related to the production process, and the macroporosity can be tuned based on the packing density of the granules. From a practical point of view, the presented modeling approach allows to determine from when onwards a particular area of the mandible including an implant composed of the studied biomaterial can be used for mastication if the mechanical loading acting onto this mandibular region is approximately known, e.g. from Finite Element simulations (Korioth et al. 1992;Meijer et al. 1993;Choi et al. 2005;Hellmich et al. 2008;Bevilacqua et al. 2011).
However, the simulation results also point out two model restrictions. On the one hand, the kind of counterintuitive observation was made that a severely cracked granule material implies that the bone tissue growing on the granule surface attracts most of themacroscopically applied stress. Thus, our model suggests that increasing the crack density leads to an increasing load-carrying capacity, owing to the fact that the employed failure criterion considers only the most unfavorably stressed hydroxyapatite needle contained in the granules, neglecting however the stress experienced by the newly formed bone tissue. On the other hand, our model does not consider that the dissolution of hydroxyapatite needles would eventually lead to morphological changes in the microporous hydroxyapatite matrix (hierarchical level I in compare Fig. 1), implying that physiological solution could enter the micropore space, facilitating there bone tissue formation. Given that the described model inadequacies become relevant only after a certain (not yet quantifiable) time span, but not directly after scaffold implantation, leads to the conclusion that our model is particularly accurate for early-age bone-scaffold conglomerates (with respect to the time instant when the granules are placed into the targeted physiological environment), while the prediction accuracy presumably diminishes over time. This restriction constitutes the basis for reasonable future research directions.
Finally, it is also important to discuss the relevance of traditional fracture mechanics approaches, typically focusing on the prediction of crack propagation, thus assuming the existence of an initial crack (Müller et al. 2002;Näser et al. 2007;Kolednik et al. 2010;Ott et al. 2010), in the context of the material studied in this paper. Actually, it seems to be a worthwhile subject of future research actitivities to extend the micromechanics-based assessment of specific, microscopically sized material constituents (as demonstrated in this paper) towards traditional fracture mechanics, see e.g. (Pichler et al. 2007;Pichler and Dormieux 2009b, a). Such extension would be particularly relevant for mature bone-scaffold conglomerates containing already a substantial amount of bone tissue, allowing for studying the effects of crack emergence and propagation in the bone tissue-given that a respective failure criterion has been formulated, see e.g. (Fritsch et al. 2009b)-as well as of biologically driven crack healing.