This chapter is organized in the following manner. Modelling of an orbital cavity using finite element method is presented in section 2. Finite element method (FEM) is one of the basic tools used for mechanical investigations of a skull, a pelvic bone, eye-socket and in reconstruction of a bony face deformed by congenital defects or injuries. For example, it was applied for modelling of a skull with gnathothisis (Boryor et al., 2008), for investigation of infant head injuries caused by impact (Roth et al, 2009) and for examination of facial skull dystosis of a child (Gautam et al, 2007). FEM served also for an radiological and mathematical analysis of facial deformation provoked by curved central axis of the skull (Iannetti et al., 2004). Furthermore, the method was utilized for modelling of orbit deformation being result of bunt injury (Al-Sukhun et al, 2006) and also for investigation of biomechanical properties of the orbit (Sander et al, 2006).
The aim of the study is to develop the numerical model of a bottom arch of an orbital cavity using a FEM. Based on the data obtained from computer tomography, the model of a healthy orbit was proposed. The results obtained from numerical analysis may serve as a basis for further investigations concerning stresses and deformations in orbital implants, including a direct implant application. Due to its geometrical complexity the whole skull (and especially its facial part) modelling presents a substancial engineering problem. To resolve it one subject the area of interest (the surface or the space) to a segmentation by means of finite number of elements averaging the physical state of the body. To generate the maps of stresses, deformations and displacements prior to calculations, it is indispensable to prepare the geometrical and material data of the model.
The main subject of section 3 is to present a model of a double layered pelvic bone and some phenomena during leg flexion, extension, adduction and abduction, using finite element method (FEM). In the musculoskeletal system, the pelvis is one of the most important bones. It is a support of whole body, transfer external and gravitational loads across sarco-illiac and hip joints (Mrozowski & Awrejcewicz, 2004). According to reference (Wall, 2002), for the population over 65 years, more than 50% of all chronic diseases are joint diseases, most often spine and hip joint. Fracture rate is significantly higher for women over 65 years with ostheoporosis – even 4.6% (Cranney et al, 2007). A properly build model of this bone should bring satisfactory results, minimizing or just eliminating a difficulty of possible harmful medical activities on the examined subject. According to (Mazurkiewicz & Topoliński, 2009), the best method for today is BMD (Bone Mineral Density), painless and noninvasive, however it does not descript the shape of the bones. To make a FEM calculations, a model of a pelvic bone was introduced on the basis of high quality computer tomography (CT) images of an anonymous person pelvic bone. The authors propose some simplification of the model. The aim of the paper was to present the algorithm. To speed up the calculations, it was decided to use simple materials properties. Bone was modeled as homogenous, isotropic, elastic material. However, it is well known that bones have non-homogenous, anisotropic properties. To make model more accurate, this should be taken into account. Model of pelvis took into account a dual structure of the bone, consisting of cortical (compact) – outside part and trabecular (spongy) – inside part bone tissues. Stress in a bone comes from the body weight and from the attached muscles generating loads during movements of lower limb. Numerical computation results are presented in form of the stress and strain maps.
Section 4 deals with oscillations of coupled DNA bases. It is well known now, that the structure of the DNA molecule is not static but dynamic and that
2. Modelling of an orbital cavity using finite element method
A starting point for this part of investigations is to obtain a reliable geometry of the examined object, as close as possible to the real one. Nowadays, the most popular and effective way to perform it is the use of computer tomography (CT) scanning.
In the presented analysis a CT scan of an orbital cavity in the form of a DICOM (Digital Imaging and Communications in Medicine) file was utilized (Fig. 1).
Based on this file, a CAD eye-socket model was created by using a Solid Edge program (Fig. 2). The last step in preparation for calculations was the conversion of the CAD file to the format accepted by Ansys Workbench program.
Next, the grid of finite elements for the model subjected to numerical analysis has been generated. Covering the model’s area by finite element network is based on the relationship between stress and strain (Hooke’s law):
where σ and ε denote stress and strain tensors, respectively, an
Cortical bone material properties used in numerical calculation process was adopted from reference (Furusu et al., 2001), see Table 1. Loading, in the form of a concentrated force, was applied to the external part of a lower arch of zygomatic bone, in the place being statistically the most exposed to impact or pressure (Fig. 3).
The same procedure as in the case of lower orbital arch was applied to the whole orbital cavity. A CAD model prepared in a Solid Edge program was exported to Ansys Workbench program in which geometry of the whole eye-socket was generated (Fig. 5).
Then, according to the aforementioned procedure, the whole model of eye-socket has been meshed using tetrahedron elements (Fig. 6).
As this was the case for a lower arch bottom, a force located in a lower arch of a zygomatic bone has been applied to the model. The force value as well as the bone material properties (Young modulus and Poisson ratio) was the same as in previous case.
The results of numerical analysis were presented in the form of the map of stresses reduced according to the Huber-von Mises hypothesis, based on the strain energy of distortion. In conformity with this hypothesis the reduced stress can be calculated as follows:
where are the normal stresses and are the shear stresses.
For the assumed boundary conditions and the place of force application the distribution of reduced stresses for a lower orbital arch was obtained (Fig. 7). It is readily to notice that the highest stress occur in its paranasal part. In case of an impact this is a potential region of bone fracture due to exceeding of allowable stress values.
The results of the analysis of the whole orbital cavity display a slightly different map of reduced stresses (Fig. 8). The stress concentration occurred here, in the external part of eye-socket. This is a confirmation of characteristics of real isolated fractures of orbital bottom.
The last stage of a whole orbital cavity investigation was to calculate the total deformation. It may be concluded that the greatest deformations occur in the middle part of zygomatic bone (Fig. 9).
3. FEM pelvic bone modelling
The very first step was to prepare a model of a bone. For this purpose a wireframe model was created by using AutoCad. A set of 54 CT scan slices of heavy woman (100 kg, 50 years old), made with 5mm intervals along Z axis were used (see an example in Figure 10). Measuring of the bone layers thickness was made on the basis of CT scans and then compared with literature. Then, consequently slice images were imported to corresponding layers and one after another, closed boundary contours in XY plane, were sketched with straight lines. After that, all lines of adjacent contours were connected with straight lines in order to create triangles. The resultant wireframe was exported to ANSYS and meshed using triangular shell elements for cortical tissue (4 node triangular version of Shell 63 element type with real constant of 3mm for thickness and default element size of 5mm) and tetrahedral shape solid elements for trabecular tissue (10 node Solid 92 element type in tetrahedral shape, default element size was 5mm). In the end it gave 156738 of elements and 181980 nodes.
Because of strain concentration in acetabulum, it was decided to make there more precisely mesh, which gives 39286 surface elements and 194154 volume elements. Material properties of cortical bone and trabelacural bone were taken from (Kutz, 2003) and amounts: Young’s modulus of 20000 MPa and Poisson’s ratio of 0.3 for compact bone, and Young’s modulus of 200 MPa and Poisson’s ratio of 0.4 for spongy bone part. The thickness of the bone layers was taken from manual measurements compared with (Dalstra & Huiskes, 1995) and amounts 3 mm for the compact bone – the rest is the spongy one. The stress-strain state of the real pelvic bone is an effect of loads coming from the muscles attached to the bone through tendons and the weight of the upper part of the body. According to that, following boundary conditions were used (see Table 2 (forces) and Table 3 (movements restrictions)).
It should be noted that the used values were based on maximum forces determined by multiplication of crossectional area of each muscle in its physiological state and a value of maximum stress which can be induced in the muscle σmax=0.5MPa.
In the assumed model we applied the following schemes of supports and movements restrictions (see Table 3):
Values of forces acting on the pelvic bone were taken from literature. Static forces values make calculations faster. During flexion following muscles take a part: rectus femoris, sartorius, iliacus, gracilis, pectineus, during extension: gluteus maximus, semitendinous, semimembranous, biceps femorisgluteus medius, during adduction: semitendinous, semimebranous, adductor magnus, adductor longus, adductor brevis, during abduction: tensor fasciale lateae, gluteus medius sartorius. An assumption was made that all the movements are made while standing on one leg only. From (Phillips, 2007) and (Huston, 2007) mass percentage of lower limb was estimated, equal to 19%. The resultant weight of 800 N coming from the upper body and single lower limb, was directed downwards along the Z axis and applied on nodes related with area of connection between pelvic and sacrum bones. For every motion a separate model was prepared with its boundary conditions. An example for flexion and extension can be seen in Figure 10.
3.1. Results and discussion
For all four load variants related to basic motions of flexion, extension, adduction and abduction the results were analyzed by plotting the three principal stresses, the Huber von Mises equivalent stress, total displacement and the Huber von Mises equivalent strain. Plot in Figure 11 shows the Huber von Mises equivalent strain for flexion. The biggest value is placed in the middle part of acetabulum and amounts about 0.005.
From 1st principal stresses it can be observed that the biggest tensile stresses occur during extension at acetabulum and are equal to 32.6 MPa (Fig. 12). Greatest compressive stress of -31 MPa comes from the 3rd principal stress and is found also for extension on the surface of acetabulum (Fig. 13).
The highest value of the Huber von Mises equivalent stress, as it was easy to expect, is also characteristic for acetabulum surface during extension movement and is equal to 28.6 MPa (Fig. 14). Here it can be seen from the crossectional view that the highest stresses are placed at the external parts of the cortical layer, where stresses for trabecular layer are much lower. Similar effect can be observed in models of other authors, even for full solid models (see (John, 2004) to compare).
Regarding the total displacements the biggest and equal to 0.28 mm are observed for extension (Fig. 15), located at ischial tuberosity and for abduction equal to 0.23 mm at front part of iliac crest (Fig. 19).
For all of the presented motion variants the greatest Huber von Mises strains occur in acetabulum. Similar conclusions can be found in the works (Kutz, 2003) and (Phillips, 2007). The highest strain is found for extension and amounts to 0.009 (Fig. 16) compared to adduction 0.006 (Fig. 17) and abduction 0.004 (Fig. 20).
Analysing the stress distribution for rest of the movements the maximum principal and equivalent stresses occur on acetabulum, except adduction for which greatest stresses are placed on the inside surface of superior ramus of ischium, where equivalent stress is equal to 14.8 MPa (Fig. 18).
The smallest values of displacements are found on constraint surface of acetabulum for all motion variants.
The above observations prove that the acetabulum of pelvic bone is a main subject of pathology of the whole pelvis girdle. As a place of biggest stresses and strains concentration, acetabulum is exposed to a high risk of degeneration.
Earlier works on modeling pelvic bone, as an example (Dalstra & Huiskes, 1995) had not enough accurately mesh (inter alia, because of the possibility of computing) or models of the pelvic bone were simplified, for example authors were used only surface model (only trabercular bone), without the spongy one (Phillips, 2007; John, 2004). That is the reason why there is still a place for improvement of model and loads (Phillips, 2007) or models. This simulation was made for static situation – one leg standing. It is well known that the hip joint during normal walking is around 300% body weight (Bergmann et al, 2001) (we use 100% of body weight – mass of the leg + muscle forces). FEM method can give consistent results with medical observation even for more complicated situation as for example sideways fall (Majumder et al, 2007).
4. Modelling of DNA bases
The analogy permits to reduce the problem of oscillations of complementary DNA bases forming the base pairs (AT or GC) to the mechanical problem of oscillations of two non-identical coupled pendulums shown schematically in Fig. 21. In reference (Yakushevich et al, 2009), linearized (particular) version of the problem has been investigated. Here we consider the nonlinear (general) version.
The function of Lagrange of the model imitating oscillations of two non-identical coupled pendulums shown in Fig. 21 has the form
where , and are the lengths of the first and of the second pendulums, respectively;
4.1. Solutions of model equations and trajectories in the configuration space
We solved numerically equations (6)-(7) with the help of the MAPLE. To construct corresponding graphs, we used the values of the coefficients
Graphs of the solutions of equations (6)–(7) obtained for AT and for GC base pairs and for three cases: (a) ideal case when effects of dissipation and external fields are neglected; (b) non-ideal case when effects of dissipation are taken into account but effects of external fields are neglected; (c) non-ideal case when both effects of dissipation and effects of external fields are taken into account, are presented in Fig. 22. The trajectories of the considered dynamical system are presented in Fig. 23.
When constructing the graphs of the solutions and trajectories, we used initial conditions:
which are closed to the point (0.2753, 0.3347, 0.0, 0.0) in the phase space, and we shall show in the next section that this point is stable in the case of AT base pair.
Summarizing the results presented in Fig. 22-23, we can state that effects of dissipation lead to decreasing the amplitudes of base oscillations. Effects of external field lead to displacement of the equilibrium positions, displacement in the case of GC base being less than displacement in the case of AT base pair. So, we can state that polynucleotide chains saturated by AT base pairs are more sensitive to the action of external field than the chains saturated by GC base pairs.
4.2. Analysis of stability and diagram of bifurcation
To analyze the stability of the system (6)-(7), we firstly transformed the system of two differential equations of the second order (6)-(7) to the system of four differential equations of the first order
Since system of equations is nonlinear, we cannot find solution analytically but we can do it numerically. We found that for each base pair there are four equilibrium points. Only one of them being stable, as all eigenvalues of the Jacobian of the right side of (8) – (11) have negative real parts. It is P4AT = (0.2753, 0.3347, 0.0, 0.0) point in the case of AT base pair and P4GC = (0.1658, 0.1658, 0.0, 0.0) point in the case of GC base pair.
Below we present the diagrams of bifurcation of the analyzed system. In general, diagrams of that type show possible equilibrium points or periodic orbits of a system as a function of one of the parameters of the system. For every value of that parameter we can choose many initial conditions and find the trajectories. If trajectory is attracted to the stable equilibrium, we obtain only one point in the diagram, if trajectory is approaching to periodic (or quasi-periodic) orbit we obtain two (or more) points in the diagram. If each trajectory is attracted by different point, we obtain a lot of points in the diagram which mean a chaotic behavior of the system. Bifurcation occurs when the change of the value of parameter causes the changes in the behavior of the analyzed system.
In Fig. 24 (a) - (b) we present the bifurcation diagrams as a function of k01 parameter, describing the external force. Because we assumed that the external force is constant, the increase of value of k01 parameter does not change the dynamical properties of the system (in both, AT and GC, cases) and we do not observe the bifurcation. It is worth to notice, that bigger value of k01 corresponds to the bigger displacement of the bases from the (0.0,0.0) equilibrium position
In Fig. 24 (c) - (d) we present the bifurcation diagrams as a function of b01 parameter, describing the parameter of dissipation. For very small values of b
In fact, trajectories are very slowly approaching to the stable point, and in the bifurcation diagrams we observe transient state of the system for small values of b
One of the most dangerous consequences of the road accidents is the facial skull injury. Especially, the eye-socket damage can cause the serious vision problems and brain dysfunctions. Hence we have presented the results of stress and strain analysis of an orbital cavity. Thanks to the numerical analysis in Ansys Workbench the concentration of the reduced stresses in eye-socket had been confirmed. In case of isolated fractures of an orbital bottom, the lower edge of an orbital cavity due to its elasticity does not break but only bends transmitting the force to a brittle bottom of an orbital which is prone to fracture (Smolarczyk-Wanyura, 2005). Further work is required in case of numerical model optimization but also in the fact that simulation must include various loading cases. The study may provide an assistance for surgeons performing bony face operations.
In the case of a pelvic bone our results show, that most of the highest stresses and strains occurred in acetabulum. The comparison with medical data cases (Schmalzried et al, 2004), indicates that degeneration and wear of acetabular regions is the main case of pathology in the area of pelvic girdle. It is also consistent with the results presented in (Będziński, 1997). Results are comparable also with papers (Dalstra & Huiskes, 1995) and (Phillips, 2007), i.e. the same magnitude of results is reached. This method can be used for examination of human pelvic bone and can help, after some additional measurements regarding the bone density estimation, to predict when the failure (and if) could happen. It is also possible to test prototypes of pelvic bone implants.The presented method of creating a 3D model from images produced with Computer Tomography, although based on many simplifications, can be used for many other analyses, where the examined subject has an irregular shape. The presented model is not an ideal one. Each pelvic bone is different. The shape of it depends on sex and age. The presented algorithm serves as a proposal of a method of pelvic bone modeling. In future, if enough cases will be modeled, it will be possible to found a database of ready to use models. Minor changes will be needed for a data fitting. There is still a lot of space for model improvements like building a much more detailed 3D geometry model, creation of much finer mesh, considering anisotropic properties for material or modeling phenomena during walking, running, jumping (like in (Schmalzried et al, 2004)) or even skiing, when because of higher forces value, hip fracture is more likely. Creation of more accurate model and denser mesh will be possible, provided that more CT scans are available – depending on computational resources.This will allow to avoid mistakes such that bone contacts only in the fascia lunata and do not in the center where the ligamentum capitis femoris placed. Due to technical limitations adopted element are of size 5mm and constant thickness amounts 3mm. Authors of paper (John & Wysota, 2010) propose a method based on CT images analysis of tissue density and determination of Young modulus. They took into account the dual structure of the bones.
In section 4 we have studied rotational nonlinear oscillations of two coupled DNA bases that form a base pair: adenine-thymine (AT) or guanine-cytosine (GC). We have presented the graphs of solutions and trajectories in configuration space. We have also found stable equilibrium points for each base pair and constructed diagrams of bifurcations. This work part can be considered as a base for more complicated model of DNA, including not only one base pair, but a fragment of DNA consisting of three base pairs.