Procedure for Estimating the Elastic Properties of Inhomogeneous Microstructures Such as Bone Tissue

Fig. 1

Illustration of the RVE generating algorithm. Material is agglomerated randomly around existing material

While the volume fraction is directly regulated, the RVE structure develops indirectly by the number of initial cells. A rather rough cluster with high material agglomerations emerges from a small number of initial cells, whereas a fine dispersed cluster emerges from many initial cells. Figure 2 shows three different RVEs, each with 50 elements per edge and equal volume fraction but varying number of initial cells.

Fig. 2

Agglomerated material of different structures is emerged from varying number of initial cells. Left 0.00625 %, middle 0.2 % and right 6.4 % initial cells, all with 25 % volume fraction in a grid of $50\times 50\times 50$

2.2 Continuum Mechanics Approach

The continuum mechanics approach is used to calculate three-dimensional material deformations. For static considerations the momentum balance of the current configuration is reduced to a time and mass invariant equilibrium that can be expressed by the divergence of the Cauchy stress tensor.

$\begin{aligned} div\varvec{\sigma } = 0 \end{aligned}$

(1)

This expression is under-constrained and additional definitions are required. First of all, the continuity of the field quantities is postulated, meaning that all deformations are physically objective and infinitesimal small material particles are not allowed to penetrate each other or fluctuate. This uniqueness is obtained by the definition of the deformation tensor.

$\begin{aligned} \mathbf{{F}} = \frac{{dx}}{{dX}} \end{aligned}$

(2)

In terms of the physical objectivity, the material behavior of elastic bodies (also denoted as Cauchy elasticity) now demand tensor compatibility of stress and deformation.

$\begin{aligned} {\varvec{\sigma }} = f(\mathbf{{F}}) \end{aligned}$

(3)

Such compatibility is given by the Rivlin-Ericksen theorem for isotropic behavior [11].

$\begin{aligned} \varvec{\sigma } = a_0 + a_1 \mathbf{{F}}^T \mathbf{{F}} + a_2 (\mathbf{{F}}^T \mathbf{{F}})^2 \end{aligned}$

(4)

Thereby, the coefficients

and

are arbitrary scalar functions of the invariants of the deformation tensor. Usually, this dependence is formulated in relation to the strain energy density $\psi$ . Thus, the Green or hyper-elastic material behavior can be defined by a pure scalar function $\psi =f(I_1,I_2,J)$ , with the constitutive law:

$\begin{aligned} \varvec{\sigma } = J^{ - 1} \frac{{d\psi }}{{d\mathbf{{F}}}}\mathbf{{F}}^T \end{aligned}$

(5)

However, for this assumption isotropy is presumed, meaning that not every material can be modeled by this constitutive law. In order to consider anisotropic effects as well, the constitutive law is expressed by using the geometric linearized form of the right Cauchy-Green tensor.

$\begin{aligned} \varvec{\varepsilon } = \frac{1}{2}(\mathbf{{F}}^T \mathbf{{F}} - \mathbf {1}) \end{aligned}$

(6)

This modified tensor can be used to derive a simplified relation between stress and deformation. In analogy to (4) an equalization of $f(F) \sim f(\varepsilon )$ leads to Hooke’s law in continuum mechanics

$\begin{aligned} \varvec{\sigma } = \frac{1}{2}(\lambda tr\varvec{\varepsilon }) + 2G\varvec{\varepsilon } \end{aligned}$

(7)

with the assumption $\sigma (\varepsilon =0)=0$ and the Lame’s constants $\lambda =E\nu /((1+\nu )(1-2\nu ))$ and $G=E/(2(1+\nu ))$ [12].

The physical validity of this equation expires with increasing deformation due to the linearization and does not describe the stress decrease of real material. However, in practical (5) and (7) only differ for deformations higher than 5 % strain.

The vectorization of the stress and strain tensors by use of the Voigt notation leads to the well-known matrix equation:

$\begin{aligned} \left[ \begin{matrix} {\sigma _{11} } \\ {\sigma _{22} } \\ {\sigma _{33} } \\ {\sigma _{12} } \\ {\sigma _{13} } \\ {\sigma _{23} } \\ \end{matrix} \right] = \frac{E}{(1+\nu )(1-2\nu )} \left[ \begin{matrix} {1 - \nu } &{} \nu &{} \nu &{} 0 &{} 0 &{} 0 \\ \nu &{} {1 - \nu } &{} \nu &{} 0 &{} 0 &{} 0 \\ \nu &{} \nu &{} {1 - \nu } &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} \frac{1 - 2\nu }{2} &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} \frac{1 - 2\nu }{2} &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} \frac{1 - 2\nu }{2} \\ \end{matrix} \right] \left[ \begin{matrix} {\varepsilon _{11} } \\ {\varepsilon _{22} } \\ {\varepsilon _{33} } \\ {\varepsilon _{12} } \\ {\varepsilon _{13} } \\ {\varepsilon _{23} } \\ \end{matrix} \right] \end{aligned}$

(8)

This linearized equation still describes pure isotropic material behavior. It can be converted into the generalized Hooke’s law by a phenomenological motivated consideration, so that in principle all 36 coefficients can be chosen independently.

$\begin{aligned} \left[ {\begin{array}{*{20}c} {\sigma _{11} } \\ {\sigma _{22} } \\ {\sigma _{33} } \\ {\sigma _{12} } \\ {\sigma _{13} } \\ {\sigma _{23} } \\ \end{array}} \right] = \left[ {\begin{array}{*{20}c} {C_{11} } &{} {C_{12} } &{} {C_{13} } &{} {C_{14} } &{} {C_{15} } &{} {C_{16} } \\ {C_{21} } &{} {C_{22} } &{} {C_{23} } &{} {C_{24} } &{} {C_{25} } &{} {C_{26} } \\ {C_{31} } &{} {C_{32} } &{} {C_{33} } &{} {C_{34} } &{} {C_{35} } &{} {C_{36} } \\ {C_{41} } &{} {C_{42} } &{} {C_{43} } &{} {C_{44} } &{} {C_{45} } &{} {C_{46} } \\ {C_{51} } &{} {C_{52} } &{} {C_{53} } &{} {C_{54} } &{} {C_{55} } &{} {C_{56} } \\ {C_{61} } &{} {C_{62} } &{} {C_{63} } &{} {C_{64} } &{} {C_{65} } &{} {C_{66} } \\ \end{array}} \right] \left[ {\begin{array}{*{20}c} {\varepsilon _{11} } \\ {\varepsilon _{22} } \\ {\varepsilon _{33} } \\ {\varepsilon _{12} } \\ {\varepsilon _{13} } \\ {\varepsilon _{23} } \\ \end{array}} \right] \end{aligned}$

(9)

However, both the strain tensor $\varvec{\varepsilon }$ and the stress tensor $\varvec{\sigma }$ are symmetric, so the stiffness matrix $\varvec{C}$ and the respective compliance matrix $\varvec{N}=\varvec{C}^{-1}$ have to be symmetric as well. Consequently, only 21 independent coefficients remain, which have to provide a positive determinant.

One should consider that this kind of anisotropic modeling is only valid for homogeneous bodies. However, anisotropic behavior in general is caused by inhomogeneity, so this is a crude assumption with very limited validity. For example, a structure causing momentums cannot be homogenized with the constitutive law (9).

However, the estimated stiffness is strongly influenced by the numerical process even for suitable structures.

2.3 Homogenization Approach

In general an analogical homogeneous constitutive law for the underlying inhomogeneous microstructure should fulfill the Hill condition [13].

$$$\begin{aligned} <{\!}\varvec{\sigma }{\!}> <{\!}\varvec{\varepsilon }{\!}> = <{\!}\varvec{\sigma \varepsilon }{\!}> \end{aligned}$$” src=”/wp-content/uploads/2016/10/A324942_1_En_1_Chapter_Equ10.gif”></DIV></DIV><br /> <DIV class=EquationNumber>(10)</DIV></DIV>This condition states that the average strain and stress of the heterogeneous RVE correlate with the quantities of an analogous homogeneous RVE. Strictly speaking, a perfect homogenized equation (<SPAN class=InternalRef><A href=$ 9) of the inhomogeneous microstructure is found if this condition can be fulfilled. This equation is valid for both the analog homogeneous RVE and an arbitrary material point within the continuum. This so-called sequentially homogenization procedure is illustrated in Fig. 3.

Fig. 3

Sequentially homogenization of a heterogeneous microstructure

The boundary value problem is solved by the finite element method (FEM). Any desired structures can be modelled by this procedure. Thereby, the differential equation (1) is converted to its weak form by use of the variational principle.

$\begin{aligned} \int \limits _\Omega {\left( {\sigma _{ij} n_j - t_i } \right) } \delta u_i d\Omega - \int \limits _V {\sigma _{ij,j} } \delta u_i dV = 0 \end{aligned}$

(11)

This equation is numerically solved by discretization [14]. The field quantities are calculated discretely at the nodes of the elements and converge towards the real solution with increasing mesh refinement.

Fig. 4

Left Inhomogeneous strain distribution in cross sections. Right normalized strain distribution over the whole RVE

Figure 4 exemplary shows the inhomogeneous strain distribution in two cross sections. Additionally, the normalized distribution of all strains is plotted. Two peaks arise in the distribution, that represent the different strain levels of both materials. The macroscopic obtained strain, however, lays in between.

2.4 Window Size and Boundary Conditions

Three boundary conditions in terms of displacement or stress can be defined on each surface of the RVE, compare Figs. 5 and 6.

Fig. 5

Left Definitions of RVE surfaces. Right FEM visualization

Fig. 6

Boundary condition sets. Left Pure kinematic constraints. Middle Mixed constraints proposed by Pahr and Zysset [8]. Right Pure stress constraints

The Hill condition demands a compatibility of strain and (!) stress, which would require a definition of six boundary conditions on each surface, meaning that three strain and three stress components have to be applied. Consequently, the Hill condition cannot be completely fulfilled for boundary value problems like the FEM. The obligatory chosen conditions have a strong influence on the stiffness estimation. The structure is estimated too stiff, if pure displacement conditions are applied. Contrary to this, the structure is estimated too soft, if pure stress conditions are applied. In principle, only an apparent solution can be expected. Fortunately, the boundary cells, which cause discontinuity, become less important with increasing RVE size. The apparent solution finally converges towards the effective solution [15]. However, computer resources are limited and size is a very important factor. Doubling the RVE size means eight-times more memory and tripling the RVE size even means 27-times more memory. The stochastic distributed stiffness further reduces the efficiency of common sparse-solver. For example, the calculation of six load cases needs about 2 h on one core and 25 min on six cores (each Xeon E7-4830) for a cube with 50 elements per edge. A simulation with 100 elements per edge runs 27 h on 8 cores and needs 100 GB memory instead of 6.4 GB for 50 elements per edge.

Table 1

KUBC: Kinematic uniform boundary conditions (3 kinematic constraints on each surface)

Load case	Top	Bottom	North	South
	$u_1=x\cdot u_0$	$u_1=x\cdot u_0$	$u_1=x\cdot u_0$	$u_1=x\cdot u_0$
x-tension


xy-shear	$u_2=x\cdot u_0$	$u_2=x \cdot u_0$	$u_2=x\cdot u_0$	$u_2=x\cdot u_0$

Table 2

PMUBC: Periodic mixed uniform boundary conditions (combination of kinematic and stress constraints)

Load case	Top	Bottom	East	West	North	South

x-tension

		Only gold members can continue reading. Log In or Register to continue Share this: Click to share on Twitter (Opens in new window) Click to share on Facebook (Opens in new window) Like this: Like Loading... Related Related posts: Related Infections Test Models for Implant-Associated Inflammation and Infections of a Model of the Electrically Stimulated Cochlea Mechanically Stimulated Fracture Healing Model Using a Finite Element Framework Stay updated, free articles. Join our Telegram channel Tags: Biomedical Technology Oct 21, 2016 \| Posted by admin in BIOCHEMISTRY \| Comments Off on Procedure for Estimating the Elastic Properties of Inhomogeneous Microstructures Such as Bone Tissue Full access? Get Clinical Tree Get Clinical Tree app for offline access Get Clinical Tree app for offline access