Homogenization of heterogeneous tissue scaffold : A comparison of mechanics , asymptotic homogenization , and finite element approach

Actual prediction of the effective mechanical properties of tissue scaffolds is very important for tissue engineering applications. Currently common homogenization methods are based on three available approaches: standard mechanics modeling, homogenization theory, and finite element methods. Each of these methods has advantages and limitations. This paper presents comparisons and applications of these approaches for the prediction of the effective properties of a tissue scaffold. Derivations and formulations of mechanics, homogenization, and finite element approach as they relate to tissue engineering are described. The process for the development of a computational algorithm, finite element implementation, and numerical solution for calculating the effective mechanical properties of porous tissue scaffolds are also given. A comparison of the results based upon these different approaches is presented. Parametric analyses using the homogenization approach to study the effects of different scaffold materials and pore shapes on the properties of the scaffold are conducted, and the results of the analyses are also presented.


INTRODUCTION
Tissue engineering is an emerging interdisciplinary field that applies the principles of biology as well as engineering toward the development of viable substitutes that can restore, maintain, or improve the function of human tissues (Langer 1994(Langer , 1999)).The underlying concept of tissue engineering is that cells can be isolated from a patient, expanded in a culture, and then seeded onto a scaffold prepared from a specific building material (eg extracellular matrix, biodegradable polymer) to form a scaffold/biological three-dimensional tissue construct.The construct can then be grafted into the same patient to function as replacement tissue.Blood vessels attach themselves to the new tissue, the scaffold gradually dissolves, and the newly grown tissue is integrated into its surroundings.Scaffolds for tissue engineering should be designed to meet multiple biological and bioengineering design criteria such as structural integrity, stability, degradability, and transport properties (Hutmacher 2000;Hollister et al 2002).In this regard, porous, three-dimensional (3D) tissue scaffolds play an important role in cell attachment, proliferation, and guidance of new tissue formation.Performance of various functions of the tissue structure depends upon porous scaffold microstructures with specific porosity characteristics.Studies show that pore size and overall porosity affect cellular adhesion, viability, ingrowth, distribution, and the formation of an extracellular matrix (ECM) by specific cell types (Zeltinger et al 2001;Nikolovski and Mooney 2000).The internal architectures of porous hydroxyapatite (HA) implants control the degree of bone regeneration (Magan and Ripamonti 1996;Chang et al 2000), influence the path of bone regeneration (Liu 1996), and determine the mechanical properties of the implants (Daculsi and Passuti 1990;Le Huec et al 1995).Therefore, the ability to quantify the structural heterogeneity and the mechanical properties of a scaffold with a designed microstructure is essential for the tissue construct applications.Available homogenization methods for characterization of the effective properties of porous scaffolds or heterogeneous tissues are primarily based on experimental approaches (Hing et al 1999;Bose et al 2002), standard mechanics modeling (Hill 1963;Hashin 1983), homogenization theory (Suquet 1987;Hollister et al 1991), and finite element method (Beaupr and Hayes 1985;Williams and Lewis 1982).There are some inherent limitations in each of the above approaches.For example, the mechanics modeling and finite element method are both based on the averaged field theory and use a representative volume element (RVE) to predict the effective mechanical properties.The results predicted by these two approaches are sensitive to the size and applied boundary conditions of the used RVE.The asymptotic expansion homogenization theory, on the other hand, may reduce the effect of the RVE boundary and the size as reported in Hollister andKikuchi (1992, 1994) and Aoubiza et al (1996), but the utilization of the theory requires computational algorithm for numerical implementation.So far, there is no available computational algorithm for scaffold characterization and for a general application to tissue engineering.
The objectives of this paper are as follows: (1) compare the available mechanics, homogenization, and finite element approach, including the development of homogenization modeling, its formulations and solution procedures; (2) present our study on developing a computational algorithm for the finite element implementation of asymptotic homogenization theory; and (3) apply the developed algorithm for the prediction of the effective mechanical and structural heterogeneity of tissue scaffolds.The formulations and the solution procedures for available mechanics, homogenization, and finite element approach are introduced in the first section.The computational algorithm developed for the finite element implementation of asymptotic homogenization theory is described in the second section.The applications of the theory, modeling and algorithm, and the results of the predictions, their comparisons, and the parametric analyses are presented in the third section, followed by the summary and conclusions in the last section.
Consider a heterogeneous body ε with a periodic microstructure represented by a porous unit cell (representative volume element) as shown in Figure 1.In general, ε is subjected to a body force f , a traction t at the traction boundary t , and a displacement d at the displacement boundary d .It also contains multiple network channels (or voids) with a pressure p applied on the internal surfaces.There are two length scales defined for ε : a global length scale D, which is on the order of the body size with a global coordinate system x i , and a local length scale d which is proportional to the wavelength of the variation of the unit cell with a local coordinate system y i .The size of the unit cell is assumed to be much smaller compared to the size of the body.The relation between the global coordinate system x i and the local coordinate system y i is defined as where ε is a small scaling parameter.In general, it is equal to the ratio between the size of the unit cell and the size of the macroscopic region to which it belongs.

Formulation of standard mechanics homogenization
Composite materials such as tissue grafts are frequently used to fabricate large structural components.Yet the behavior of these components depends on the composite microstructure.Analyzing large structures at a microstructural level, however, is very difficult due to the structural heterogeneity.Therefore, analysis methods have sought to approximate composite structural mechanics by analyzing a representative section of the composite microstructure, commonly called a representative volume element (RVE).Hill (1963) stated that the RVE was (1) structurally entirely typical of the composite material on average, and (2) contained a sufficient number of irregularities or variety in the structure such that the apparent modulus should be independent of the RVE boundary displacements or tractions.Hashin (1983) also emphasized the nature of RVE analysis stating that stress and strain fields in the RVE should be statistically homogeneous when subjected to homogeneous boundary conditions except in a layer near the external surface.In RVE analysis, representative sections of a structural level are analyzed under assumed boundary conditions to calculate the average or the effective properties of that section.The effective stiffness predictions will change as the number of repeating units contained in RVE changes.This is known as the local analysis.The local analysis provides a matrix relating local strains {ε i } to the global strains {ε i −1 }: where the subscript i denotes the current structural level being analyzed, i -1 denotes the next largest macroscopic level and [M i ] denotes the local structure matrix (Magan and Ripamonti 1996).The average or effective stiffness may then be calculated from the local structure matrix: where [C i-1 ] is the stiffness of the i -1 structural level, [C i ] is the stiffness of the ith structural level, V RVE is the RVE volume and [M i ] is the ith level local structure matrix.
The subsequent macroscopic level can then be analyzed using the previously calculated effective properties.This is known as the global analysis and provides the global strains using Eq.(3).The RVE concept may be applied to estimate strains on multiple structural levels.However, the in situ RVE boundary conditions are not known, and the assumed boundary conditions can only provide an estimate of both effective properties and local strains.In the standard mechanics approach a chosen RVE is generally analyzed using either uniform traction or uniform displacement boundary conditions.These boundary conditions are chosen so as to produce an average strain (if displacements are used) or an average stress (if tractions are applied) within a homogeneous material of the same size as the RVE.The accuracy of these estimates depends significantly on the assumed boundary conditions.Using the minimum potential energy and the minimum complementary energy principles of mechanics, it can be shown that the displacement boundary conditions will provide an upper bound on the effective stiffness.
The relationship between the average strain, the average stress, and the boundary conditions over a RVE can be written using the divergence theorem as where ε i j is the local strain, σ i j is the local stress in the RVE, εi j is the average strain, σi j is the average stress, u i is the displacement, t i is the traction imposed on the RVE boundary, n i is the normal vector to the RVE boundary, y i are the local coordinates of the RVE boundary and S RVE is the RVE boundary.This relation can be further modified in terms of the density ρ.Assuming that the mass of the RVE (for example, the tissue graft) is constant, we can interchangeably use either density or volume.Hence the above two equations can be written as We would like to point out that, in general, the mass of the scaffold will be a time-dependent variable due to the degradation of the scaffold and the new tissue ingrowth.However, in the above equations, we assume it is a constant.
It is important to note that there is no unique relationship between the average stress or strain and the boundary tractions or displacements in the two-or three-dimensional case.In other words, a number of different boundary displacements integrated over the boundary may produce the same average strain.
Once the boundary conditions are chosen, the standard weak form of the equilibrium equations is solved to calculate the local RVE strain.The weak form of the RVE equilibrium equations for the case of applied boundary tractions is solved once for applied stress states corresponding to each component kl of the stress tensor.
where t kl is the boundary traction which would produce an average stress σkl in a homogeneous material, C i jmn is the stiffness tensor of the RVE material components, ε i j (v) is the virtual strain, ε kl mn (u) is the total microstructural strain for the klth traction, and v i is the virtual displacement.
The weak form of the RVE equilibrium equations for the case of applied boundary displacements is where the boundary displacements are implemented using a penalty method and λ is the penalty parameter, u kl i is the boundary displacement, g kl i is the specified displacement which would produce a uniform average strain εkl in a homogeneous material, and the other terms are all defined previously.
Applying three unit displacements in each orthogonal direction in Eq. ( 8) or ( 9) to determine the local stiffness, and then to further determine uniform strain.Once the three strain states are determined, the local structure tensor M i j pm , which relates the average strain εkl pm and the local or microstructural total strain ε kl i j , can be calculated from The local structure tensor M i j pm has minor symmetries such that M i j pm = M ji pm = M i jm p but in general does not have major symmetry M i j pm = M pmi j .Once M i jkl is determined, the local strain at any point within RVE may be calculated from an arbitrary average strain as The effective stiffness tensor Ci jkl , which relates the average strain to the average stress can also be calculated from M i jkl .Starting from Hooke's law at the microscopic level Both sides are integrated over RVE and divided by the total RVE volume to give Substituting for ε kl from Eq. ( 11) and recalling Eqs. ( 4) and ( 5) gives From this the effective stiffness tensor may be defined as

Formulation and finite element implementation of asymptotic homogenization theory
Asymptotic homogenization theory is developed from studies of partial differential equations with rapidly varying coefficients.Two explicit assumptions are made in homogenization theory (Suquet 1987;Hollister et al 1991;Beaupr and Hayes 1985;Williams and Lewis 1982).First, it is assumed that fields vary on multiple spatial scales due to the existence of a microstructure.Second, it is assumed that the microstructure is spatially periodic.To derive the formulation of the asymptotic homogenization theory, let the stress-strain and strain-displacement relations in the heterogeneous body as shown in Figure 1 be expressed as where σ ε i j , e ε i j , E ε i jkl and u i are the stress, strain, stiffness matrix, and displacement within ε , respectively.
Applying the principle of virtual work, one obtains where v i denotes the virtual displacement, and S ε denotes the boundary of all void domains.

Formulation of asymptotic homogenization theory
Consider an asymptotic expansion of u ε i as follows: where the functions u 0 i , u 1 i , u 2 i ,. . .are Y-periodic with respect to the local coordinate y which satisfy the relation . ., with Y being the period of the unit cell.u ε i is the exact value of the field variable, u 0 i is the macroscopic or average value of the field variable, u 1 i , u 2 i , etc. are perturbations in the field variables due to the microstructure.In elasticity theory u 0 i would be the continuum level displacements while u 1 i would be the microstructural displacements.
Derivatives of any function v ε i (x) = v i x, y = x ε with respect to x are written using the chain rule as Substituting Eq. ( 19) into Eq.( 18), and with the help of Eq. ( 20), we obtain For a given ϕ(y) function with Y-periodic, when ε → 0, the relationship of the average of ϕ in the global and in the local domain Y can be defined as where is the global homogenized domain which has the same outside boundary as ε but without the local geometrical details and voids.Y c and S c are the domain and boundary of the solid portion of the unit cell excluding the void and |Y| is the volume of the entire unit cell with void.
Applying Eqs. ( 22) and ( 23) to Eq. ( 21), and by equating the terms with the same power of ε, we obtain Because v i is an arbitrary virtual displacement, it may be chosen to vary only on the macroscopic or microscopic level.In this derivation, we chose v i to vary only on the microscopic level, ie v i = v(y).So on the basis of Eq. ( 24), we found that the first term in Eq. ( 19) is only the function of the global coordinate x, ie We introduce two functions χ kl (χ kl ∈ Y c ) and ϕ (ϕ ∈ Y c ) and let them satisfy the following two equations: Then substituting Eqs. ( 13) and ( 14) into Eq.( 10), we derive where ũ1 i (x) are the arbitrary constants of integration in Y. Inserting Eq. ( 30) into Eq.( 26) and setting Equation ( 30) can be further rewritten as where Here E H i jkl is defined as the homogenized effective stiffness constants calculated from the local unit cell.b i (x) are the averaged body forces and τ i j (x) are the averaged residual stresses caused by the void pressure p i (transport properties) within the unit cell.In the current analysis, we do not consider the effect of the residual stresses, ie we assume that the residual stresses τ i j (x) are zero in the unit cell.
Equation ( 32) is an expression of the principle of virtual work over the homogenized domain with E H i jkl being the effective elastic constants.E H i jkl calculated by Eq. ( 33) is not dependent on the unit cell size due to the periodicity assumption.This differs from the standard mechanics approach in which the effective mechanical properties are usually dependent on the size of the representative volume element due to the St. Venant effect caused by applied displacement or traction boundary conditions.

Finite element implementation
The solution to the homogenized effective elastic constants E H i jkl can be achieved by (a) determining χ kl p from Eq. ( 28) and (b) determining E H i jkl from Eq. ( 33).In general, the solution to Eq. ( 28) is sought through an analogized finite element approach, for example, by expressing E i jkl in a compact matrix form D i j , and denoting where superscripts k, l and i in Eq. ( 37) are defined to take the general compact notations: For example, for a given local coordinate y i and in the case k = l = 1 on the left-hand side of Eq. ( 37) and i = 1 on the right-hand side of Eq. ( 37), the differential operator ε( f ) is defined as Therefore, Eq. ( 28) can then be rewritten as We would like to point out that the form of the expression of Eq. ( 41) is similar to a generalized displacement based finite element formulation if we analogize D i j in Eq. ( 41) as a stiffness matrix, ε(u i ) and ε(v) as strains.Similarly Eq. ( 33) can be rewritten as In doing so, we can directly apply the finite element solution algorithm to solve Eqs. ( 41) and (42) for ε(u i ), and then to determine χ kl p and E H i jkl through the following finite element formulations.For example, after dividing the unit cell into a finite number of elements, we can apply the following shape functions: where v and ûi are the nodal virtual displacements and nodal displacements, respectively; N is the number of nodes in the discrete domain and Ñg m is the global shape function associated with the node m and the element type.
Using the differential operator Eq. ( 40) to Eq. ( 43), we obtain where L is a linear differential operation matrix defined as and is the associated global strain matrix.Now, substituting Eqs. ( 47) and (49) into Eq.( 41), we obtain vT By eliminating v from both sides of Eq. ( 50), since it is a nodal virtual displacement, we obtain The above equation can be rewritten as a standard finite element stiffness equation with K = Y B T D B dY as an analogized element stiffness matrix and f i = Y B T d i dY as analogized nodal forces.In the asymptotic homogenization theory, this nodal 'force' vector is equivalent to the initially applied strain loading.Similarly, Eq. ( 42) can be written as Therefore, we can determine the overall mechanical properties D H i j from Eq. (53).The computational numerical solution algorithm that was developed will be described in the section on computational algorithm.

Homogenization by finite element approach
Due to their complex geometries, most tissue structures are analyzed using numerical techniques, of which the most widely used is the finite element method (FEM).However, standard finite element analysis codes running even on modern, large capacity supercomputers cannot do an analysis at every tissue organizational level.Therefore, FEM also analyzes tissue structural levels by using RVEs.The following section describes a detailed procedure in using FEM to calculate the effective Young's modulus, shear modulus, and Poisson's ratios for the scaffold.

Numerical calculation of Young's modulus
For a given RVE of a scaffold with open pore architecture, a schematic illustration of boundary and loading conditions for calculating the effective constant E XX is shown in Figure 2. As shown in the figure, a uniform displacement field U X = 0.001L X , which is equivalent to a strain ε X = 0.001, is prescribed on the surface S X1 , and a constraint, U X = 0.0, is prescribed on the surface S X2 .
The effective constant E XX can be calculated from the following equation: where A SX1 and A SX2 are the areas of the surface S X1 and S X2 , respectively; L X is the dimension of the RVE in the X direction; and R X is the X-component of the average reaction force produced on the surface S X2 caused by the prescribed displacement U X on the surface S X1 .
It can be determined as Similarly, the above derivation can also be applied to calculate the effective Young's modulus E YY and E ZZ : with U Y = 0.001L X on S Y1 and U Y = 0.0 on S Y2 .
with U Z = 0.001L Z on S Z1 and U Z = 0.0 on S Z2 .
In the above equations A SY2 is the area of the surface S Y2 in the Y direction, A SZ2 is the area of the surface S Z2 in the Z direction, L Y and L Z are the unit cell dimension in Y and in Z, and R Y and R Z are the Y-component and Zcomponent of the average resultant reaction force produced on the S Y2 surface and S Z2 surface, respectively.

Numerical calculation of shear modulus
A schematic illustration of boundary and loading conditions for calculating the effective shear modulus G XY is shown in Figure 3.As shown in the figure, we apply a uniform displacement field U X = 0.001L Y , which is equivalent to γ XY = 0.1%, on the top surface S Y1 , and constrain U X = 0.0 on the bottom surface S Y2 .In order to ensure that a pure shear deformation can be produced in the model, we also apply a constraint U Y = 0.0 on both left and right surfaces S X1 and S X2 , respectively.
The effective constant G XY can be calculated from the following equation: where A SY2 is the area of the surface S Y2 , and R X is the X-component of the average reaction force produced on the surface S Y2 caused by the applied displacement U X on the surface S Y1 .
Similarly, the above derivation can also be applied to calculate the other effective shear moduli, G YZ and G XZ : with U Z = 0.001L Y on S Y1 ; U Z = 0.0 on S Y2 ; and U Y = 0.0 on both S Z1 and S Z2 ; and with U X = 0.001L Z on S Z1 ;U X = 0.0 on S Z2 ; and U Z = 0.0 on both S X1 and S X2 .The reaction forces in Eqs. ( 64) and ( 65) can be determined by

Numerical calculation of Poisson ratios
Calculations for Poisson ratios are based on the following equations: where ε i j , i j = X, Y, and Z, are the results obtained from the FEA models.

COMPUTATIONAL ALGORITHM FOR NUMERICAL SOLUTION
A numerical solution algorithm has to be developed for the finite element implementation of the asymptotic homogenization theory.For a designed porous/cellular RVE, or unit cell, the solution procedure consists of converting a CAD-based unit cell model into a finite element pre-process model, including the model discretization, generating a mesh, and applying boundary conditions and scaffold material properties.Then, the finite element implementation is carried out through our in-house developed computational algorithm: 3DHOMOG.3DHO-MOG is three-dimensional in nature and is capable of analyzing a unit cell with anisotropic material properties.
The functional flowchart of 3DHOMOG is illustrated in Figure 4.As shown in the flowchart, the solution algorithm takes the modeling information, such as the elements, nodes, material, and boundary conditions obtained from the FEM model, computes local stiffness matrices, and then reassembles them into a global stiffness matrix used for the calculation of the effective properties.The global stiffness matrix only needs to be assembled once for all cases (six cases in total for a three-dimensional model).
The local and global force vectors are calculated and assembled for each case.After the boundary conditions are imposed and the global stiffness matrix is inputted into Eq.( 30), the local-level displacements for each of the given cases can be calculated using Gaussian elimination.Subsequently, the overall effective mechanical properties can be determined and be used for further studies such as biomechanical compatibility study between the tissue and the scaffold (Charriere et al 2003).
Store and retrieve stiffness matrix K for all 6 case calculation It can be observed that the effective constants of the scaffolds are, in general, a function of the scaffold materials and the overall porosity of the scaffold structure.The stiffer material can get higher effective material properties compared to the softer material if the microstructure is the same.It can be observed that with the increase in porosity, these three elastic constants decrease.For the special case where the unit cell becomes a fully porous structure, ie the porosity approaches unity, the effective elastic constants approache zero, which agrees with the results reported in Daculsi and Passuti (1990), Charriere et al (2003) and Yan et al (2003).

Effect of pore geometry on the effective properties
To not consider the effect of the scaffolding material, we use normalized effective constants by introducing the following definition: where  As shown in Figures 6 and 7, the mechanical constants generally decrease with the increase in the scaffold porosity for all three scaffold materials.The results in Figure 7 also show that there is no significant difference between the relative effective mechanical properties of the scaffold in terms of the difference in the pore shapes.This indicates that the pore shape may not be a critical parameter when compared to scaffold materials and to the overall scaffold porosity.However, different geometry of the pore shape will result in a different mechanical environment at the microstructure level, as shown in Figure 8 for stress distribution along the inside the pore surface for a square and a circular shape pore.Therefore, pore shape geometry in a scaffold is a very important parameter since it affects the  deformation and stress at the local microstructure level and cells will be influenced by the surrounding environment for ingrowth and migration.

COMPARISON WITH FINITE ELEMENT METHOD
As pointed out earlier, the FEM approach predicts effective properties based on the selected representative unit cell element and on the assumed boundary conditions, and uses the predicted values to represent the overall properties of the entire structure.Since the representative unit cell element is isolated from the entire scaffold the exact boundary conditions of the unit cell are unknown, and uniform strain/stress boundary conditions of the unit cell are usually applied in using FEM to predict the effective properties.Therefore, there are two inherent limitations associated in using this approach: (1) size of the unit cell representative element, and (2) the boundary effect, ie the upper and lower bounds of the predicted effective properties depend on the applied stress or strain boundary condition.In general, to avoid boundary effects, one should choose a representative unit cell element that is as large as possible.Conversely, the difference between the average stress under admissible and actual boundary displacement becomes greater when the representative unit cell element size is bigger.The presented asymptotic expansion homogenization and the computational algorithm which are based on the double scale techniques to establish stress and strain relations at both micro-and macrostructural level, provide an alternative approach to predict the effective mechanical properties of heterogeneous porous structures with periodic or quasi-periodic microstructure.In addition to its advantage of being able to analyze at both local and global levels, the asymptotic homogenization theory does not have the bound error and the boundary effect when predicting the effective constant unlike conventional analytical and FEM approaches.The effect of the number of elements when using homogenization theory and using the FEM approach on the prediction of the effective modulus for a designed biocompatible hydroxyapatite tissue scaffold is presented in Figure 9.As shown in the figure, the accuracy of the FEM result depended on the number of elements used; for example, in the case shown in the figure, the FEM results tended to be stable when the number of elements was greater than 468, while the results obtained from the asymptotic homogenization model was much less sensitive to the number of elements used in its finite element implementation.
Figure 10 shows the results of the comparison of the 'size' effect for the same scaffold shown in Figure 9.There were 1, 8, and 27 cells used in the standard FEM model, and only 1 cell was used in the homogenization analysis because of the periodicity assumption.A 'cell' here refers to the base unit cell repeating in the scaffold.To minimize the effect of the number of the elements, we retained the same number of elements and mesh for each base cell (about 700 eight-node brick elements).The effective constant E XX was calculated using Eq. ( 54).Results of the prediction clearly showed that the FEA results depended on the number of base cells used, and the FEA prediction tended to be stable only when enough base cells were used.A comparison of the results predicted by using the FEM approach and by using the homogenization approach is presented in Figure 10.

CONCLUSIONS AND DISCUSSIONS
A comparison of three available homogenization approaches, ie standard mechanics modeling, homogenization theory, and finite element method, for their application in the prediction of the effective properties of heterogeneous and porous tissue scaffold was presented.A computational algorithm that was developed in-house for the finite element implementation of the asymptotic homogenization theory was introduced.Details of the homogenization formulations, solution procedure, and discussions on the advantages and limitations of the approaches were described.Application of the computational algorithm in the characterization of the effective mechanical properties of tissue scaffolds with variable design parameters and the results of the parametric analyses were presented.
Results of the predictions show that the effective constants of the scaffolds are, in general, a function of the scaffold materials and the overall porosity of the scaffold structure.The effective constants decrease with the increase of porosity for all three scaffold biomaterials, as shown in Figure 6.The results of the analyses also show that the pore shape may not be a critical parameter when considering the overall scaffold structural properties.However, different pore shapes will result in different mechanical stress distributions at the microstructure level, where it may have a critical effect upon the cells that will interact with the scaffold at microscopic levels.
As indicated earlier, both standard mechanics modeling and finite element methods have inherent limitations in their application-for example, the upper and lower bound differences resulting from the size and boundary effect of the representative volume element (or unit cell) used in the calculation of the unit cell modules.The computational algorithm that was described and the asymptotic expansion theory-based homogenization may reduce these effects.This has been shown by comparing the algorithm's predictions and the finite element model's predictions.It can be observed from Figures 9 and 10 that the accuracy of the FEA result relied on the number of elements.The results of the developed algorithm and the asymptotic homogenization theory were much less sensitive to the number of the elements used, as compared to the FEA results which were more sensitive.In addition, since the asymptotic homogenization theory applies asymptotic expansion and double scale techniques to establish stress and strain relations at both micro-and macrostructural level in the calculation of effective mechanical properties, the developed computational algorithm may provide an important tool to allow the future study of biological and biomechanical properties at both the global (scaffoldtissue) level and local (cell-scaffold) level.

Figure 3
Figure 3 Loading and boundary condition in determining G XY .

Figure 6
Figure 6 Effect of the scaffolding materials and porosity (square shape pore) on (a) Young's modulus, (b) shear modulus and (c) Poisson's ratio.

Figure 8
Figure 8 Comparison of local stress distribution along inside (a) square pore surface and (b) circular pore surface [29].

Figure 9
Figure 9 Effect of the number of elements on model predictions: (a) unit cell (scaffold) model and (b) model predictions.

Figure 10
Figure 10 Effect of the number of base cells on the predictions.
E r , G r , υ r are the dimensionless relative effective Young's modulus, shear modulus, and Poisson's ratio, respectively.E * , G * , υ * are the predicted effective constants, and E s , G s , υ s are the scaffold material constants for Young's modulus, shear modulus, and Poisson's ratio, respectively.Comparisons of the predicted effective mechanical properties of the scaffolds with square pores and