Direct image-based micro finite element modelling of bone tissue

In the paper, a method for the direct image-based creation of the finite element model from images is presented. Image information is taken into account during the calculation of the element stiffness matrix. In this case, material heterogeneity can be included directly in the finite element model. For this purpose, the hypothesis about the correlation between pixel values and elastic properties was used. Four nodes plane element was built. The element can be used with the quantitative phase or scanning electron images and computed tomography data. Simulation for bone data performed. The influence of pixel on the error estimate was studied. The method to parallelize the calculation of the stiffness matrix is presented. As an example, a slice of bone was used in the calculation. Results for average stress distribution for the origin and improved mesh are presented.


Introduction
Nowadays, numerical methods became very popular in different science compartments [1][2][3][4][5][6]. For mechanical simulation, it is important to know mechanical parameters, which can get from experiments. On the other hand, using CT-scans can solve this problem. Value of CT data expresses in the Hounsfield scale. Nowadays it is common knowledge that Hounsfield values related to different kinds of biological materials. It was shown that Hounsfield values related to bone mineral density, elasticity, and critical stress. Homogenized properties can get from the material structure, using a quantitative phase or microscopy images or μCT. Of course, in this case, the properties of components should be known. This approach allows getting information about complex anisotropic, fine-grained or porous materials, as well as other structures characterized by heterogeneity and can be used for homogenization of the media. There are various methods to get homogeneous models using CT data.
Various methods of homogenization include representative volume elements [7][8][9][10][11][12][13], but in this case number of calculations should high. Another side there are methods to describe the structure of media using e.g. fabric tensor [14][15][16][17], in this case, specific functions for material should be founded. The development of numerical technologies allows improving classical FE methods for such tasks [17][18][19][20]. The purpose of the study is to present a finite element for describing the mechanical behavior of objects using it's CT data.

Materials and methods
The study is to get the mechanical behavior of a system filling some continuous region Ω and connected discrete data Ω'. Studied object V lies in the Ω. Corresponding CT data exists. The CT data consists of a number of micro-volumes with values of radiodensity, which can be visualized using colors. To create discrete data V' from Ω' the region V should be split by a large number of micro-areas (pixels). The size of a pixel is described by ∆x 2 , coordinates of a pixel are described by xk, yk. Pixel size depends on measuring device resolution. On figure 1 the example is illustrated. V consists of material and pores are supposed. The material is solid and isotropic, but the presence of pores makes the material anisotropic on a macro scale. CT should be binarized, for this purpose different methods can be used (e.g. Otsu method), but it is important to understand the correlation between Hounsfield scale and material parameters [21][22][23]. There are continuous region Ω and connected discrete data Ω' shown in figure 1. A full description of the method of the three-dimensional case has presented the paper [23]. Let's give key points of the method in case of plane strain. Consider the FE realization principle of a heterogeneous medium. Well known method [24][25][26][27][28][29] of constructing a fournode plane isoparametric FE with a linear geometry approximation and the displacement field in local coordinates was used. For the approximation of the radius-vector r (geometry) and the displacement vector ϑ such functions were used: where Nn(ξ, η) = (1+ ξnξ)•(1+ ηnη)/4 -shape functions, ξn, ηn -local coordinates of element nodes, u, v -the displacement vector projections on the orts of the global Cartesian coordinate system x, y. The relations can be rewritten in matrix form: where [N] -matrix of approximating functions, ϑe -vector of nodal displacements. The deformation of the medium can be described by the components of linear deformations and shear deformations, which are represented as a reduced vector ε and expressed through displacements by known relations: xx = , ε yy = , ε xy = + , ε zz =ε yz =ε zx = 0, (4) which can also be written in matrix form: The stress state is determined by the stress tensor through the components of normal σxx, σyy, σzz, tangential σxy, σxz=σzy=0 stresses. The stress tensor, as well as the elastic strain tensor, can be written as a reduced vector σ. Hooke's law relating the reduced stress and strain vectors are presented in the form: The well-known relations are used to calculate the local stiffness matrix [29]: where V e volume of FE in Ω, [D(r)] is the elasticity matrix of the body, and [B(r)] is the matrix connecting the reduced vector of elastic deformations and the vector of nodal displacements [23,29]: Using the hypothesis that anisotropy of media appears because of the porous, local stiffness matrix of FE in Ω' can be calculated by means of the equation: where V e volume of FE in Ω', ω(r) -function Ω → Ω' and [D(r)] is the elasticity matrix of the isotropic material and can be described by Young's modulus and the Poisson's ratio.
In subsequent equations integration in z-direction will be put down. Using to the local coordinates, the equation can be rewritten as: where |J(ξ,η)| is the Jacobian determinant of coordinate transformation, where h edimension of FE in z-direction.
To integrate local stiffness matrix the method of central rectangles was used. The FE of a heterogeneous medium is a convex quadrilateral with linear edges. The quadrature points for such element are the coordinates (local coordinates) of the pixel centers from the image data, then: where: ξi, ηilocal coordinates of integration points, ∆ξ, ∆η -subinterval's step value in two directions in local coordinates, ω(ξi,ηi) -weights of the quadrature formula determined by the image data value at the integration point, I, J -the set of quadrature points for each local coordinate in the element.
To check is pixel inside FE the following inequality was used: where rA -radius vector of A point, rI -radius vector of I point, ni -normal vector to edge number I (i from 1 to 4) of the FE, δ -accuracy value. On figure 2 digital data and FE are shown. IJKL are nodes of the element. In case when quadrature point A is fully inside FE (green squares on figure 2) δ can be equal to zero. Some problems occur during integration because of the discrete origin of the image data. It is usual problems like for quadrature points on edges (blue and orange squares in figure 2). For these cases, the neighborhood of the point was studied. Distance between quadrature point B and FE's edge can be calculated and pixel's volume can be recalculated. Figure 2 illustrated options with the following criteria: if the center of a pixel inside FE it included in integration (orange squares) otherwise not (blue squares). The problem of the meshing area goes beyond the paper, but nowadays there are some methods to create mesh on CT data [30][31][32].

Fig. 2. Schematic view of binarized CT data (left) and finite element and it's quadrature points (right): green squares certainly lies inside FE and they are included in integration, orange squarescenters of squares lies inside FE it included in integration too, blue squares -centers of squares do not lie inside FE and they are not included in integration.
In the assumptions stress results in node are invalid. For analyzing stress state average stress should be used. Average stress can be calculated by expression: where V -an averaging volume of FE in Ω, V' -a volume of FE in Ω'. Averaging can be performed in the whole FE or in subvolumes of FE. This procedure can improve understanding of the stress-strain state in macro volume, which is important in problems of bone tissue growth and prediction fracture [33][34].

Results and discussion
The described method was implemented for four-node FE with two degrees of freedom (in the corresponding Cartesian coordinate axes) in each node. C++ program was implemented for calculations to be done. For solving problems computer was used with the following complete set: processor -AMD Ryzen 7 1700 with 8 physical and 16 logical cores with a frequency of 3.7 GHz., memory -G. Skill Aegis 16GB DDR4 3000 16GISB K2 C16 with a frequency of 2933 MHz., motherboard -MSI B350M MORTAR. Because of the high labor costs of integrating the local stiffness matrix, the algorithm was parallelized using OpenMP technology: local stiffness matrix integration for each FE was performed on different threads. The results were visualized in the ParaView software.
For numerical problems, CT data of rat's femur diaphysis was used. Micro/nanofocal Xray inspection system for CT and 2D inspections of Phoenix V | tome | X S240 was used for scanning. The system is equipped with two X-ray tubes: microfocus with a maximum accelerating voltage of 240 kV power of 320 W and nanofocus with a maximum accelerating voltage of 180 kV power of 15 W. For primary data processing and creating a volume (voxel) model of the sample based on x-ray images, the datos|x reconstruction software was used. The sample fixed in the holder was placed on the rotating table of the X-ray computed tomography chamber at the optimum distance from the X-ray source. The survey was carried out at an accelerating voltage of 90-100 kV and current 140-150 mA, the dimension of the voxel was 6.747 μm.

Fig. 3. Origin CT data and FE mesh.
For the test problem, all degrees of freedom for nodes on the inner perimeter were fixed and compressive pressure was applied on the outer perimeter. In figure 3 FE mesh and CT data are shown. The following mechanical properties of bone tissue was used: Young's modulus equal to 20 MPa and the Poisson's ratio equal to 0.3. Mesh was generated with three elements in cross-section. On figure 4 results for average stress are shown.  The nodal approximation is worse on the rough mesh, that's why averaging stress not smooth (see figure 4). To specify the location of maximum stress averaging can be carried out on fine mesh. Grinding the mesh improves the smoothness of stress distribution. On figure 5 average stress distribution is shown. Voxel density in the FE influence not significantly on stress result and this fact allows improving mesh for aver-aging without losing the accuracy. It can be done easily by remeshing every FE. Figure 5 results on the fine mesh are shown. It can be noticed that the maximum values of stress differ (8% -25%). It can be explained by decreasing averaging volume and it should be understood as clarification of stress distribution.

Conclusion
This paper presents a version of the finite element method with associated CT data and is developed to solve CT data-based structural problems. It is based on a special technique of integration of the FE local stiffness matrix and allows it to takes into account heterogeneity of the material. The method is based on the assumption that anisotropy of material is based on heterogeneity and that correlation between pixels values and material properties exists. The additional value of the method is averaging FE the results in subvolumes, which can increase the quality of results distribution and simplify analysis of the stress-strain state.
The work was partly supported by the Russian Foundation for Basic Research within the scientific projects № 20-01-00535.