## Introduction

Human tibia is located in the lower front portion of the leg of a person and carries a significant part of the body weight. It consists of soft spongy ends called the trabecular regions with compact intermediate bone referred to as the cortical region. The bone is subjected to a variety of loads and performs all its functions due to its mechanical properties that vary gradually from point to point everywhere inside the bone. Due to the deposition of various elements at different rates in the bone matrix, the bone material becomes highly inhomogeneous and anisotropic. Considering the huge variation in the material distribution inside the bone, exact point-to-point determination of mechanical properties and assigning those properties in the analysis of the bone is a cumbersome task. However, there have been some experimental studies wherein the mechanical properties of tibia have been determined in different locations in the tibia [112]. Some computational studies [1320] have incorporated this variation in the mechanical properties in determining the stress distribution in the bone. However, these studies involve jumps in mechanical properties at interfaces which could lead to irregular stress behavior around these interfaces (as will be shown in this paper) affecting the accuracy of the stress prediction. In this study, we have developed smooth functionally graded mechanical properties using the discrete data available in the literature and used them in a finite element analysis (FEA) of the human tibia under a static load.

All the previous computational studies related to human tibia involve development of a 3D geometric model of the bone, experimental determination of the material properties, and performing FEA of the model with the help of experimentally obtained properties. To develop a 3D geometric model of bone, computed tomography (CT) is used extensively as a diagnostic tool in biomechanical research. There are several studies that have used CT scans for bone modeling [2128]. In our work, we have followed the approach of the previous studies and have used CT images to develop a 3D geometric model of the human tibia.

As far as the material properties are concerned, one of the earliest studies for the determination of material properties of bone was performed by Dempster and Liddicoat [1] in which they experimentally established that compact bone is a nonisotropic or orthotropic material. Vincentelli and Evans [2] studied the relation between the mechanical properties (like ultimate tensile stress and strain, modulus of elasticity) and the collagen fiber orientation and degree of calcification in the cortical region of human tibia. Currey [3] and Reiliy and Burstein [4] presented a detailed review of studies on the mechanical properties of bone including tibia. Knets et al. [5] measured the material parameters of the human tibial cortical bone through mechanical experimentation by the application of tensile loading at different locations of the bone which established the cross-sectional inhomogeneity in the cortical region in addition to the axial inhomogeneity because of the presence of the trabecular and cortical regions. Pfafrod et al. [6] reported another remarkable work on cross-sectional inhomogeneity by describing the relationship between the biochemical composition and torsional strength in three orthogonal directions in six zones of the cortical region. Lindahl [7] performed compression test on spongy bone and calculated various parameters including the limit of proportionality and modulus of elasticity. Williams and Lewis [8] determined the stiffness of cancellous region in the proximal human tibia in three orthogonal directions and suggested that cancellous bone has transversely isotropic nature. Further, Goldstein [9] presented a review of studies on the physical properties of trabecular bone. In a study carried out by Ciarelle et al. [10], the orthogonal mechanical properties (elastic moduli and density) of human trabecular bone were evaluated with the help of materials testing and CT images. Ashman et al. [11] and Rho et al. [12] determined mechanical properties of cancellous bone using an ultrasound transmission technique and established correlations for predicting elastic modulus from density and CT numbers in human bone. In our study, we will incorporate all these features viz. axial as well as cross-sectional inhomogeneity and orthotropy in the mechanical properties.

In the studies related to FEA, Hayes et al. [13] developed an axisymmetric finite element (FE) model of the lateral tibial plateau to predict stresses and displacements in the proximal tibia (upper trabecular region). The model predicted maximum compressive stresses beneath the center of the contact region and maximum shear stresses beneath the edge of the contact region for trabecular bone. Little et al. [14] proposed a three-dimensional FE model of the proximal tibia and analyzed the bone by assuming it to be linearly elastic and isotropic for the evaluation of strength. In another study, Mehta and Rajani [15] developed one material model, three material models, and four material models of tibia, respectively. In these different models, the effect of varying material properties on bone stresses in human tibia and fibula were observed. Muller-Karger et al. [16] also developed similar models by considering cases such as whole tibia (with and without the intramedullary canal) and only the proximal tibia. They further studied these models assuming bone to be transversely isotropic. Hull coworkers [17] studied the contact behavior of the knee joint by developing an FE model in which the cancellous region of tibia was considered isotropic, whereas the cortical part was assigned orthotropic properties. In recent studies by Ionescu et al. [18] and Verim coworkers [19], the tibia was analyzed by dividing the bone into three parts (proximal part, shaft, and distal part) and each part was assigned a different set of elastic properties, thereby incorporating the axial inhomogeneity. However, the cross-sectional inhomogeneity existing in the cortical region was not considered in these FEA studies. Kerem and Calik [20] further included the cross-sectional inhomogeneity as reported by Knets et al. [5] and Pfafrod et al. [6]. However, all these previous studies related to FEA of human tibia were performed by assigning material properties to trabecular and cortical regions of the bone in a discrete or nonfunctionally graded manner. In other words, a gradual change in the material properties of the bone was not considered previously. Although these types of analyses give a very good insight, they also lead to some unrealistically spiked stress values at the interfaces of different material regions inside the bone. However, if the material properties are assigned to the bone in a continuous manner, the stresses can be predicted more appropriately. This type of assignment of material properties can be done if one functionally grades the material properties of the bone. Against this backdrop, we have carried out FEA of tibia by considering it a functionally graded structure.

In this work, the material properties of human tibia were obtained from the literature and functional grading was done by smoothly varying the material properties in the form of analytic functions. These smooth functions for the material properties were subsequently used for analyzing the stress distribution in the tibia model under static loading, which is akin to the one encountered during a normal gait. We first compare the stress distribution of functionally graded tibia with the nonfunctionally graded one for isotropic and cross-sectionally homogeneous model to highlight the importance of functional grading in correctly ascertaining the stresses. We then study the effect of the orthotropy and the variation of mechanical properties along the cross section on the developed stresses.

The rest of the paper is organized as follows. In Sec. 2, we present the methodology used for the analysis including the geometric modeling, development of functionally graded material properties, and the development of the finite element model appropriate for the further analysis. Detailed results are presented in Sec. 3 along with discussions. We finally conclude in Sec. 4.

## Methodology

The methodology involved developing an approximate 3D model of human tibia and analyzing the model using finite element method (FEM) in comsolmultiphysics. In the analysis, following four types of models have been considered for analysis: (i) a nonfunctionally graded isotropic model (NFIM), (ii) a longitudinal functionally graded isotropic model (LFIM), (iii) a longitudinal functionally graded orthotropic model (LFOM), and (iv) a three-dimensional functionally graded orthotropic model (TFOM) of tibia. Inhomogeneity along the length of the bone has been included in all the models. The steps involved to perform the analysis are discussed in Sec. 2.1.

### Geometric Model Formation.

The geometric model for the study was formed by performing CT scan on the legs of a person (a healthy 24 years old male student of IIT Kanpur having BMI 25.3, height 154 cm, and weight 60 kg at the time of CT scan with no bone deformity). CT imaging (image resolution 1280 × 1024, pixel size 0.29 mm) was performed on Siemens made CT scanner (Somatom Emotion, 16-slice model, Germany) at Ratan MRI center, Kanpur, Uttar Pradesh. Written consent from the person for CT scan was taken in advance. 41 slices covering the total length of the tibia (approximately 320 mm) containing its transverse cross section were taken at an interval of 8 mm. A typical CT slice of the human leg is shown in Fig. 1(a), where two distinct bony structures, viz., the tibia and the fibula can be identified. A zoomed portion of Fig. 1(a) around the tibia is shown in Fig. 1(b) and the solid model obtained by stacking such zoomed views of each cross section is shown in Fig. 1(c). The 3D model was developed by importing all the 41 CT images in Autodesk Inventor, extracting the contours of the cross section and lofting them. The top surface of the 3D model thus obtained has been split into two parts by a centroidal line so that one part (the right one in Fig. 1(c)) represents the medial half and the other (the left half in Fig. 1(c)) represents the lateral half. This differentiation is necessary in our analysis for the application of the load to be consistent with the one recorded in a normal gait cycle as will be discussed later in Sec. 2.3. We will use the Cartesian coordinate system XYZ with the XY-plane representing the cross-sectional plane and the Z-axis parallel to the bone axis as shown in Fig. 1(c). The origin of the coordinate system is such that $z=0$ represents the top surface, while $z=320$ represents the bottom surface.

### Material Properties.

In the analysis, material properties namely Young's modulus (E), Poisson's ratio (ν), and shear modulus (G) were taken from literature [5,6,11,1820]. As discussed earlier, the material properties in different studies were assigned to the model in four different manners. Since the medullary cavity does not contribute significantly toward the bone strength [18], no material properties were assigned to it in all the studies. As mentioned earlier, the bone was assumed to be isotropic and the analysis was carried out without grading the material properties functionally in the first study (i.e., NFIM). The material properties assigned to the trabecular and cortical regions (as shown in Fig. 1(c)) of NFIM are given in Table 1 [20]. The second study was carried out by functionally grading the isotropic bone structure that is discussed later in this section. Contrary to the first two studies, the bone was considered to be orthotropic in studies 3 and 4. While the third study assumed the bone to be orthotropic and cross-sectionally homogeneous, the fourth study assumed the bone to be orthotropic and three-dimensionally inhomogeneous. In all the studies, the locations of interfaces were taken at $z=64 mm$ in the proximal tibia and $z=256 mm$ in the distal tibia. The reason for this choice is that although there is no exact demarcation between the trabecular and cortical regions in the bone, some researchers have considered 20% of the total length of the bone to be trabecular in both the proximal and distal tibia [25,29,30]. It is important to mention here that the studies 2, 3, and 4 were carried out after functionally grading the material properties.

We note that the mechanical properties of the trabecular and cortical regions of the bone are typically measured by isolating a small specimen from an interior region and assuming that these are distinct regions with sharp boundaries within which the properties remain constant. Since the mechanical properties depend on the local level of calcification [2], which changes gradually in space, these material properties should also vary gradually as we move from one region to the other leading to a diffused interface instead of a sharp interface. To capture this feature, the generalized functional form of the material properties (E, ν, and G) used in these studies (i.e., study 2, 3, and 4) is chosen to be
$Ptibiax,y,z=Ptrabecular+0.5 * Pcorticalx,y,z−Ptrabecular * 1+tanh0.1 * z−64+0.5 * Ptrabecular−Pcortical x,y,z * 1+tanh0.1 * z−256$
(1)

where P represents the material property to be graded (i.e., E, ν, or G). With this functional form, there is a gradual change in the properties from the trabecular to the cortical region along the axial direction without significantly affecting the properties in the interior of these regions. We note that there is no data related to cross-sectional inhomogeneity in the trabecular region, and accordingly, $Ptrabecular$ is a constant depending on the material parameter under consideration. The various assigned values for these constants are given in Tables 1 [20] and 2 [11,1820] according to the type of study being carried out. For example, in study 1 and 2, material property Young's modulus of trabecular region (i.e., $Etrabecular$) has been assigned a value of 9.64 GPa, whereas in study 3 and 4, the values assigned to $Etrabecular$ are 4.48 GPa, 4.48 GPa, and 9.64 GPa along the X-axis, Y-axis, and Z-axis, respectively. In the second and third study, we will consider a cross-sectionally homogeneous model for the cortical region as well. Hence, we will assume $Pcorticalx,y,z=K1$, where $K1$ is a constant depending on the material property under consideration (again reported in Tables 1 and 2). In the fourth study, cross-sectional inhomogeneity in the cortical region has also been incorporated. For this, we first grade the material properties (E, ν, and G) by forming functions for their angular/circumferential variation in the cortical region as per the best functional fit to the discrete data reported in literature (see Table 3) [5,6,20]. The procedure followed to obtain this functional fit is described next.

#### Procedure for Modeling Cross-Sectional Inhomogeneity in the Cortical Region.

Since cross-sectional inhomogeneity has been reported only in the cortical part of the tibia, we select a typical cross section in the cortical part of the bone as shown in Fig. 2 (the cross section has centroid at point C (xc, yc)). This cross section is further marked with six material regions according to the literature [5,6,20] such that each material region covers an angular range of 60 deg (i.e., region 3 from 0 deg to 60 deg, region 2 from 60 deg to 120 deg, region 1 from 120 deg to 180 deg, region 6 from 180 deg to 240 deg, region 5 from 240 deg to 300 deg, and region 4 from 300 deg to 360 deg (or 0 deg) as shown in Fig. 2.

The material properties of these six regions taken from Refs. [5,6], and [20] are reported in Table 3. We note that the coordinate system reported in Refs. [5,6], and [20] and used in this work is different which necessitated certain transformations. As a result, our Table 3 is slightly different from a similar table in Ref. [20]. In order to grade these properties circumferentially in the form of analytical functions, we have assumed that they represent the properties at the centerline of each region, i.e., lines at angular locations of 30 deg, 90 deg, 150 deg, 210 deg, 270 deg, and 330 deg from the X-axis in the material regions 3, 2, 1, 6, 5, and 4, respectively, as shown by dotted lines in Fig. 2. Six representative material points namely P1, P2, P3, P4, P5, and P6 on these angular positions are shown in Fig. 2. Using the angular positions (θ) and corresponding values of material properties (E, ν, and G) as given in Table 2, the best fit material property functions in the form of
$Pcorticalθ=a+b cos w1θ+c sin w2θ$
(2)
is obtained using the curve fitting toolbar of matlab. In the above, $Pcorticalθ$ is the property (E, ν, or G) to be graded and a, b, w1, and w2 are the fitted parameters. Depending on the property under consideration, we find that w1 and w2 come out to be 1, 2, or 3 as per the best fit. Accordingly, the terms cos w1θ and sin w2θ in Eq. (2) become cos θ, sin θ, cos 2θ, sin 2θ, cos 3θ, and sin 3θ.
For a particular case when $w1=w2=1$, Eq. (2) can be written in terms of x, y, z by replacing sine and cosine terms in Cartesian coordinates as
$Pcorticalx, y= a+b * x−xcx−xc2+y−yc2+c * y−ycx−xc2+y−yc2$
(3)
In Eq. (3), xc and yc represent the cross-sectional centroids at a given axial location (z). The functional dependence of xc and yc on z is obtained by locating the x and y coordinates of the centroids of all the cross section of cortical region of the bone and observing their variation along the length of the bone. This variation fits very well the linear forms
$xc= 0.0692 * z+159$
(4)

$and yc=−0.05334 * z−110.9$
(5)

The linear fit is based on studies reported earlier for intramedullary canal axis of human tibia [31,32].

Substituting $xc$ and $yc$ in Eq. (3), circumferential variation of material properties (i.e., inhomogeneity) in the cortical region of tibia can be expressed as
$Pcorticalx, y, z=a+b * x−0.0692 * z+159x−0.0692 * z+1592+y−−0.05334 * z−110.92+c * y−−0.05334 * z−110.9x−0.0692 * z+1592+y−−0.05334 * z−110.92$
(6)
Using Eq. (6), Eq. (1) for the property to be graded can be written as
$Ptibiax, y, z=Ptrabecular+0.5 * (a+b * x−0.0692 * z+159x−0.0692 * z+1592+y−−0.05334 * z−110.92+c * y−−0.05334 * z−110.9x−0.0692 * z+1592+y−−0.05334 * z−110.92) * tanh0.1 * z−64−tanh0.1 * z−256+0.5 * Ptrabecular * (tanh0.1 * z−256−tanh0.1 * z−64)$
(7)
Such a function expressed in Eq. (7) takes into account the overall spatial variation of the material property (cross-sectional variation in the cortical part and the axial variation throughout the length) of human tibia. The functional forms of all the relevant material properties (including those where $w1, w2=2 and 3$) formed in this manner are used in this work. These functional descriptions of the material properties were applied to the human tibia model in comsolmultiphysics for the FEA. As an example, the best fit functional form for the Young's modulus used in TFOM along Y-axis (Ey) in the cortical part of tibia is
$Ey corticalθ=6.786E9+0.693E9 cos 3θ−1.456E9 sin θ$
(8)
which can be further written as
$Ey corticalθ= 6.786E9+0.693E9 4 cos3θ−3 sin θ−1.456E9 sin θ$
(9)
Replacing sine and cosine terms in the above equation by Cartesian coordinates, we get
$Ey corticalx, y= 6.786E9−3.535E9y−ycx−xc2+y−yc2+2.772E9x−xcx−xc2+y−yc23$
(10)
Using the linear fit for the variation of xc and yc along the Z-axis from Eqs. (4) and (5), we get the following form of Ey in the cortical region:
$Ey corticalx, y, z= 6.786E9−3.535E9y−−0.05334 * z−110.9x−0.0692 * z+1592+y−−0.05334 * z−110.92+2.772E9x−0.0692 * z+159x−0.0692 * z+1592+y−−0.05334 * z−110.923$
(11)
Substituting $Ey cortical$ (x, y, z) in Eq. (1) and after simplifying we get
$Eyx,y, z=4.48E9+(3.393E9−1.767E9y−−0.05334 * z−110.9x−0.0692 * z+1592+y−−0.05334 * z−110.92+1.386E9x−0.0692 * z+159x−0.0692 * z+1592+y−−0.05334 * z−110.923−2.24E9)tanh0.1z−64−tanh0.1z−256$
(12)

This functional form of Young's modulus takes into account the overall spatial variation of Ey (cross-sectional variation in the cortical region and the axial variation throughout the length) of human tibia. Functional forms for the variation of the other material properties (i.e., Ex, Ez, νxy, νyz, νxz, Gxy, Gyz, and Gxz) were also formed in a similar manner and are reported in the Appendix.

### Boundary Conditions.

The computational study was done using a 64-bit Intel(R) Xeon(R) processor (3.30 GHz) with 16 GB RAM. We have adopted the boundary and loading conditions from the literature [20,33] and accordingly we have kept the distal end ($z=320 mm$) fixed and prescribe the load at the proximal end ($z=0 mm$) in all the studies. A total load of 2450 N has been applied which is a typical value encountered during normal gait at the stance phase in near full extension (taken from Harrington [33]). Since the knee contact force acts in such a manner that 60% of the load is borne by the medial condyle whereas, 40% by the lateral condyle [20,34,35], we have split the total load on the top surface such that a uniformly distributed normal compressive traction totalling 1470 N acts on the medial half shown in Fig. 1(c) and another uniformly distributed normal compressive traction with a net force of 980 N acts on the lateral half as shown in Fig. 1(c).

### Discretization and Convergence Study.

After assigning the material properties, FEA is performed on the bone model in comsolmultiphysics after suitable discretization. We have used the ten-noded tetrahedral element for the discretization in all the four studies. In order to ascertain the convergence of the FEM results, we have used the equivalent von Mises stress values (typically used in the failure analysis) as the measure. In this work, we will be reporting the maximum von Mises stress in the entire tibia as well as the distribution of the von Mises stress in several distinct transverse planes and the longitudinal YZ plane at $x=175 mm$ situated inside the bone as shown in Figs. 3(a) and 3(b). Convergence of the results implies that both the maximum value as well as the stress distribution is captured properly by the given discretization. However, checking convergence of stress distribution over an entire plane is rather cumbersome, and therefore, we will also report the stress values along a line joining the points (175, −112, 0) and (175, −120, 320) lying on the YZ plane as shown in Fig. 3(c). Separate convergence tests were performed on the NFIM and the TFOM to get the optimum mesh size in each case. The convergence plots are shown in Figs. 4 and 5, wherein the von Mises stress along the line shown in Fig. 3(c) is plotted for the NFIM and LFIM, respectively. It can be noted from Figs. 4 and 5 that for both the cases, the stress values for the mesh with 508,006 and 771,033 nodes are almost the same. Therefore, we will employ a mesh consisting of 508,006 nodes for the further analysis.

## Results and Discussion

Prediction of accurate von Mises stress behavior in the bone is quite significant for assessing failure and the associated fracture risk [36]. In this work, the distribution of von Mises stresses in the nonfunctionally graded model and different functionally graded models of human tibia were analyzed and compared with each other. An important feature of this work is that we have analyzed the results both qualitatively as well as quantitatively. As mentioned earlier, we have investigated the stress variation along the cut line, the longitudinal plane ($x=175 mm$) and at several transverse planes including the distal interface ($z=256 mm$) of tibia in all the bone models. For the qualitative analysis, we have observed the pattern and distribution of stresses in the line plots and the contour plots of von Mises stress obtained in the four studies as shown in Figs. 69. For quantitative analysis, we have determined the maximum von Mises stress for the entire tibia as well as along the cut line, on the longitudinal plane and the distal interface of tibia, in the four studies and compared them. These stress values are given in Table 4.

Figure 6 clearly shows that the NFIM of tibia has a jump in the von Mises stress at the two interfaces. Since there is a discontinuity in the material properties at these interfaces, there could have been jumps in the strain field but the stress field should have been continuous (from static force equilibrium). Hence, these jumps are erroneous and are artifacts of the approximations involved in the FEA. Also, there is a sharp gradient in the stresses on either side of the jumps. Depending on the magnitude of the jump and the location of the interface, one could wrongly conclude the critical locations to be lying on the interface and also the magnitude of the maximum von Mises stress could be wrongly predicted. In addition, there is an uneven distribution of stresses at the distal interface of the nonfunctionally graded model as is clear from the contour plots shown in Fig. 9(a). These features are undesirable for the correct biomechanical analysis. Contrary to this, in the LFIM, the line plot of the stresses (see Fig. 6) is continuous and free from any such artifact. Also, there is an even distribution of stresses at the distal interface as is clear from the contour plot shown in Fig. 9(b). This smoothening happens because functional grading leads to a smooth and continuous variation in the material properties describing the model which is reflected as a smooth variation in the results from the model. Therefore, functionally graded models are most suitable for the correct analysis of bone. For the LFIM, a reduction in the maximum von Mises stress in the entire tibia by 31.24%, in the cut line by 0.19%, in the longitudinal plane by 2.51% and at the distal interface by 49.82% was observed in comparison to the NFIM as shown in Table 4.

The comparative study presented in Fig. 6 highlights the significance of functional grading of material properties in bone analysis but is not a true reflection of the actual bone behavior. We have extended the concept of functional grading to develop a bone model which mimics as closely the actual bone as possible. This is done by considering LFOM and TFOM, respectively. In LFOM considered in study 3, functional grading has been done after considering orthotropic nature of the tibia while ignoring the cross-sectional variations in the material properties. The variation of von Mises stresses along the cut line in this study too is continuous all along the length as shown in Fig. 7. However, the maximum stress values in the entire tibia and in the longitudinal plane increased by 4.4% and 6.19%, respectively, in comparison to LFIM. Whereas, in the cut line and at the distal interface stress values reduced by 1.87% and 4.71%, respectively, as compared to LFIM (see Table 4). Therefore, analysis of the LFOM suggests that inclusion of orthotropy increases the maximum stress values in the entire tibia.

In the studies 1–3, we have assumed the bone to be cross-sectionally homogeneous. However, as mentioned before, the bone is inhomogeneous in nature with significant change in material properties in the circumferential direction. Therefore, in TFOM, circumferential inhomogeneity of the cortical region has also been incorporated along with homogeneous trabecular regions.2 It can be seen from the line plots in Fig. 8 that incorporation of circumferential inhomogeneity decreases the von Mises stress in the cortical region along this line.

However, it is evident from Table 4 that the maximum stress values in the entire tibia is almost similar to LFOM but a slight increase of 0.44%, 0.65%, and 6.27% in the cut line, longitudinal plane and the distal interface, respectively, was observed in comparison to LFOM. The contour plots showing variation of the von Mises stress across the distal interface for TFOM in Fig. 9(d) hardly shows any variation from LFOM (Fig. 9(c)). Hence, to better understand the effect of circumferential inhomogeneity on the stress distribution, we have further investigated the variation in stress values in different transverse planes (shown in Fig. 3(a)) situated inside the bone. The maximum von Mises stress values obtained at these different transverse planes for the LFOM and TFOM models are given in Table 5. It can be concluded from the table that the primary effect of the inhomogeneity is to increase the stress values on the transverse planes. To better understand the stress distribution pattern, we further compare the stress contours at other transverse planes, in particular, the plane $z=192$ mm, for which a comparison between the LFOM and TFOM shown in Figs. 10(a) and 10(b) indicates that although the maximum von Mises stress increases marginally in TFOM, there is a smaller region of high stresses as compared to LFOM. Hence, circumferential inhomogeneity effectively leads to a more even distribution of the stresses across the transverse cross section.

## Conclusions

In this work, we have studied the influence of functional grading of material properties on the von Mises stress distribution in a typical human tibia. Moreover, apart from predicting maximum stress, the location and the manner in which they are distributed inside the bone have been analyzed in this work. Our results indicate that nonfunctionally graded models are not appropriate for the analysis due to presence of artificial interfaces. Moreover, it was observed that functionally graded bone model can predict different von Mises stress values in the bone as compared to the nonfunctionally graded model. Therefore, functional grading of the bone has to be taken into account to obtain a reasonable prediction of stress distribution in the bone, which could be important in the studies related to analysis of bone strength. Furthermore, our study also highlights the importance of orthotropy and cross-sectional inhomogeneity of the material properties in the correct distribution of the stresses throughout the human tibia. We find that orthotropic nature of the bone tends to increase the maximum von Mises stress in the entire tibia. Further, the inclusion of cross-sectional inhomogeneity typically increases the stresses across normal cross section. However, it makes the stress distribution more even. Hence orthotropy and cross-sectional inhomogeneity help in better identification of the critical zones in the tibia from the fracture point of view under a given loading. Accordingly, we conclude that both orthotropy as well as cross-sectional inhomogeneity should be included to correctly capture the stress distribution.

## Acknowledgment

The first author would like to thank Kapil Manoharan, Shashank Sharma and Ajay Bhandari for their meaningful suggestions in the work. This work was supported through the grant given to Niraj Sinha by Science and Engineering Research Board, India (Grant No. SR/FTP/ETA-464/2012).

## Funding Data

• Science and Engineering Research Board, India (Grant No. SR/FTP/ETA-464/2012; Funder ID: 10.13039/501100001843).

## Nomenclature

• E =

Young's modulus of the bone

•
• Ex =

Young's modulus of the bone along the X-axis

•
• Ey =

Young's modulus of the bone along the Y-axis

•
• Ez =

Young's modulus of the bone along the Z-axis

•
• G =

shear modulus of the bone

•
• Gxy =

shear modulus of the bone in XY plane

•
• Gxz =

shear modulus of the bone in XZ plane

•
• Gyz =

shear modulus of the bone in YZ plane

•
• ν =

Poisson's ratio of the bone

•
• νxy =

Poisson's ratio of the bone in XY plane

•
• νxz =

Poisson's ratio of the bone in XZ plane

•
• νyz =

Poisson's ratio of the bone in YZ plane

2

Inhomogeneity in the trabecular region has not been included due to unavailability of data for the same.

### Appendix: Global Functional Forms Used forMaterial Properties in Different Studies of Human Tibia

##### Study 2: Longitudinal Functionally Graded Isotropic Model.

In study 2, the functionally graded Young's modulus, Poisson's ratio, and shear modulus of tibia were assigned in the following functional forms:

Young's modulus (Pa)
$Ezz=9.64E9+4.38E9tanh0.1z−64−tanh0.1z−256$
(A1)
Poisson's ratio
$νz=0.35−0.11tanh0.1z−64−tanh0.1z−256$
(A2)
##### Study 3: Longitudinal Functionally Graded Orthotropic Model.

In study 3, the functionally graded Young's modulus, Poisson's ratio, and shear modulus of tibia were assigned in the following functional forms:

Young's modulus (Pa)
$Exz=4.48E9+2.015E9tanh0.1z−64−tanh0.1z−256$
(A3)

$Eyz=4.48E9+1.215E9tanh0.1z−64−tanh0.1z−256$
(A4)

$Ezz=9.64E9+4.38E9tanh0.1z−64−tanh0.1z−256$
(A5)
Poisson's ratio
$νxyz=0.35+0.07tanh0.1z−64−tanh0.1z−256$
(A6)

$νyzz=0.12$
(A7)

$νxzz=0.12+0.01tanh0.1z−64−tanh0.1z−256$
(A8)
Shear modulus (N/m2)
$Gxyz=1.41E9+E9tanh0.1z−64−tanh0.1z−256$
(A9)

$Gyzz=1.28E9+1.14E9tanh0.1z−64−tanh0.1z−256$
(A10)

$Gxzz=1.28E9+1.815E9tanh0.1z−64−tanh0.1z−256$
(A11)
##### Study 4: Three-Dimensional Functionally Graded Orthotropic Model.

In study 4, the functionally graded Young's modulus, Poisson's ratio, and shear modulus of tibia were assigned in the following functional forms:

Young's modulus (Pa)
$EXx,y, z=4.48E9+(4.219E9−1.668E9y−−0.05334 * z−110.9x−0.0692 * z+1592+y−−0.05334 * z−110.92+1.53E9x−0.0692 * z+159x−0.0692 * z+1592+y−−0.05334 * z−110.923−2.24E9)tanh0.1z−64−tanh0.1z−256$
(A12)

$Eyx,y, z=4.48E9+(3.393E9−1.767E9y−−0.05334 * z−110.9x−0.0692 * z+1592+y−−0.05334 * z−110.92+1.386E9x−0.0692 * z+159x−0.0692 * z+1592+y−−0.05334 * z−110.923−2.24E9)tanh0.1z−64−tanh0.1z−256$
(A13)

$EZx,y, z=9.64E9+(8.024E9+x−0.0692 * z+159x−0.0692 * z+1592+y−−0.05334 * z−110.92(0.562E9x−0.0692 * z+159x−0.0692 * z+1592+y−−0.05334 * z−110.92−3.321E9y−−0.05334 * z−110.9x−0.0692 * z+1592+y−−0.05334 * z−110.92)−4.82E9)tanh0.1z−64−tanh0.1z−256$
(A14)
Poisson's ratio
$νxyx,y, z=0.35+(0.23+0.0098x−0.0692 * z+159x−0.0692 * z+1592+y−−0.05334 * z−110.92−0.058y−−0.05334 * z−110.9x−0.0692 * z+1592+y−−0.05334 * z−110.92x−0.0692 * z+159x−0.0692 * z+1592+y−−0.05334 * z−110.92−0.175)tanh0.1z−64−tanh0.1z−256$
(A15)

$νxzx,y, z=0.12+(0.08+0.006x−0.0692 * z+159x−0.0692 * z+1592+y−−0.05334 * z−110.92+0.037y−−0.05334 * z−110.9x−0.0692 * z+1592+y−−0.05334 * z−110.92x−0.0692 * z+159x−0.0692 * z+1592+y−−0.05334 * z−110.92−0.06)tanh0.1z−64−tanh0.1z−256$
(A16)

$νyzx,y, z=0.12+(0.068+0.011x−0.0692 * z+159x−0.0692 * z+1592+y−−0.05334 * z−110.92+0.031y−−0.05334 * z−110.9x−0.0692 * z+1592+y−−0.05334 * z−110.92x−0.0692 * z+159x−0.0692 * z+1592+y−−0.05334 * z−110.92−0.06)tanh0.1z−64−tanh0.1z−256$
(A17)
Shear modulus (N/m2)
$Gxyx,y, z=1.41E9+(1.247E9−0.205E9y−−0.05334 * z−110.9x−0.0692 * z+1592+y−−0.05334 * z−110.92+0.415E9x−0.0692 * z+159x−0.0692 * z+1592+y−−0.05334 * z−110.923−0.705E9)tanh0.1z−64−tanh0.1z−256$
(A18)

$Gxzx,y, z=1.28E9+(2.453E9−0.853E9y−−0.05334 * z−110.9x−0.0692 * z+1592+y−−0.05334 * z−110.92+0.87E9x−0.0692 * z+159x−0.0692 * z+1592+y−−0.05334 * z−110.923−0.64E9)tanh0.1z−64−tanh0.1z−256$
(A19)

$Gyzx,y, z=1.28E9+(1.758E9−1.662E9y−−0.05334 * z−110.9x−0.0692 * z+1592+y−−0.05334 * z−110.92+1.646E9x−0.0692 * z+159x−0.0692 * z+1592+y−−0.05334 * z−110.923−0.64E9)tanh0.1z−64−tanh0.1z−256$
(A20)

## References

References
1.
Dempster
,
W. T.
, and
Liddicoat
,
R. T.
,
1952
, “
Compact Bone as a Non-Isotropic Material
,”
Am. J. Anat.
,
91
(
3
), pp.
331
362
.
2.
Vincentelli
,
R.
, and
Evans
,
F. G.
,
1971
, “
Relations Among Mechanical Properties, Collagen Fibers, and Calcification in Adult Human Cortical Bone
,”
J. Biomech.
,
4
(
3
), pp.
193
201
.
3.
Currey
,
J. D.
,
1970
, “
The Mechanical Properties of Bone
,”
Clin. Orthop. Relat. Res.
,
73
, pp.
210
231
.
4.
Reilly
,
D. T.
, and
Burstein
,
A. H.
,
1974
, “
The Mechanical Properties of Cortical Bone
,”
J. Bone Jt. Surg.
,
56
(
5
), pp.
1001
1022
.
5.
Knets
,
I. V.
,
Saulgozis
,
Y. Z.
, and
Yanson
,
K. A.
,
1975
, “
Deformability and Strength of Compact Bone Tissue During Tensioning
,”
Polym. Mech.
,
10
(
3
), pp.
419
423
.
6.
Pfafrod
,
G. O.
,
Slutskii
,
L. I.
,
Moorlat
,
P. A.
, and
Knets
,
I. V.
,
1977
, “
Evaluation of Some Deformation and Strength Characteristics of Compact Bone Tissue According to the Data of a Biochemical Study
,”
Polym. Mech.
,
12
(
6
), pp.
934
942
.
7.
Lindahl
,
O.
,
1976
, “
Mechanical Properties of Dried Defatted Spongy Bone
,”
Acta Orthop. Scand.
,
47
(
1
), pp.
11
19
.
8.
Williams
,
J. L.
, and
Lewis
,
J. L.
,
1982
, “
Properties and an Anisotropic Model of Cancellous Bone From the Proximal Tibial Epiphysis
,”
ASME J. Biomech. Eng.
,
104
(
1
), p.
50
.
9.
Goldstein
,
S. A.
,
1987
, “
The Mechanical Properties of Trabecular Bone: Dependence on Anatomic Location and Function
,”
J. Biomech.
,
20
(
11–12
), pp.
1055
1061
.
10.
Ciarelli
,
M. J.
,
Goldstein
,
S. A.
,
Kuhn
,
J. L.
,
Cody
,
D. D.
, and
Brown
,
M. B.
,
1991
, “
Evaluation of Orthogonal Mechanical Properties and Density of Human Trabecular Bone From the Major Metaphyseal Regions With Materials Testing and Computed Tomography
,”
J. Orthop. Res.
,
9
(
5
), pp.
674
682
.
11.
Ashman
,
R. B.
,
Rho
,
J. Y.
, and
Turner
,
C. H.
,
1989
, “
Anatomical Variation of Orthotropic Elastic Moduli of the Proximal Human Tibia
,”
J. Biomech.
,
22
(
8–9
), pp.
895
900
.
12.
Rho
,
J. Y.
,
Hobatho
,
M. C.
, and
Ashman
,
R. B.
,
1995
, “
Relations of Mechanical Properties to Density and CT Numbers in Human Bone
,”
Med. Eng. Phys.
,
17
(
5
), pp.
347
355
.
13.
Hayes
,
W. C.
,
Swenson
,
L. W.
, and
Schurman
,
D. J.
,
1978
, “
Axisymmetric Finite Element Analysis of the Lateral Tibial Plateau
,”
J. Biomech.
,
11
(
1–2
), pp.
21
33
.
14.
Little
,
R. B.
,
Wevers
,
H. W.
,
Siu
,
D.
, and
Cooke
,
T. D. V.
,
1986
, “
A Three-Dimensional Finite Element Analysis of the Upper Tibia
,”
ASME J. Biomech. Eng.
,
108
(
2
), pp.
111
119
.
15.
Mehta
,
B. V.
, and
Rajani
,
S.
,
1995
, “
Finite Element Analysis of the Human Tibia
,”
Trans. Biomed. Health
,
2
, pp.
309
316
.
16.
Muller-Karger
,
C. M.
,
Gonzalez
,
C.
,
,
M. H.
, and
Cerrolaza
,
M.
,
2001
, “
Three Dimensional BEM and FEM Stress Analysis of the Human Tibia Under Pathological Conditions
,”
CMES
,
2
(
1
), pp.
1
13
.https://www.researchgate.net/publication/262485756_Three_dimensional_BEM_and_FEM_stress_analysis_of_the_human_tibia_under_pathological_conditions
17.
Haut Donahue
,
T. L.
,
Hull
,
M. L.
,
Rashid
,
M. R.
, and
Jacobs
,
C. R.
,
2002
, “
A Finite Element Model of the Human Knee Joint for the Study of Tibio-Femoral Contact
,”
ASME J. Biomech. Eng.
,
124
(
3
), p.
273
.
18.
Ionescu
,
I.
,
Conway
,
T.
,
Schonning
,
A.
,
Almutairi
,
M.
, and
Nicholson
,
D. W.
,
2003
, “
Solid Modeling and Static Finite Element Analysis of the Human Tibia
,”
Summer Bioengineering Conference
,
Key Biscayne, FL
,
June 25–29
, pp.
889
890
.http://www.tulane.edu/~sbc2003/pdfdocs/0889.PDF
19.
Tasgetiren
,
S.
,
Verim
,
Ö.
,
Songur
,
A.
, and
Akcer
,
S.
,
2011
, “
3D Modeling and Static Finite Element Analysis of Human Tibia and Fibula Bones
,”
,
Elazig, Turkey
,
May 16–18
, pp.
98
101
.
20.
Ün
,
K.
, and
Çalik
,
A.
,
2016
, “
Relevance of Inhomogeneous–Anisotropic Models of Human Cortical Bone: A Tibia Study Using the Finite Element Method
,”
Biotechnol. Biotechnol. Equip.
,
30
(
3
), pp.
538
547
.
21.
Cattaneo
,
P. M.
,
Dalstra
,
M.
, and
Frich
,
L. H.
,
2001
, “
A Three-Dimensional Finite Element Model From Computed Tomography Data: A Semi-Automated Method
,”
Proc. Inst. Mech. Eng., Part H
,
215
(
2
), pp.
203
212
.
22.
Mayott
,
C. W.
,
Langrana
,
N. A.
,
Alexander
,
H.
, and
Curtis
,
G.
,
1983
, “
Geometric and Mass Properties of Bone as Measured Using Computer-Aided Tomography
,”
,
American Society of Mechanical Engineering
,
New York
, pp.
28
31
.
23.
Woolson
,
S.
,
Parvati
,
D.
,
Linda
,
F.
, and
Arthur
,
V.
,
1986
, “
Three-Dimensional Imaging of Bone From Computerized Tomography
,”
Clin. Orthop. Relat. Res.
, pp.
239
248
.
24.
Marom
,
S. A.
, and
Linden
,
M. J.
,
1990
, “
Computer Aided Stress Analysis of Long Bones Utilizing Computed Tomography
,”
J. Biomech.
,
23
(
5
), pp.
399
404
.
25.
Gosman
,
J. H.
,
Hubbell
,
Z. R.
,
Shaw
,
C. N.
, and
Ryan
,
T. M.
,
2013
, “
Development of Cortical Bone Geometry in the Human Femoral and Tibial Diaphysis
,”
Anat. Rec.
,
296
(
5
), pp.
774
787
.
26.
Viceconti
,
M.
,
Zannoni
,
C.
, and
Pierotti
,
L.
,
1998
, “
TRI2SOLID: An Application of Reverse Engineering Methods to the Creation of CAD Models of Bone Segments
,”
Comput. Methods Programs Biomed.
,
56
(
3
), pp.
211
220
.
27.
Stojkovic
,
M.
,
Trajanovic
,
M.
,
Vitkovic
,
N.
,
Milovanovic
,
J.
,
Arsic
,
S.
, and
Mitkovic
,
M.
,
2009
, “
Referential Geometrical Entities for Reverse Modeling of Geometry of Femur
,”
Second Thematic Conference on Computational Vision and Medical Image Processing
(
VIPIMAGE
), Porto, Portugal, Oct. 14–16, pp.
189
194
28.
Stojkovic
,
M.
,
Milovanovic
,
J.
,
Vitkovic
,
N.
,
Trajanovic
,
M.
,
Grujovic
,
N.
,
Milivojevic
,
V.
,
Milisavljevic
,
S.
, and
Mrvic
,
S.
,
2010
, “
Reverse Modeling and Solid Free-Form Fabrication of Sternum Implant
,”
Australas. Phys. Eng. Sci. Med.
,
33
(
3
), pp.
243
250
.
29.
Marchi
,
D.
,
Maria
,
S.
, and
Tarli
,
B.
,
2004
, “
Cross-Sectional Geometry of the Limb Bones of the Hominoidea by Biplanar Radiography and Moulding Techniques
,”
J. Anthropol. Sci.
,
82
, pp.
89
102
.
30.
Macintosh
,
A. A.
,
Davies
,
T. G.
,
Ryan
,
T. M.
,
Shaw
,
C. N.
, and
Stock
,
J. T.
,
2013
, “
Periosteal Versus True Cross-Sectional Geometry: A Comparison Along Humeral, Femoral, and Tibial Diaphyses
,”
Am. J. Phys. Anthropol.
,
150
(
3
), pp.
442
452
.
31.
Song
,
S. J.
, and
Jeong
,
B. O.
,
2010
, “
Three-Dimensional Analysis of the Intramedullary Canal Axis of Tibia: Clinical Relevance to Tibia Intramedullary Nailing
,”
Arch. Orthop., Trauma Surg.
,
130
(
7
), pp.
903
907
.
32.
Kwak
,
D. S.
,
Han
,
C. W.
, and
Han
,
S. H.
,
2010
, “
Tibial Intramedullary Canal Axis and Its Influence on the Intramedullary Alignment System Entry Point in Koreans
,”
Anat. Cell Biol.
,
43
(
3
), pp.
260
–26
7
.
33.
Harrington
,
I. J.
,
1976
, “
A Bioengineering Analysis of Force Actions at the Knee in Normal and Pathological Gait
,”
ASME J. Biomech. Eng.
,
11
(
5
), pp.
167
172
.
34.
Morrison
,
J. B.
,
1970
, “
The Mechanics of the Knee Joint in Relation to Normal Walking
,”
J. Biomech.
,
3
(
1
), pp.
51
61
.
35.
Duda
,
G. N.
,
Mandruzzato
,
F.
,
Heller
,
M.
,
Goldhahn
,
J.
,
Moser
,
R.
,
Hehli
,
M.
,
Claes
,
L.
, and
Haas
,
N. P.
,
2001
, “
Mechanical Boundary Conditions of Fracture Healing: Borderline Indications in the Treatment of Unreamed Tibial Nailing
,”
J. Biomech.
,
34
(
5
), pp.
639
650
.
36.
Nalla
,
R. K.
,
Stölken
,
J. S.
,
Kinney
,
J. H.
, and
Ritchie
,
R. O.
,
2005
, “
Fracture in Human Cortical Bone: Local Fracture Criteria and Toughening Mechanisms
,”
J. Biomech.
,
38
(
7
), pp.
1517
1525
.