© The Author(s), Published by EDP Sciences, 2024

1 Introduction

The problem of determining bone health in vivo has been the subject of intensive research for decades. The majority of this work has focused on assessing the impact of osteoporosis (OP). OP is a bone condition which describes the weakening of bone to the point of fracturing from low trauma falls i.e. from standing height or lifting a mildly heavy item. OP is a bone condition which is defined as “a disease characterized by low bone mass and microarchitectural deterioration of bone tissue, leading to enhanced bone fragility and a consequent increase in fracture risk” [1]. The condition costs healthcare services around the world several billion dollars each year owing to treatment mainly on hip fractures and the co-morbidities [2–4]. The World Health Organisation (WHO) guidance on diagnosing OP uses Dual X-Ray Absorptiometry (DXA/DEXA) to measure the Bone Mineral Density (BMD) of an area of interest, such as the hip or spine [5]. This BMD is then compared to a reference value derived from a population of White women under the age of 30 [6]. The severity of OP is determined by using standard deviations from the reference BMD. A value of less than −2.5 SD away is considered the diagnostic threshold for OP. Values between −1.5 and −2.5 SD are given the diagnosis of Osteopenia (OA) which is considered a less severe form of OP. But it should be stated that the majority of fragility fractures occur to the population within this bracket [7], thereby highlighting the need for early diagnosis.

Bone tissue is found in two forms: cortical bone and trabecular bone (also sometimes referred to as compact and cancellous bone respectfully). The dense and stiff compact bone is found on the external surface of all bones but is particularly concentrated along the midshaft of the limb bones. This provides most of the long bone’s twisting and bending stiffness, while the porous and elastic trabecular bone provides some bending stiffness while reducing mass [8]. The trabecular structure also is a source of strength through the controlled production and repair of microcracks to absorb energy input [9]. There is also collagen and bone marrow which fill inside the pores of the trabecular bone and in the midshaft which contribute various structural benefits but are not the focus of this investigation [8]. OP affects the bone reproduction process by increasing bone reabsorption and reducing bone material deposits, creating a net loss of bone mineral [10]. This appears in trabecular bone as fewer cross-strut connections and in thinner structures, while compact bone becomes more porous and the interior surface layer of bone erodes [10, 11]. This therefore results in a fall in BMD and total mass, but also a fall in stiffness from the weakened structure and loss of material. It has yet to be determined if the fall in mass or the fall in stiffness is more prominent as OP progresses, but the research tends towards the belief that it is mass which is directly reduced and stiffness is affected less though compensation in remodelling [9, 10, 12, 13].

There have been previous suggested alternative methods to determine the health of bone. Quantitative Ultrasound (QUS) has been suggested as a suitable means of recording the BMD and a direct measurement of the Young’s modulus through the attenuation of ultrasonic waves through the trabecular microstructure [14, 15]. Yet adoption is limited owing to clinical barriers plus a limited improvement over the DXA method [16]. Mechanical Response Tissue Analysis (MRTA) and its successor Cortical Bone Mechanics Technology^TM (CBMT) [17, 18] uses two-point bending which is a non-destructive method to measure bone stiffness [19]. It has been well researched and found some good results, but again has faced barriers to adoption owing to the competition with DXA and poor repeatability, among other issues [20]. The final method is to use vibro-acoustic excitation to excite the bone’s modal frequencies. This method uses an impact hammer and an accelerometer (or other transducer combination) placed on one of the long bones, typically the tibia. The modal frequencies can then be analysed for their relationship to the stiffness of the bone. This method has been researched on and off since the 1970s [21–25] but it too has faced challenges both from the practical and theoretical aspects, particularly in the use of a metric which is dependent on both stiffness and density which are confounding factors in the determination of the resonance frequencies [26].

There have been several attempts to calculate the eigenmodes and resonance frequencies of long bones. Much of the work in modelling bone has focused on simplifying the vibration behaviour of bone using analytical and numeric mathematical models. The first attempts used simple slender continuous models which attempted to simplify the structure, material, and shape properties; other works go further and include the use of Timoshenko models or FEM which account for the anisotropic properties of bones including the changing shape moment of inertia [27]. Several authors have claimed to achieve an acceptable model that describes the long bone’s vibration response [25, 28, 29]. But there have been doubts raised such as using the fundamental frequency change to determine the bone health [22, 24, 30]. Others have questioned the removal of physical details from the models such as rotational inertia and changing moment of inertia which appear to have an important impact on the modal response [27, 28, 31]. The boundary conditions have the largest impact on the vibration response, but in vivo bone is surrounded by several physical constrictions and influences from muscle tissue, the fibula, and tendon structures which makes modelling extremely difficult [32]. Therefore, bones are often modelled in vitro with free-free conditions, a convention we continue in this work.

Focusing on the work involving FE methods, several authors have attempted to include more realistic shape and material properties [27, 32–37]. Thomsen observed that a simpler model with only a few crucial features included such as shear deformations is sufficiently close to experimental results, suggesting that the tibia can be modelled more uniformly despite the geometry and material distribution [27]. Hobatho and co-workers investigated the use of FEM to evaluate how the response to vibration is related to bone strength, showing good agreement between numeric and experimental results but expressed limitations with using resonance frequencies as a measure of bone strength [36, 37]. Recent work in this field has continued this research to include accurate 3D scans with an interesting study on bone shape variation by Campoli and co-authors [33]. But there is a gap in the literature for a more thorough investigation into the ways that the cross-sectional shape and other features can be isolated, simplified, and examined more closely.

While the analytical continuous models are not the focus of this work, there are some useful insights relevant to the studies which are mentioned here. The most relevant equations include the modal frequencies of a free-free Timoshenko model [38]:

$${ F}_n= \frac{{{\gamma }_n}^2}{2\pi L^2}\sqrt{\frac{\mathrm{EI}}{\rho A}}$$(1)

where L is the length of the beam, γ_n are the roots of the solution to the modal equation [39], EI is the stiffness term (bending rigidity), E is the Young’s/elastic modulus and I is the mass moment of inertia, also known as the 2nd moment of inertia. ρ is the density of the material and A the cross-sectional area. The stiffness is itself influenced by changes to elastic modulus and the 2nd moments of inertia I, defined as:

$$I_y={\iint }_R^{ }{y}^2\mathrm{d}A$$(2)

$$I_z={\iint }_R^{ }{z}^2\mathrm{d}A$$(3)

where R is the region representing the cross-section. The area moment of inertia is typically different in the y and z directions and depends on the choice of reference point. In this paper, the reference point is the centroid of the cross-sectional shape R. The modal frequencies are therefore directly influenced by changes in the stiffness and the density of the material. It should be mentioned that other wave types such as torsional and longitudinal waves been modelled in other literature and will be mentioned here [28]; trabecular bone can transmit Biot’s fast and slow waves though its solid structure and fluid pores but these waves are not the main focus of our work here [15].

This paper investigates the use of an FE model of a tibia and femur bone to generate modal frequency data and is an extension of previous reported work on the subject [40]. The key contribution if this paper is the investigation to quantify how geometrical features of the bone models affect the modal frequencies. The full bone shape of a tibia and femur is analysed first, then the thicker ends of the bone are removed parametrically to assess their effect on the modal frequencies. Then the investigation will look at experimenting with several cross-sectional shapes, suggested from the literature and by observation of the bone cross-sections, to quantify how this simplification will change the modal frequencies. Two features of the bone shape are investigated: twist along the length of the model plus scaling the ends larger than the midshaft. A discussion and analysis of the results follow, then some conclusions are drawn from the results and discussion.

2 FEM experiments – “In Silico”

This section describes the use of FEM to explore the effect of separate features of bone shape on the modal frequencies. FEM simulations of a bone shape with typical bone material properties were deemed as adequate for the purposes of the investigation. The change in modal frequencies, either in absolute values in Hertz or relatively as a percentage, is the main metric to judge the extent of these simplifications.

2.1 COMSOL model setup

Free and publicly available 3D scans of a tibia and a femur bone were imported into COMSOL 6.1 as .stl files [41, 42]. The Solid Mechanics physics module was used with an eigenfrequency study. A single isotropic material was used to represent both the compact and cancellous bone. The 3D model did not contain any internal layering to the geometry; the distribution of cortical and trabecular bone in different parts of the tibia and femur requires a more complex structure which was not possible to generate from this geometry [27]. The model material properties were taken from selected literature to allow comparison with previous experimental and model research [34]. The details of the parameters used in the studies and other data in the model are given in Table 1.

2.1.1 Full bone

The eigenmodes and modal frequencies of the full bone shapes were studied to give a reference point for later comparisons with cross-sectional shapes and features. The eigenfrequencies are listed in Table 2.

As seen from previous works, the modal response contains two series of bending modes in the two principal directions, plus torsional and longitudinal modes in between [31, 34]. The categorisation of the modes was made from visual inspection of the mode shapes generated by COMSOL The simplification of the material, differences in the bone dimensions, and differences in methods for the categorisation of modes are likely the contributing factors to the mismatch of the modal frequencies. The anisotropic material properties will make bending difficult in certain directions but not others, thereby making bending modes in one direction stiffer and therefore greater in frequency. But their pattern is similar to that reported in literature [31, 33, 34]. Further analysis of the whole bone’s modal response has been extensively covered previously both for tibiae and femurs [23, 35, 43]. A basic sensitivity analysis was carried out using FEM to understand the influence of the fundamental properties on the modal response. The elastic modulus, the density, the length and the “radius” i.e. the lateral dimensions were selected as variables which are varied by ±10% each individually. Cross-modulation effects are to be not investigated here. The results are given in Table 3.

The dimensions of the bone were the major factors in influencing the modal frequencies, with a 10% change in length producing a 20% change in the modal frequencies. An increase in the radial component of the bone was directly proportional to the increase in the resonance frequencies of the bending modes. Increasing the density or Young’s modulus produced a 1:1/2 ratio change of the modal frequencies but in opposite directions. So, a 1% increase in density resulted in a 0.5% decrease in the modal frequencies, while a 1% increase in elastic modulus produced 0.5% increase in resonance frequency. With reference to equation (1), the length has a strong influence on the modal frequencies owing to the 1/L² term. The limited influence of the radius is owing to the confounding influence of the inertia term in equations (2)–(3) and the cross-sectional area – this will become clearer with the simpler cross-sections but the principle holds here. This limited sensitivity study is useful for understanding which parameters require the strongest control and highest accuracy for measurement. It would indicate that, in modelling terms, the Young’s modulus and density do not have to be highly accurate (within 20% of the real value) before the change the modal frequencies by 10% is achieved. Length on the other hand needs to be very accurately controlled to make proper comparisons between bones.

2.1.2 Cut bone

In “thin beam theory” using the Euler Bernoulli model, it is assumed that the depth of the beam is significantly smaller than the length of the beam, expressed as λ << d where λ is the wavelength of the vibration and d is the depth of the beam. As the length of the beam bounds the wavelengths of the resonance modes, they are given as equivalent terms. This is to avoid considering cross-dimensional modes and the effects of rotational inertia and shear deformation which is considered in “thick beam theory” of Timoshenko–Ehrenfest beams [39]. In long bones, the distal and proximal ends are typically thicker to accommodate the joints and provide attachments to muscles and tendons. It has not yet been found in the literature to date to what extent these thicker ends influence the vibration response of bones. Ex vitro studies have removed soft tissue and other influences but with virtual models we can also dissect the bone in many different ways and locations without requiring several real bones [22, 26, 44]. It also has the added benefit of being able to artificially equalise the length and other parameters to control their influence on the modal frequencies. This allows us to just focus on the influence of shape on the modal response.

The cross-section of the two bones changes considerably over the whole length of the bone, but over smaller lengths the changes are more gradual. It therefore may be difficult to specify a point at which the midshaft or thicker ends begin. To investigate this, two cut planes are created at the distal and proximal ends of the bone model. The distances are parameterised and measured from the 0 point of the COMSOL model at the distal end of the bone. This is demonstrated in Figure 1, showing the two thicker ends of the bone highlighted in blue.

Figure 1 Tibia model in COMSOL with the “thick” ends highlighted in blue. The distances at which designate these ends are captioned Cut L1 for the distal [ankle end, left] and Cut L2 for the proximal [knee end, right].

Cutting into the bone shape at regular intervals, the changing cross-section can be observed (Figs. 2 and 3). Starting at the proximal end, the cross-section starts asymmetric and large, compared to its midshaft where it rapidly reduces in size and changes shape to become more rounded and closer to being symmetrical on one axis. Towards the distal end it begins to become more symmetric before expanding again into a larger, squarer shape. The femur, likewise, has the unique shape of the femoral head and the trochanter, then becoming more circular and shifting its centroid along the midshaft. It then steadily increases in area towards the distal end where the cross-section is symmetric but much larger and with the medial and lateral condyles for the knee joint.

Figure 2 Cross-sections at selected points along the tibia, starting top left (a) at the distal end, tow points along the midshaft (b) and (c) then the proximal cross-section (d).

Figure 3 Cross-sections at selected points along the femur, starting top left (a) at the knee joint, then moving along the midshaft (b) and (c) then the cut section of the femoral head (d).

The tibia bone and femur bone were dissected with these cut planes and the lengths swept with the values in Table 4. The bone was then scaled in the X direction to preserve the overall length. Figures 4 and 5 show the first and second pairs of bending modes respectively of the femur model. The change in Cut L1 is represented as the separate line colours and Cut L2 is plotted along the x-axis.

Figure 4 Resonance frequencies of femur varying with Cut L1 length (see legend) and Cut L2 length (along x-axis). Bending mode 1 and 2 captioned.

Figure 5 Resonance frequencies of femur varying with Cut L1 length (see legend) and Cut L2 length (along x-axis). Bending mode 3 and 4 captioned.

It can be observed that cutting the proximal end of the bone (where the trochanter is) has a major impact on the modal frequencies, while the distal end has a more subtle effect. Whereas the decreasing length of Cut L2 (i.e. including more of the proximal end) leads to a rapid decrease of the modal frequencies, some non-linearity in this dependence is observed for F2 bending mode. This is obvious given the unusual shape and the much greater size of the femoral head and the trochanter compared to the midshaft. The effect of the proximal end becomes apparent by the change in the gradient around 320 mm from the distal end for the first 2 bending modes. The distal end has an effect but only if less than 35 mm cut – the change in the frequencies comes much smaller for all but the first bending mode. A summary of the relative changes of the removal of the distal and proximal ends is given in Table 5. The results for the tibia are plotted in Figures 6a–6d showing each bending mode frequency with changing Cut L1 and Cut L2 as described previously.

Figure 6 (a) Resonance frequency of bending mode 1 of the tibia varying with Cut L1 length (see legend) and Cut L2 length (along x-axis). (b) Resonance frequency of bending mode 2 of the tibia varying with Cut L1 length and Cut L2 length. (c) Resonance frequency of bending mode 3 of the tibia varying with Cut L1 length and Cut L2 length. (d) Resonance frequency of bending mode 4 of the tibia varying with Cut L1 length and Cut L2 length.

For the tibia, the results are more unintuitive: removing more of the distal end again has a predictable increase in frequencies, while the proximal end appears to decrease the bending modes in the anterior direction (F1 and F2) but increases the frequencies in the lateral direction (F2 and F4). Removing over 25 mm to 30 mm from the distal end reduces the change in modal frequency, indicating that beyond this region the midshaft is fairly uniform. A summary of results is given in Table 6 by calculating the relative changes to the modal frequencies using the extents of the ranges chosen.

2.2 Shapes

The cross section area and shape of long bones change over their length. The area and shape are important in the determination of modal frequencies both directly and through the moments of inertia. It is seen in Figures 2 and 3 that the midshafts of both the femur and tibia are rounded but asymmetric cross sections, with more circular and triangular features respectively. The proximal and distal ends of the long bones are much larger in area and less regular in shape, and not possible to simplify by a more basic geometry. Several authors have made simplifications to the cross-sectional shape to make calculation of the modal frequencies easier [25, 28, 45]. These often reduce the shape to a circular cylinder, but other shapes can be used which describe aspects of the cross sections found in the dissections. This section explores the use of some selected cross-sections to simplify the modelling of bone. These shapes are inspired from their use in the literature and from observation of the cross-sections of the bone models and from previous authors [27, 34].

2.2.1 Circle

The most common simplification is to assume a circular cross section. This appears most often in analytical solutions as it simplifies the mathematics of the moment of inertia and calculating the cross section. Using a circular cross section gives us only control of the radius as a shape parameter. The radius is swept over a range of values to investigate if using the thick or midshaft dimensions are required to accurately match the modal frequencies of the real bone. The results are given in Table 7.

The cross section has two axes of symmetry which results in repeated modes along perpendicular directions. This therefore cannot model the separate series of modes in two different directions as seen with the real bones. This is not mentioned by the authors who have used cylindrical or circular beams in their models before, though most have focused on just matching the fundamental bending mode which reduces the requirements for the model’s complexity [28, 45]. As found with the real bone sensitivity tests, the increase in radius produces an increase in modal frequencies except the torsional resonances which remain mostly static. As the circular cross-section has symmetrical moments of inertia where I = πr⁴/4 while the cross-sectional area is A = πr², these terms result in an r² term within $\sqrt{\mathrm{EI}/\rho A}$ which explains the direct proportional relationship with the bending modes.

2.2.2 Ellipse

To capture the difference in bending directions from moments of inertia, the circle can be morphed into an ellipse. This is characterised by two dimensions, A and B, producing an A:B ratio which can be parameterised. By changing the A:B ratio, the effect of change in moments of inertia on the eigenfrequencies can be investigated. The ratio is varied from a range of 1.0–2.0 in 0.1 steps, and the results are given in Table 8.

The frequencies of bending modes across the two principal axes begin to show the split as seen in the real bone, with the gap between the “A” and “B” modes increasing with increasing ratio. By normalising the modal frequencies to their fundamental, the modal pattern becomes apparent as the ratio between the A and B directions increases. Yet the absolute values for the modal frequencies do not match the tibia and the femur frequencies – adjusting the ratio to better match one of the resonance frequencies worsens the error at other modes. A:B ratios of 1.8 and 1.9 gives an error for ±2.7% respectively from the first modal frequency of the tibia, but is over 150% off from the second bending mode frequency and is −12.5% and −17.1% from the third bending mode. Taking a ratio of 1.6 achieves a closer match to the third bending mode with an error of −1.9% but the error for the fundamental bending modes in the A and B directions are 15.3% and 130.5%. As it was with the circle radius sweep, the shifting of the modal frequencies re-orders the appearance of the bending and torsional modes, which makes identifying modes from the eigenfrequencies alone very difficult without a knowledge of the modes those frequencies are connected to. Comparing to the real bone, though the two series of bending modes have now been represented but the absolute values of the frequencies are off compared to the tibia or femur.

2.2.3 Triangle

Focusing on the tibia cross-sections in Figure 2, it may be observed that the midshaft has three vertices and three sides, two with flat surfaces. This can be simplified into a saline triangle which can be parameterised by the internal angle at one of the vertices. Very few authors have attempted to use this interpretation; only Collier and authors analysed the modal response of the tibia with a triangular prism cross-section and claimed to achieve a fundamental frequency of 110 Hz in vivo [31]. A triangle with the same cross-sectional area as the ellipse and circle models previously is generated and one of the internal angles is varied from 60° to 30° to produce and equilateral, saline, and right-angled triangle [40]. The eigenfrequencies using these shapes are recorded in Table 9.

The equilateral triangle cross-section has a different modal frequency series compared to the circle or ellipse, but because of its symmetry the modes are repeated just like for the circular cross-section. By changing the internal angle to create a saline triangle, the change in moment of area is achieved and the modal series splits between the two principal axes of bending. The more acute the angle, generally the greater the difference between the first and second bending frequency series.

2.3 Twist

A feature of the long bones shape is the way the principal axes change along their length. The absolute position of the centroid also changes along the length of the bone, but for brevity is not investigated in this paper. This twist has been mentioned in previous papers by Thomsen as well as Hight, Piziali, and Nagel respectively but they reached conflicting conclusions on its effect on the modal frequencies [27, 32]. We have investigated this effect to confirm which is the case in our simpler models.

COMSOL implements twist using the extrude geometry function. This transforms the cross-section along a length and then offsets the vertices by a twist angle. This has the side-effect of reducing the cross-sectional area in the middle of the twisted extrusion. There has been no attempt to resolve this issue but may possibly be a compounding factor in the results. The twist angle was varied from 0° to 45° and the results for the circle and 50° triangle are in Table 10.

The effect of twist on the circular modes shows a 1.5–2.5% decrease in the modal frequencies for every 20° twist. But this depends on the type of mode: torsional modes of course are more affected with much greater rates of change compared to the bending modes.

This would indicate that Thomsen’s conclusion that “twist of the tibial shaft is almost insignificant” is correct [27]. Even with a twist of 45° the modal frequencies are decreased by less than 10% for the circular cross section, but it may be a subtle effect when combined with another shape factor.

2.4 Scale

As evidenced by cutting the proximal and distal ends of the bone, the larger cross-sectional area and change of shape contribute to the bone’s vibrational behaviour. But the question remains how we best simplify these ends in a parametric way. What is proposed is a way of simulating the changing cross-section area over a limited part of the length at the ends of the model: the cross-sections at the end of the model are scaled by an independent factor and are extruded over a defined length but while reducing the cross-sectional area. This is done separately on the other end of the model. An example of this is given in Figure 7. The scales and lengths of both ends were independently swept as per the ranges specified in Table 11.

Figure 7 A COMSOL model geometry showing the two emulated distal and proximal ends of the bone. The scale of the ends is adjusted using SD and SP [Blue Face], and the length of these ends adjusted with LD and LP [Red Lines] respectively.

From the cut bone analysis, the proximal and distal ends are about 332–745% greater in cross-sectional area than the midshaft, and removing the ends to 50 mm either end reduces the frequencies by 10–15%. The ranges set above attempt to avoid excessive computational overhead when sweeping over the results, but still give an indication of the effect over a smaller scale. The overall bone length is kept consistent but overall mass is left uncontrolled. This study uses just the circular cross-section for brevity. Abridged results are given in Table 12.

The results show that modal frequencies are reduced by 5.0% when the scale at the proximal end increases from 1 to 2.6 times the midshaft area. When the length of the proximal end is increased from 0.01 m to 0.06 m and kept at 2.6 times the midshaft area, the frequencies fall by 14.9% compared to the uniform model. When both proximal and distal ends of the bone are 2.6 times larger than the midshaft and both 0.06 m in length, the bending modes reduce substantially to −27.7% from the original model without these end scales.

3 Discussion

The tibia and femur bones have been the most studied of long bones in vitro, likely owing to their relative ease of access in vivo and their relevance in hip fractures respectively. As the bone shape includes a long slender midshaft, there is an opportunity to simplify the complex shape into a simpler model to identify the key parameter changes resulting from OP. These physical parameter changes can be identified by the change in modal frequencies, in particular the bending stiffness EI and the density ρ. As our primary interest is in the resonance frequencies themselves, the magnitude of these resonances are not considered in this work, given the number of other factors involved which are unrelated to the resonances such as input force for instance. The modal response of the tibia shows a wider frequency gap in the pairs of similar order of bending modes compared with those of the femur, likely from the more distinctive difference in moments of inertia in the two bending directions. Removing the ends of the bones changes the eigenfrequencies considerably, especially for the femur with the femoral head which can produce an 10.6% change to the fundamental bending mode if the bone is cut at 250 mm from the distal end. Cutting the distal end itself produces a smaller change of at most 2.5% when cut at 30 mm. Owing to COMSOL geometry errors, further cuts were not possible which might have shown a similar large influence on the resonances. The effect is more subtle for the tibia bone where the decrease in modal frequencies from the inclusion of the ends is more linear: the furthest extremes of the cuts to the proximal and the distal ends produced a frequency change to the first bending mode of 3–8% respectively. The second bending mode increases in frequency as the proximal end is cut while removing the distal end changed the frequencies by 15%. But returning to the sensitivity analysis and the results from the circular cross-section radius sweep gives a possible explanation to the increasing modal frequency with increased proximal end: the increasing radius of the cross-section increases the modal frequencies, and from the ellipse results the increase is in the direction of the smallest dimension, of which the F2 and F4 are bending in the smaller dimension of the tibia. This explanation is less convincing for the fall of F1 and F3 as the radius increase affects these modes too but with a negative direction. Following equations (1) and (2), the increase in one direction of the cross-section may have a larger effect on the moment of inertia than on the cross-sectional area overall, thereby increasing the modal frequencies in one direction but reducing in the other.

The cross-sectional shapes chosen in this paper were not only inspired both by the physical shape of the bone but also what appeared in previous research. Comparing the results of models using these shapes and the modal frequencies from the real bone, it is clear that there are other factors involved in the modal frequencies which should be considered. The circular and equilateral triangular shapes had repeated eigenfrequencies which is not apparent in real bone, thereby making them inappropriate as simplifications of the bone cross-section. The ellipse and triangle cross-sections had this property, due to the differences in moments of inertia of the two bending directions. Adjusting the A:B ratio showed an asymmetrical change in the frequencies of the repeated bending modes: the modes in the A direction changed less in Hertz than the bending mode in the B direction. This can therefore be a parameter to investigate further with the real bone cross-sections to see if this same relationship holds true. A similar but lesser effect was found when adjusting the internal angle of the triangular cross-section, indicating a similar influence on the moments of inertia.

The longest dimensions of the tibia’s proximal and distal cross-sections are, respectively, over 2.5 times and 1.5 times the bone’s midshaft diameter; removing these thicker ends had a considerable effect to the modal frequencies. When this effect was parameterised in a model with larger ends, the effect was similar on a simpler cross-section – the change can be as high as almost 28% when both ends are 2.6 times the midshaft area and are each 60 mm long. This suggests that this shape factor should be considered in future models as it has a non-trival influence on the modal frequencies both when excluded (as in the cut bone experiment) and included (the parameterised model).

The work covered here extends on the author’s previous reported results on using FEM for bone modelling [40]. The suggestion of a more uniform bone model is in line with the comments made by previous researchers on the subject [27, 28]. Such as with DXA, there might be a possibility in the future to generate a database of reference values to compare with a person’s unique impulse response, perhaps with some individualisation of the reference to make the comparison more relevant. This can be generated from using accurate but fast models such as these FE models to allow for many variations of dimensions without needing to record real data for every instance.

4 Conclusions

This paper investigated the use of the shape parameters of bone on the eigenfrequencies using FE models. The key research question is to examine how the bone shape can be simplified for faster simulation and identifying the key shape parameters which influence the eigenfrequencies. Of these shapes, both the circular and the equilateral triangle have repeated modes which are not apparent in the real bones as these shapes have axial symmetry in each of its bending directions. Only the ellipse and the non-equilateral triangle had the necessary differences in moments of inertia which allow for different modal frequencies for different bending directions. Changing the A:B ratio and the internal angle allowed for further fine tuning of the modal frequencies to better match the real bone eigenfrequencies, but the other factors of scale and twist will be needed to be incorporated together to achieve better closeness.

One of the goals of this research was to join the analytical modelling and FE techniques to allow for the benefits of both in the analysis of bone vibration. The key application of this can be in generating datasets for sensitivity studies into shape without needing to calculate for entire bones. Some of these features would require extensive modification of the analytical models to be accounted for, while FE methods can calculate the effects of these features much quicker and with fine variation possible. In the case of the saline triangle, this approach can inform further analytical development using FE methods to experiment with cross-sectional shapes in the first instance.

There is much more work to investigate these features further, especially focusing on the material properties and internal structures if these have a further effect on the eigenfrequencies. The focus of this work has been on the eigenfrequencies following from the decades of focus on this metric; other vibro-acoustic features in both the time and frequency domain should be investigated for how the shape further affects these features.

Conflicts of interest

The authors declare they have no conflicts of interest to report concerning any aspect of the work mentioned in this paper. The authors have full control over the primary data and welcome any review of the data if requested. The authors received no funding to complete this research.

Data availability statement

Data are available on request – please contact the corresponding author.

References

All Tables

All Figures

	Figure 1 Tibia model in COMSOL with the “thick” ends highlighted in blue. The distances at which designate these ends are captioned Cut L1 for the distal [ankle end, left] and Cut L2 for the proximal [knee end, right].
In the text

	Figure 2 Cross-sections at selected points along the tibia, starting top left (a) at the distal end, tow points along the midshaft (b) and (c) then the proximal cross-section (d).
In the text

	Figure 3 Cross-sections at selected points along the femur, starting top left (a) at the knee joint, then moving along the midshaft (b) and (c) then the cut section of the femoral head (d).
In the text

	Figure 4 Resonance frequencies of femur varying with Cut L1 length (see legend) and Cut L2 length (along x-axis). Bending mode 1 and 2 captioned.
In the text

	Figure 5 Resonance frequencies of femur varying with Cut L1 length (see legend) and Cut L2 length (along x-axis). Bending mode 3 and 4 captioned.
In the text

Figure 6 (a) Resonance frequency of bending mode 1 of the tibia varying with Cut L1 length (see legend) and Cut L2 length (along x-axis). (b) Resonance frequency of bending mode 2 of the tibia varying with Cut L1 length and Cut L2 length. (c) Resonance frequency of bending mode 3 of the tibia varying with Cut L1 length and Cut L2 length. (d) Resonance frequency of bending mode 4 of the tibia varying with Cut L1 length and Cut L2 length.

In the text

	Figure 7 A COMSOL model geometry showing the two emulated distal and proximal ends of the bone. The scale of the ends is adjusted using SD and SP [Blue Face], and the length of these ends adjusted with LD and LP [Red Lines] respectively.
In the text