An algorithm for finding and adding boundary conditions with the aim of solving the contact problem in computational mechanics

— In computational mechanics, the contact problem between moving bodies is always an interesting topic. The reason is simple: there is no unique solution for contact problems that works in any case in terms of arbitrary geometry, applied boundary conditions, external loads, material properties, etc. In addition, the process of finding parts of the outer surfaces of the bodies in contact is always a programming challenge. When we keep in mind the fact that numerical programs are written mostly in Fortran programming language, which does not have as powerful tools as some modern languages do, the level of the problem becomes clearer. Of course, there is the possibility of writing such algorithms in some up-to-date programming language with powerful libraries for geometry, searching, etc., but in that case other problems arise, related to code compatibility, data transfer, execution speed, parallelization, etc. This paper shows an algorithm for finding contact entities between two or more different bodies. The algorithm is implemented within the PAK software package for finite element analysis [ref]. The PAK software is written in Fortran and uses its own skyline solver for solving of system of linear equations, but also has MUMPS parallel sparse solver [ref] incorporated in itself.


I. INTRODUCTION
Numerical software based on finite element method can be used for solving of different types of partial differential equations: Laplace's equation, wave equation, Maxwell's equation, Helmholtz equation, Poisson's equation, etc. Those equations, with application of variational methods, are commonly used for solving of different types of scalar or vector field problems such as electrostatics, electromagnetism, sound wave propagation, diffusion, etc. Similar to that, in solid mechanics, using principle of virtual work and constitutive equations, unknown displacements vector field can be solved for prescribed initial and boundary conditions. In all cases, the system of differential equations is reduced to a system of linear equations by application of the numerical integration methods.
This concept is common and solves many different problems. However, it also imposes one significant limitation. To solve any complex problem, the domain of the problem must be fully connected. This means, for example, if we solve the problem of the propagation of a wave through some medium, it must be completely connected, that is, the waves will propagate to all the boundaries of the finite element mesh. If the considered domain consists of two separate finite elements meshes, the propagation of the waves would not manifest from one mesh to another, although they can be in geometrical contact. The equations in the nodes belonging to different meshes are not connected, although the nodes may be in the same geometric position. Similar to the previous example, if we analyze the problem of interaction between two solid bodies, where the boundary condition can be prescribed displacement, velocity, acceleration or force on a body, bodies are not aware of each other's existence. Their mutual interaction must be ensured by appropriate additional equations, boundary conditions or modification of existing equations [ref]. Since the movement of the body takes place over time, at each time point the bodies are in a different configuration and spatial position. This means that for each time point the relative positions of the various bodies should be examined. For the boundary nodes of the finite elements mesh, which are found to be in contact with the nodes of another body mesh, appropriate action should be performed: adding a new equation, adding of new boundary condition or modification of existing equations [ref].
Finding the finite element nodes of different bodies that interact with each other is one of the demanding tasks that needs to be done in order to enable the mechanism of contact between those two bodies. This paper presents one of the ways how this can be done using geometric entities such as planes, vectors, points, polygons, etc.  Figure 1 shows simple explanation about the problem. There are three different bodies represented by three independent finite element meshes: two simple plates and wired structure placed between them. This example is a An algorithm for finding and adding boundary conditions with the aim of solving the contact problem in computational mechanics Velibor Isailovic, Nenad Filipovic simple test that is used in medical stent development. The test procedure consists of compressing of the stent between two flat plates and measuring displacement of characteristic points and contact forces as a test output. However, the test cannot give any information about the state of stress or deformation within the material. For this purpose, numerical simulations are very often used to support mechanical tests such as one mentioned before. Some facts and conclusions can be determined from the simulations that cannot be seen by experiment. For example, state of internal forces, stress or deformation, etc. But, in order for such a numerical test to be possible, it is necessary to develop finite element algorithm that allows contact between the two bodies. Otherwise, the bodies can pass through each other without any limitations. Such an algorithm requires a mutual relationship of the bodies as input information, i.e. which nodes (or outer faces) of one body are in contact with the nodes (or outer faces) of another body. Finding and calculation of that information are the topic of this work.
As mentioned in abstract section, all the work was done in the Fortran programming language, which was a challenge in itself. The main reason why it was done in this way is because the rest of the code related to the finite elements is written in the same language. For the purpose of solving examples with a huge number of nodes and finite elements, the code is adapted for running on supercomputers by using MUMPS parallel solver [ref] [ref].

II. METHODS
In order to provide information about the mutual contact of the body, it is necessary to have information of the outer polygons, i.e. the polygons that make up the outer surface. Known data are geometric positions of points in space and polygons with a corresponding list of points (ID of points).
In general, there are two ways how the contact between two bodies can be recognized. The first way is to find a node of one body that penetrates the boundary surface of another body. The second way is very similar, but instead of node, the center of gravity of face is observed. In that case the contact is achieved if center of gravity of face from one body surface penetrates the boundary surface of another body. Both concepts are similar in implementation; the difference is only in set of points that should be tested. where the point C of the body 1 is placed inside the body 2. This is the simplest case of contact between two bodies. In order to define if point C is in body 2, we should examine their relative positions. For cube and point it is pretty simple: based on the angles between cube sides and coordinate axes, we can calculate angle of rotation of cube so that it would take a position parallel to the axes, rotate cube and point and based on new coordinates determine whether the point is inside the cube or not.
However, the finite elements can have arbitrary tetrahedral or hexahedral geometry. In computational mechanics hexahedral elements are used more often than tetrahedral elements due to accuracy of calculation of postprocessing variables like stress or strain. The paper deals with hexahedral elements but this is also applicable to tetrahedral elements. This kind of elements has eight nodes, and their position in space is such that the nodes of one side of element do not have to lie in the same plane ( Figure 3).

Figure 3: Hexahedral finite element with arbitrary geometry
For the element with such shape it is not possible to determine whether the point is inside based on coordinates.
To solve this problem, we must define a new curvilinear coordinate system with the origin at the center of gravity of the element. The axes of this coordinate system follow the geometric orientation of the element. For such a coordinate system, the values of the (r,s,t) coordinates are in the range from -1 to 1. The criteria for determining whether a point is within an element or not is the value of the coordinates in the local curvilinear coordinate system. If the values are in range [-1, 1] the point is inside, otherwise the point is outside. The algorithm that we used for finding the values of coordinates is multivariate Newton -Raphson method [ref]. When a point that falls inside a boundary element is found, the contact pair (node to face or center of face to face) is defined.
It is now possible to calculate the distance from the contact node to the contact face. This can be done by placing a plane through the face nodes. Since faces consist of four nodes that do not have to lie in plane, we define plane by calculation of two normal vectors by dividing a quadrilateral face into two triangles. The mean value of those two vectors is resulting normal vector. It doesn't matter which product is taken. All cross products give the vector with the same direction. Since this vector is normalized, the choice of three nodes is completely arbitrary.
In the second case, when face nodes do not lie in a plane, face normal vector is calculated as mean value of two normalized vectors, obtained from two triangles made by dividing a quadrilateral face.
The resulting vector is obtained as the normalized sum of these two vectors: The way the quadrilateral face is divided does not affect the direction of the resulting normal vector. When we have gravity center of face (the point with mean values of face nodes coordinates) and face normal vector, we can define the plane for calculation of distance from the boundary node of the same body, that penetrates the observed face, and also normal projection of mentioned node on face plane. Those data are necessary for setting up elastic supports, in order to enable contact between the two bodies. The plane parameters can be defined using calculated normal vector and center of gravity: Where: face center of gravity, any point on the plane.
A plane defined in this way allows us to calculate the depth of entry of a point of one body into another body. The depth is calculated by: where are plane parameters obtained from equation [eq] by substitution of face normal and gravity center. It is now possible to calculate the stiffness and force values necessary to return a node that has breached the boundary of another body to its boundary. Figure 3 shows the mutual relationship of the two bodies in contact. During body movement, nodes of one body may pass through the boundary surface of another body. In this simple example both bodies consist of only one finite element. The node A passed through a shaded face. When this happens, the node needs to be returned to the boundary of the body. This is done by projecting node A on the face plane defined by normal vector calculated using of equation [eq]. The criteria for defining the contact between two bodies are:

Figure 5: Contact between two bodies: The node from one body entered another body
1. Boundary node of one body is located inside a finite element of another body; 2. Projection of mentioned boundary node on the face plane of another body falls inside the face.
The first condition must be satisfied in order to test the second one. If the first condition is satisfied and the second one is not, there is the first special case, described in Figure  6. In order to give simpler explanation, the figure is given in two dimensions.

Figure 6: Relation of the nodes and elements that cause special case
As can be seen from Figure 6, it is possible for boundary node from first body to be inside an element of another body while the projection of that node on the plane of the element boundary face is completely outside the element. In that case, an adjacent face must be found for which the projection of the node on the plane of the face falls within the face itself.
In the example given in Figure 6, node N is located inside element E2 but its projection on the plane of the boundary face is outside the element. However, the projection of the node on the plane of the boundary face of the element E1 is inside that face, so it is the geometric position where the node N should be moved. The second special case is shown in the Figure 7.

Figure 7: Second special case: Node N is inside
This case appears when the node N falls inside element E2, which does not have any boundary face. If such an element is found, it is necessary to find all its surrounding elements that have boundary faces, and calculate the mean value of the normal vector of all boundary faces. In the case shown in Figure 7, the normal vector would be calculated as the mean of the normal vector of the boundary faces of the elements E1 and E3. The situation can be more complicated in case of 3D problems, because element like E2 can be surrounded by more than two elements, but the way the normal vector is calculated is similar. By visiting all boundary nodes of one body and checking whether they enter a boundary element of another body (including neighbors of boundary elements due to special cases) we can find appropriate relations and calculate directions where to move boundary nodes that broke the boundary of a neighboring body, based on calculated normal vectors (or mean value of normal vectors in special case two).

III. ELASTIC SUPPORTS
Now, for each boundary node we have two pieces of information: 1. Information whether the node is in contact with another body or not; 2. Value of normal vector in node which is in contact with another body.
Based on these information, the stiffness matrix for elastic supports can be calculated. That stiffness matrix can be calculated by using theory for one-dimensional finite elements [ref].

Figure 8: Elastic supports placed in boundary nodes that broke outer boundary of neighboring body
In some cases, depending on the geometry, it may be more appropriate to observe the relations of the center of gravity of the face to the outer faces of another body. The contact algorithm is very similar to previously explained. In this case the positions of the center of gravity of faces are calculated, based on their corner nodes. Other searches and calculations are the same as in the previously explained algorithm.

IV. RESULTS
As a result, finite element model of medical stent is presented here. The medical stents are wired structures used in treatment of cardiovascular diseases caused by narrowing of blood vessels. During stent development, a large number of experiments are performed to achieve optimal stent design. These experiments are quite expensive and require a fairly long period of time for testing. They are usually performed until complete failure of stent or by applying of cyclic load with number of cycles calculated based on heart rate and the number of years the stent should survive in the human body. Computer modeling plays a very important role here in stent development. For instance, in experiments, it is not possible to see what happens in the stent materialthe state of internal forces caused by external loading is not available information. This information can be obtained by numerical tests using finite element method with contact algorithm. In that sense, development of medical stents can be significantly improved by computer simulations.
In order to be able to make a virtual experiment, it is necessary to provide a contact algorithm to be able to simulate the interaction of the stent with the rest of the equipment for mechanical testing of the stent, i.e. the interaction of the finite element mesh of the stent with the finite element mesh of other pieces of testing equipment.
There are several different tests where stents are tested with pressure loading, cyclic loading, passing stents through tubes that simulate blood vessels, etc. For any type of test, it is possible to make a numerical simulation, which would provide a detailed insight into the state of stress or internal forces in the stent itself.
The aim of this paper was to show the way in which all the necessary data used to modify the system of finite element equations are provided. Figure 9 shows the test example. An example consists of bending of stent between three cylinders. This is the so-called three-point bending test. The boundary conditions of finite element meshes are: 1. Two fixed cylinders placed bellow the stent; 2. One cylinder with prescribed displacement placed above the stent; 3. Three contact surface pairs between stent and three supporting cylinders. Figure 9: Three point bending stent test: The stent is exposed to a bending load using two fixed cylinders and one movable cylinder. Figure 10 shows detail of moving cylinder in the different time steps during the simulation. At the beginning of simulation there are no contact elements between cylinders and stent. During the movement of the upper cylinder, the contact first appears between the upper cylinder and the stent and then between the stent and the lower cylinders. The number of contact elements increases during the simulation, which is visually shown in Figure 10 and also in Figure 11.  This paper presents an algorithm for finding contact pairs between finite element meshes of two bodies. We proposed a way how to find information about contact between different bodies modelled by finite element method. Based on that information, system of finite element equations is modified by adding elastic supports, which ensure that a body does not break the boundary of another body.
The source code is written in FORTRAN 95, to ensure maximum compatibility with numerical software written earlier. Project is implemented as a separate module, which can be applied completely independently from the remaining code. The source code is available to all interested scientists who want to use it. By calling the appropriate subroutines, it can be easily adapted to any finite element code, since the input to the module is information about finite element mesh and displacements of all nodes, and the output is the position, spatial orientation and values necessary to calculate the stiffness of elastic supports and adding them into the existing system of finite element equations.
In results section, numerical model of stent subjected to the three-point bending test is presented. This is just one example where a developed contact algorithm can be applied. There are many other examples and areas of application, for instance: numerical simulation of a car crash test, plastic deformation processing, tire rolling, gear or belt transmissions, contact of bone and cartilage, kicking the ball in various sports, etc. The presented algorithm is general and can be applied in any of the stated fields of research.
For modeling of more complex examples, it is possible to introduce various nonlinear material models, in order to make a more realistic model, with as few approximations as possible. It is also possible to use structural finite elements such as a shell, beam, plate, membrane, etc. in order to reduce the number of degrees of freedom and simplify the model.
The numerical software presented here was created as a result of research on international scientific project InSilc 1 . Software solution will be uploaded in public repository provided by European Commission. In this way, it will be available for whole research community interested in the field of contact problems modeling in computational mechanics. Since there is no unique software solution for contact problems that provides a solution for each specific problem, and especially not an open source solution, this software solution can be very interesting for further upgrade and application on problems from other engineering fields.