Contact phenomena modeling in biological structures on the example of the implant-bone

The work presents modeling of interface phenomena in biological structures. Selected ways of numerical modeling of phenomena at the boundary of two materials by means of FEM methods are discussed. The work focuses on phenomena related to biological structures and their mutual interactions. Using an example of an implant-bone system, various techniques of modeling interface phenomena are compared and referred to experimental results. The study reveals the main features of the selected modeling techniques, e.g. complexity of creating the model, time-consumption of computation, reliability of the obtained results. The obtained results proved that the most advantageous method is the one that makes it possible to regard interactions as numerical values, without an excessive generation of the finite elements (the grid) at the materials boundary. The contact modeling technique based on the contact indicator showed the lowest value of standard deviation in relation to contact modeling techniques: with the use of elastic-damping elements, based on the so-called boundary slip line and based on the so-called bio-layer. The contact indicator should be understood as a numerical value describing the interaction of the boundary surfaces of two objects cooperating with each other. The value of the indicator was determined experimentally. The value have oscillated +/0.0012mm. This is a very small value knowing that the standard deviation for the modeling technique of the contact employing damping-elastic elements gives a spread of +/0,1442 mm.


INTRODUCTION
Modern medicine relies on engineering solutions in order to create procedures of a healing process in a more precise way. The engineering solutions can be found first of all while planning a reconstruction of an injured human motor system. These works generate a number of biomechanical phenomena, which should be modeled in an engineering way before the surgery, so that the phenomena occurring in the system under consideration can be analyzed (Uklejewski et al. 2008, Schmolz et al. 2000, Rocco et al. 1999. Results of experiments related to biomechanical properties of living structures having an artificial material implanted inside, are widely used in research pertaining to surgical fixation of bone fractures as well as in processes of repairing damaged joints using bone implants (Alonso et al. 2000, Billsten et al. 2000, Childs et al. 2002, Francis et al. 2000, Harpal et al. 2014, Gregory et al. 2011. Accuracy of the obtained results is of a crucial importance with regard to a future surgery. In order to eliminate risks connected with a failure of a surgery, an initial study is performed by the way of launching computer simulation, mainly based on FEM techniques. The form and complexity of the created model directly influence accuracy of the obtained results. It is striven for a simplicity of creating the computer model of the biomechanical phenomena that occur in biological structures, including the interface phenomena; yet, at the same time, the obtained results should reconstruct the real phenomena that occur in the living structures, making it possible to simultaneously analyze particular medical cases. The aim is to make a computer model of biomechanical phenomena occurring in living structures. contact phenomena was easy to work while giving results that most closely reflect the real phenomena occurring in it. An easy description of these phenomena will allow you to accurately analyze various medical cases. The fist approach to modeling a contact of deformable bodies was undertaken by Hallquist in 1985 and then by Wriggers in 1990, using convective coordinates connected with the contact surface, after regularization of the boundary conditions by a method of enhanced Lagrangian multipliers (Neto et al. 1985). Among the problems that are of crucial importance for the development of medicine and related sciences, the contact problems belong to a group of difficult problems (Kovacs et al. 2003, Kusserow et al. 2015, Parkash Srivastava et al. 2017, Carniel and Fancello, 2017. A high degree of complication of the contact problems in the case of systems containing a living structure, results from many factors, including nonlinearity (both geometrical as well as of constitutive models), material properties of a living structure (featuring orthotropism or materials regarded viscoelastic), problems that occur at the discretisation level (Sumner et al. 1992, Smith et al. 1998, Mars et al. 2017Zaczyk, Jasińska-Choromańska, 2011). The present increase of computing power and development of such computer techniques as FEM results in the fact that solving these difficult problems became more effective and is used more and more often while dealing with biomechanical problems. It was attempted to show a practical application of computer FEM methods, using a simplified modeling of contact problems in simulations that yield results comparable with outputs of simulations based on complicated contact models regarding most of the contact conditions related to the real object, among which we can name kinematics of the contact, friction occurring between the bodies, measure of elasticity of the bodies, constitutive connections, magnitude of deformation of these bodies.

MATERIALS
The presented considerations are related to simplified techniques of modeling interface phenomena in a medium consisting of two different materials. One of them has the features of a biological structure (it has the ability to grow, remodeling and adapt because of the surrounding conditions and the ability to repair) and the other is an artificial material with constant physical properties. An example of such medium is an implant-bone system (Figure 1), where the bone has traits of the living structure and the implant corresponds to the artificial material. Another example may be a connecting element mounted within a tree (Figure 2), where the tree has traits of the living structure and the connecting element in a form of a screw or a nail plays the role of the artificial material. In the considered medium, at the boundary of two materials, there take place a lot of phenomena of mechanical and biological character. Some of them are as follows: • phenomena taking place at the boundary between an isotropic and orthotropic body, at a considerable spread of the material properties, • phenomenon of occurrence of large displacements within the medium at the boundary of different materials, • phenomenon of occurrence of large displacements, accompanied by deformation of the surface that is more deformable -in the considered case, it is the surface of the material that is mechanically weaker, • occurrence of friction phenomenon (Coulomb) at the boundary of two materials -one isotropic and the other of a orthotropic character, • interaction of a living structure with an artificial material on the border line between static and kinetic friction, • possibility of destruction of the weaker element, • possibility of delamination of the medium at the boundary of two materials, • possibility of dynamic strain and adaptation of the living structure in a biological way to interactions of the external conditions.
The aforementioned phenomena illustrate complexity of the mutual interaction of living structures with an artificial material. Despite many options while modeling interface phenomena by means of engineering tools for FEM simulation, it is difficult to find a technique, which would make it possible to precisely, yet in a simple way, regard all these phenomena taking place within the considered medium. Limitations due to the accepted technique of solving the task may generate some problems connected with the following: • large disproportion in values of so-called material constants describing strength of particular materials, the medium is composed of, • possibility of occurrence of large displacements at the boundary of finite elements and breaking continuity of the model structure, what makes it necessary to introduce a continuity index of the model in order to ensure appropriate discretisation of the whole model, • possibility of a degradation of the weaker element -in case of a contact, where the finite elements of the weaker material overlap each other, • possibility of occurrence of a local interpenetration of finite elements within the medium, • possibility of delamination of an element -in case of a contact, where a local breaking of the continuity of the model takes place, what results in appearance of a clearance, a gap between elements describing the living structure and elements describing the artificial material.
A limitation in FEM techniques is also the assumption that in the contact algorithm, both contact surfaces are treated equally. In the calculations, the contact is treated as a reference of one surface to another, eg the wall of the first component acts on the wall of the second component or vice versa. A contact that takes place between the biological structure and the artificial material is a problem, where inequality conditions appear both at the level of geometry, material properties, as well as the structure of the contact surfaces. Because of searching for the most simple form of modeling the contact between living structures and an artificial material, it was accepted that the most convenient form was a contact realized between the surfaces, disregarding at the same time a contact at the nodes level. Creation of the contact at the level of the nodes requires to perform a toilsome work that consists in isolation of the boundary nodes and assigning for them appropriate mating parameters. The following table presents selected techniques of modeling the contact that were analyzed. Contact with the Contact employing damping-elastic elements allows you to model surface contact by limiting the number of degrees of freedom and giving damping in the tangential direction of the friction force. Contact with the on socalled slide lines is characterized by the possibility of defining the friction in the tangential direction and the interaction of the element subordinated to the main element. Contact based on the socalled the bio-layer defines the contact by means of the friction force in the tangent direction assuming an even interaction between the contact elements by means of Gauss points. The contact based on the contact indicator is characterized by the definition of the coefficient for the tangent of the friction force at the boundary of the two elements. This condition assumes that the contact elements are in point contact. Detailed descriptions are included in Table 1.  The techniques of interface modeling presented above were selected in such a way, as to make them applicable in any engineering tool containing a FEM module, either standard or premium one.

METHOD
An analysis related to capabilities of a simplified modeling of interface phenomena for biological structures being in contact with an artificial material is presented. An aspect of the analysis concerning the modeling techniques was complication of proceedings while modeling as correlated to quality and reliability of the obtained results. The analysis consisted in a qualitative evaluation of realization of the proceedings versus amount of the necessary input data. A complement of the performed simulations was an experiment on real objects, where the considered contact took place. The analyses and the experimental study were performed for interface phenomena taking place at the boundary between the implant mandrel and the bone, accepting that the bone is the living structure and the implant is the artificial material. Accepting the contact at the boundary between the implant mandrel and the bone seems to be justified since it is one of the most difficult to be modeled. The difficulties result from an asymmetry of the osseous structure, diversified material traits and the related biological phenomena. Material properties of particular elements of the model are listed in Table 2.  The works consisted in elaboration of a geometry of a 3D model in a form of a "solid", having a structure of the cortical and trabecular bone specified, and in simplification of the implant structure. Then, a grid of finite elements of "hibrid" type was generated over the model prepared in this way. Such approach eliminates many errors while creating the grid in objects featuring a complex geometry. The boundary conditions as well as material properties for the models were the same in all the considered cases.
Boundary conditions of the created model corresponded to a compression force applied to the artificial-materialbiological-structure system during its standard operation. The structure system are shown in Figure 3. A lot of computer simulations were launched for the model created in such way, where only the technique of modeling the contact at the boundary of the living structure and the artificial material was changed, applying a modeling technique according to Table 1. Results of the performed experiment were accepted as the reference conditions for the analyses carried out, which were related to various techniques of modeling interface phenomena using FEM. The experiment consisted in selecting such biological structure, whose geometry, mechanical properties and character of operation were similar to the proposed contact model. An effect of the selection was preparation of a real object for the artificial-material − biological-structure system, with particular attention paid to creation of the contact at the interface of particular elements of the studied object and determination of their mechanical properties, e.g. Young's modulus and Poisson ratio. The biological-structure system are shown in Figure 4. The elaborated physical object was subjected to a quasi-static compression test. During the test, there were directly recorded displacements, which occur between the artificial element and the biological structure of the tested sample, along with the respective applied loads.

RESULTS
The performed computer FEM simulations as well as experimental study were focused on determining value of displacement of the artificial material (i.e. fragment of the mandrel) with respect to an adjacent fragment of the living structure (i.e. bone). Such decision resulted from the fact that it was possible to determine this parameter directly, both in the computer simulations as well as in the experimental study. Such approach made it possible to eliminate errors resulting from determination of parameters indirectly. The following table 3 presents values of displacements obtained owing to simulations based on models implementing different techniques of modeling the interface phenomena at the boundary between a living structure and an artificial material.

STATISTICAL ANALYSIS
All the simulation studies were launched using the same object featuring the same number of finite elements within the model, being equal to 5552, and number of nodes of 12508. These data were obtained from the preprocessor of the FEM simulation tool. The original contact surface at the boundary of the living structure and the artificial material was the same in each model. For each of the created models, a number of simulation was launched, changing only value of the load in the boundary conditions and recording displacements of the elements describing the biological structure with respect the artificial material. Parameters accepted in the model (described in Table 1) were determined experimentally, investigating properties of the bone structure on the basis of density of the osseous tissue. In order to determine properties of the living structures, 18 real objects were studied for this purpose. The real density of the cortical bone was of ca. 1.9 g/cm 3 , and density of the trabecular bone from 0.14 to 1.1 g/cm 3 . On the basis of the determined real density of the osseous tissue, an replacement Young modulus was determined during the experiment. The equivalent Young modulus of the osseous preparation was determined using Keller equation. The Keller equation enables the substitution of Young's modulus for these structures based on the actual density (g / cm3) of the bone tissue (Rho et al. 1995).
A basis for determination of the other properties of the bone structure were experiments performed by M. Nordin. For particular measured quantities, their standard deviations were scaled down/up. Results of these measurements were processed using the Mann-Whitney U test, in order to determine the absolute value of parameters, which were used for determination of material conditions while creating the simulation models.

CONCLUSIONS
The performed analyses and experimental study prove that despite the fact that there are so many possibilities of modeling the interface phenomena, yet universal tools implemented in the software for FEM simulation do not reconstruct in a satisfactory way the phenomena that take place in the real objects. The considered description of the contact at the boundary of a biological structure and an artificial material is connected with significant spread of results, whereas universal models reveal an increasing spread with respect to the reference point, created for the models. It turned out that among the analyzed contact models, the most convenient one was the model based on biolayer and the model with an assigned value of the contact index of a living structure and an artificial material. Models based on elastic elements and so-called slide lines revealed the biggest deviations with respect to reference values obtained experimentally. The obtained results show that using special contacts brings many benefits. It eliminates the need to compact the FEM mesh and the problems associated with it. extended calculation time. The degree of confidence in these solutions can be assumed to be high. The contact modeling technique based on the contact indicator showed the smallest value of the standard deviation. The value have oscillated +/-0.0012mm. This is a very small value knowing that the standard deviation for the modeling technique of the contact employing damping-elastic elements gives a spread of +/-0,1442 mm. The conducted research has shown that the standard deviation for the outcomes from simulations based on special contacts (contact based on socalled bio-layer and contact based on contact index) constitutes 3% of the values sought.

DISCUSSION
One attempted to find a simplified method of modeling the interface phenomena at the boundary of a living structure mating with an artificial material, concluding that it is possible to use the existing algorithms implemented in the software designed for modeling the contact. However, only such algorithms should be selected, which implement such description of the contact that is the closest to the phenomena taking place within the real object. Biological structures themselves are problematic while modeled because of their complex geometry, variable local biomechanical properties and a possibility of manifold reaction to the conditions of external interactions (Neto et al. 1985, Uklejewski et al. 2008, Sumner et al. 1992, Zaczyk, Jasińska-Choromańska 2011, Zaoui, Hadjuoui. 2016. The accepted simplifications in the geometry result in smaller errors than in the case of an inappropriate contact model. Thus, a key issue is description of the contact itself. As results from the performed analyses and experiments, dedicated techniques of modeling the contact in the case of biological structures mating with an artificial material, by means of such methods as bio-layers algorithm or the contact index determining interactions of the material, increase reliability of the obtained results. A hardship related to this modeling method is a necessity of experimental determination of parameters describing mutual interaction in order to determine contact indexes or a relevant algorithm. Despite the difficulty, an advantage of the methods is a possibility of eliminating a compression of the grid at the boundary of the contact and thus shortening the computation time, while the process of creating the model is aimed at reconstruction of the geometry and setting simplified boundary conditions. Despite the simplification while modeling many phenomena of a biomechanical character related to biological structures, it is still possible to obtain correct results using universal tools for modeling the contact.