The porous cantilever beam as a model for spinal implants: Experimental, analytical and finite element analysis of dynamic properties

1.   Introduction

The human body is frequently exposed to various types of whole-body vibration (WBV) in daily life, e.g., from vibrating transport vehicles or during walking, running, and jumping ^[1]. However, long-term exposure to vibration may cause damage to the spine ^[2], and low back pain and degenerative spinal disorders have been found to be more common among vehicle drivers and other operators of vibrating machines ^[3]. Therefore, many researchers have studied the dynamic response of the spine under the influence of vibration. For instance, in a study by Marini et al. ^[4], human lumbar discs were tested by applying a sinusoidal displacement to the bottom of the disc and a preload to the top. The results showed that the intervertebral disc exhibited non-linear and asymmetric dynamic properties during a continuous frequency sweep and the system showed abrupt changes in vibration amplitude at certain frequencies. Vertebral endplate failure was evident in many specimens after exposure to the frequency sweep. Guo et al. ^[3] used the finite element method to determine the modal vibrational modes of the spine at resonant frequencies. The results showed that vertical oscillation is the predominant mode of human body vibration, with a small amount of motion in the anteroposterior direction. Matsumoto et al. ^[5,6], when investigating vibration of a seated human body, found the bending modes of the entire spine, with the main resonant frequencies corresponding to the vertical motion of the upper part of the trunk relative to the pelvis.

Spinal fusion is a common treatment for intervertebral disc degeneration that decompresses and stabilizes degenerated segments to eliminate the pain. In general, patients are implanted with an intervertebral fusion cage that bears direct axial load, maintains the height of intervertebral and foraminal space, and eventually helps to fuse the adjacent vertebrae together through osseointegration and, ideally, osteoconduction and osteogenesis. These fusion procedures are usually complemented by bilateral pedicle screw instrumentation to increase segmental stability and prevent subsidence and pseudoarthrosis. This instrumentation typically consists of metallic screws interconnected with rods/plates of varying stiffness, depending on material and geometry. Wei et al. ^[7] found that this fixation system significantly reduced the dynamic response of the associated intervertebral disc to the vertical vibration, compared to a spine without fixation, suggesting that the fixation absorbed a large amount of vibration energy after lumbar interbody fusion. However, the disadvantage of fusion surgery is that mobility is limited in these regions of the spine ^[8,9]. As an alternative, dynamic stabilization systems consisting of semi-rigid screws connected to flexible rods have been proposed to align and stabilize the segment, but retain mobility. In recent years, these dynamic stabilization systems have received much attention as they not only maintain some motion but also balance the load distribution between the anterior and posterior columns ^[10].

Given the substantial load-sharing potential of spinal implants, it is important to study the dynamic properties of the implant itself, as its dynamic response can strongly influence the loads transmitted to the reconstructed natural tissues ^[4,11]. Vibration is a basic dynamic response of a structure described as a mechanical oscillation about a stable reference position ^[12]. Structural vibration analysis is a non-destructive testing method to provide dynamic material properties through time-deformation curves, such as the eigenfrequency of the structure. An applied periodic force at the natural frequency can lead to resonance, and a substantial increase in oscillation amplitude. Therefore, it is essential to understand the unique vibration characteristics of an implant to avoid resonance, and the related excessive stresses in the stabilized segment caused under dynamic loading ^[13]. Many studies have shown that porosity, mass and geometric imperfections have an effect on the frequency response of the system ^{[14,15,16,17]}. For instance, Khaniki et al. ^[18] studied the nonlinear forced vibrations of a porous-hyperelastic beam by simultaneously solving the axial and transverse nonlinear coupled equations using a dynamic equilibrium technique and a Galerkin scheme. The results showed that increasing the porosity of a uniformly porous model shifted the resonance peak to lower frequencies and increased the maximum amplitude. However, the type of porosity (uniform or functional along the length) also has a significant effect on changing the nonlinear frequency response of a system, i.e., the stiffness softening behavior could turn into a combination of hardening and softening behavior for the first transverse coordinate with increased porosity.

In addition to the natural frequency, the damping ratio is another important dynamic characteristic of a structure that quantifies the rate of decay and cessation of vibrations. Damping can describe the energy dissipation of a material during vibration or cyclic deformation ^[19]. In many cases, damping is a favorable property of a system or structure. However, it is difficult to control the damping ratio of medical implants by adding additional damping elements (tuned mass dampers), which are widely used in other fields, due to space and complexity constraints. Therefore, possible approaches to adjust the intrinsic damping of the implant include choosing a different material (material damping) or optimizing the geometry of the structure (structural damping). By controlling the damping ratio, unwanted vibrations can be prevented.

In this study, we propose that the introduction of porosity would alter the dynamic response of a spinal implant. The cantilever beam model has been used to investigate the dynamic responses of different porous structures, since the cantilever beam is the most commonly used model for investigating dynamic response under bending and its shape is similar to that of a spinal plate or a rod. This study aims to experimentally measure the natural frequencies and damping of cantilever beams with open holes, as a representative model of bending beam implants, in order to validate an analytical and a numerical model. This model is then used to predict the response of more complex beams with stochastically distributed pores.

2.   Materials and methods

2.1.   Preparation of the PC beam with open holes

For the cantilever beam model, a polycarbonate (PC) test beam with a length (L) of 210 mm, a width (B) of 30 mm, and a thickness (H) of 2 mm was defined, as shown in Figure 1. The modulus and density of PC are 2380 MPa and 1200 kg/m³, respectively. For the experimental comparison, Epraform® PC isotropic sheets (Eriks, The Netherlands) were machined by Burmak AG (Switzerland) into beams with a length of 240 mm including a clamped length of 30 mm at the constrained end. To evaluate the influence of the holes, 10 open holes, i.e., full-thickness, with radii of 2.5, 5.0, and 7.5 mm were drilled through the beam and then water-jetted with abrasive sand, with a solid beam as control object.

2.2.   Analytical modelling

2.2.1.   Free vibration behavior of a cantilever beam

The cantilever beam system can be treated as an idealized single degree of freedom (DOF) system compromising spring, damping, and mass ^[20,21]. We assume that the cantilever beam is a viscous damping system. Based on D'Alembert's principle, we can add the inertial force to the accelerating body and obtain an equivalent static system.

Based on these assumptions, the kinematic equilibrium can be described as follows ^[22]:

where m, c, and k are the mass, damping, and stiffness of the system, respectively, in units of kg, N·s/m, and N/m. And x represents the displacement of the mass point; its first- and second-order derivatives denote the velocity and acceleration.

Since ${w}_{n}^{2} = \frac{k}{m}$ and ${\rm{ \mathsf{ ζ} }} = \frac{c}{2\sqrt{m}\cdot k}$, where ${w}_{n}$ is the undamped natural frequency and ζ is the damping ratio, we can get:

Assuming it is an underdamped system, the solution of this differential equation is:

where A is the amplitude of the movement, ${w}_{d}$ is damping frequency of the system, and $\varphi$ is a phase angle ^[22].

2.2.2.   Mathematical derivation of the natural frequency

In this section, an analytical formulation is presented to describe the effect of open-hole size on the change in natural frequency of a cantilever beam with open hole model. The relationship between the size of the open hole and the dynamic behavior in terms of natural frequency is determined by an analytical approach.

The undamped natural frequency (f) can be derived from the stiffness (k) and the mass (m), as shown in Eq (4). The deflection of a beam (D) is given in Eq (5), where EI is the flexural rigidity and F represents the force applied at the free end ^[23].

Stiffness, which is the force per unit displacement, is given by Eq (6).

Substituting Eq (6) into Eq (4) yields the expression for the undamped natural frequency, Eq (7) ^[23].

However, the cantilever beams in the present study have several open holes that affect the natural frequency by reducing the equivalent stiffness and mass. An equivalent stiffness is introduced by integrating the local stiffness along the hole radius (r) to account for the influence of the holes.

The moment of inertia at the section inside the hole area, ${I}_{c}$, can be expressed by Eq (8). On the other hand, the moment of inertia at the section in the area without a hole, ${I}_{s}$, is determined by Eq (9).

Then the average moment of inertia, ${I}_{A}$, is calculated as follows.

Finally, we obtained the analytical formulation for the natural frequency of the first mode. Based on classical Euler-Bernoulli beam theory and continuous uniform cantilever beam theory, the frequency constant, $\beta$, for the first mode equals to $0.24267$ ^[24,25].

The above formulation can be simplified by Eq (12).

2.3.   Free decay vibration test

The fundamental vibration frequency of the beam was measured experimentally using a non-contact vibration measurement technique, as described in our previous study ^[26]. Briefly, the beam specimen was clamped at one end and the other end was free. Initially, a weight of 40 g was placed on the free end of the beam; a sudden removal of the weight then excited the free vibration of the beam, during which, a high-speed camera (MotionPro Y8-S3, Integrated Design Tools, Ltd, UK) was used to capture the movement of the free end at 500 frames/seconds. The vertical movement of the free end, denoted as y, was extracted using an open software Tracker (open source physics, comPADRE). The equation, $\mathrm{y} = \mathrm{A}{e}^{-{\rm{ \mathsf{ α} }}t}\mathrm{cos}({w}_{d}t+\varphi)$, was used to describe the vibration of the beam, including the exponential decay to the peaks of y data and periodical change. Then, the data was fitted using Curve Fitting function in MATLAB (R2018a, MathWorks, Massachusetts, USA). The exponent, ${\rm{ \mathsf{ α} }}$, represents the decay factor of the free vibration amplitude, which is directly related to the damping ratio of the cantilever beam ^[27]. The captured data was used to determine the damped frequency and the damping ratio of the cantilever beam using the following Eq (13).

where $\alpha$ is the decay factor and ${w}_{d}$ is the damped frequency.

2.4.   Finite element simulation

2.4.1.   Eigenvalue problem for natural frequencies and associated mode shapes

ABAQUS software (ABAQUS 6.4.1, Hibbit, Karlsson and Sorenson Inc., Pawtucket, RI, USA) was used to perform an eigenvalue extraction step to extract the eigenvalues and eigenfrequencies of the PC cantilever beam containing open holes. One end of the beam is fixed (boundary condition) and then the system is subjected to vibration. The model uses a continuum 3D solid element (C3D8, 8-node linear brick). A uniform mesh with an element edge size of about 1 mm was used for the entire structure after a convergence check. The eigenvalue extraction was then performed in ABAQUS using Lanczos method ^[28] to calculate the natural frequencies and the corresponding mode shapes of the system.

2.4.2.   Damped vibration simulation of the open-hole cantilever beam

The concept of Rayleigh damping was utilized to model the damping properties of the vibrating cantilever beam in the associate finite element model. Firstly, the experimentally determined value of the damping ratio of the cantilever beam without open holes was transferred into Rayleigh damping parameters in the finite element model to predict the free vibration decay curve.

Rayleigh damping introduces damping into the vibrating structure in the form of a damping matrix [C], which is a linear combination of the mass matrix [M] and the stiffness matrix [K] of the system ^[29,30], that is:

where α and β are proportional damping coefficients.

The values of α and β are calculated from modal damping ratio (${\xi }_{i}$) and satisfy the following relation ^[31,32]:

where ${w}_{i}$ is the natural frequency of i-th mode. ${\rm{ \mathsf{ α} }}$ is predominant in the low frequency response of the system, while β mainly influences the high frequency phase.

Structural vibration mode is orthogonal on the mass matrix and the stiffness matrix, therefore Rayleigh damping as a linear combination of the mass matrix and stiffness matrix must meet the orthogonal condition ^[31]. The matrix form of the real mode ($\left[\varphi \right]$) orthogonal relation is

And

Then,

Thus,

The method of least squares was used to find the optimal parameter values by minimizing the sum, ${\rm{ \mathsf{ δ} }}$, of squared residuals for each mode shape i ^[33]:

Taking the partial derivative with respect to ${\rm{ \mathsf{ α} }}$ and ${\rm{ \mathsf{ β} }}$ and set them to 0:

Then the solution for α and β is:

We can calculate ${\rm{ \mathsf{ α} }}$ and ${\rm{ \mathsf{ β} }}$ from the above equation using the values of damping ratio and natural frequency. In the formulation of Rayleigh damping, it is generally assumed that the mass-proportional damping effect dominates at the lower frequencies, while the stiffness-proportional damping dominates at the higher frequencies. Since the vibrating cantilever beams in our study corresponded to the first mode and their vibrating frequencies were low, we set the ${\rm{ \mathsf{ β} }}$ value to 0. From the experimental result of the cantilever beam without open holes, we calculated ${\rm{ \mathsf{ α} }}\mathrm{ }\ \mathrm{a}\mathrm{s}$ 0.94. With this parameter, we then chose two different simulation algorithms in ABAQUS, modal dynamic analysis and dynamic explicit, to simulate the vibration behavior of a cantilever beam without and with open holes. The simulation processes are shown in Figure 2.

For the modal dynamic analysis, the first step was to perform a frequency analysis under the option of linear perturbation, where we chose the Lanczos method to calculate the frequency. After Step-1, modal dynamics was applied, where we fixed the Rayleigh parameters. The boundary condition was that one end was fixed and a concentrated force (0.4 N) was applied to the other end at the beginning and then released over 0.3 s.

In the dynamic explicit method, the boundary condition for the initial step was set to be the fixation of one end. In Step-1, a boundary condition for the displacement of the free end of 20 mm was applied over 0.1 s. In Step-2, the displacement boundary from Step-1 was inactivated, allowing the cantilever beam to oscillate freely. The displacement of the free end over time was calculated.

2.4.3.   Damped vibration simulation of the cantilever beam with random pores

To evaluate the dynamic response of more complex structures and the influence of the structure on this response, complex beams with stochastically distributed pores were created based on Voronoi tessellation. Each beam model with random pores has the same porosity as the corresponding beam model with circular open holes. For example, the beam model with 10 open holes of 5 mm diameter (Figure 3(A)) has a porosity (3.12%) which is maintained for the beam model with 10 random pores (Figure 3(B)).

The dynamic explicit method was used to simulate the free vibration response of a cantilever beam model with random pores. The steps specified in the dynamic explicit method were the same as in Section 2.4.2. The mesh element type was from the explicit 3D stress family and the element shape was a 10-node modified quadratic tetrahedron. The approximate global size was 1 mm.

A static analysis was also applied to analyze the deflection and stress distribution of a cantilever beam with random pores, with the same force as in the modal dynamic analysis (0.4 N) applied at the free end.

2.5.   Statistical analysis

For the experimental tests, three samples were created and three separate experiments were carried out in each group to collect the data. The data were expressed as the mean values $\pm$standard deviation (SD) and analyzed by one-way analysis of variance (ANOVA). The level of statistical significance was set at p $\le$0.05. GraphPad Prism 8.2.0 software (GraphPad Software Inc, California, USA) was used for statistical calculations. All error bars correspond to the SD of the mean to indicate the uncertainty for each measurement.

3.   Results

3.1.   Natural frequency

Figure 4 shows the mode shape of a PC cantilever beam with open holes, predicted by the eigenvalue extraction in the finite element simulation. For example, the first vibration mode of the cantilever beam with open holes of 5 mm radius was the bending mode in the horizontal direction. In this mode of vibration, the frequency was 9.76 Hz. The beam tended to bend around the minimum moment of inertia at the root. The second, fifth and sixth vibration modes were all bending modes with natural frequencies of 61.26,171.93, and 338.03 Hz, respectively. However, the individual vibration modes had different bending peak values. The third mode of vibration was twisting about the root with a frequency of 127.36 Hz. The fourth mode of vibration was in-plane shaking with a frequency of 158.09 Hz. Since the first mode is relatively important, we focused on this in our study and investigated the natural frequency of the first mode in detail using analytical, numerical and experimental methods.

The natural frequencies of the first mode of the cantilever beam without open holes derived from the analytical model, the finite element method and the experimental tests were 10.32, 10.50, and 10.38$\pm$0.02 Hz, respectively. It should be noted that the frequencies observed in the experimental tests were damped frequencies, which were slightly lower than the natural frequencies. However, the difference between the natural and damped frequencies can be neglected due to the relatively low damping ratio. For example, the average damped frequency was 10.3838 Hz, while that of the natural frequency computed from the measured damping ratio was 10.3841 Hz. Overall, the natural frequencies obtained by the three methods were in agreement, indicating that both the analytical method and the finite element method can accurately predict the natural frequency of the first mode of the cantilever beams without open holes. For cantilever beams with open holes, the results of each method are shown in Table 1. The differences between each method were calculated. The difference between the analytical, numerical and experimental results was within 2% when the radius of the open hole was less than 5 mm. For a cantilever beam with open holes of 7.5 mm radius, the results of the finite element method still showed a high correspondence with the experimental results. However, the result of the analytical model was 4.5% lower than the experimental results and 5.17% lower than the finite element prediction. The finite element prediction is generally in good agreement with the experimental result.

To observe the variation of the natural frequency resulting from the open hole, the analytical, numerical and experimental results for different sizes of open hole were plotted in Figure 5 using the interpolation method. It can be clearly seen that the natural frequency decreased as the size of the open hole increases.

3.2.   Damping ratio

The damping ratio was calculated from the free vibration decay response of the cantilever beam. Figure 6(A) shows the comparison of the free vibration response of the cantilever beam with open holes of different radii. The initial displacement amplitudes of the different groups of cantilever beams are different because the deflection of the cantilever beam with a constant weight increases with porosity. The envelope curve represents the decay rate of the free vibration. The free vibration amplitudes decayed faster for the cantilever beam with larger open holes. Additionally, the exponent of the envelope curve was also higher for cantilever beams with larger open holes, reflecting the increase in damping ratio. The damping ratios are given in Figure 6(B) with values of 0.73%$\pm$0.010%, 0.76%$\pm$0.025%, 0.78%$\pm$0.007% and 0.83%$\pm$0.003% for cantilever beams with open holes radii of 0, 2.5, 5, and 7.5 mm, respectively. The damping ratio of the cantilever beam increased with the increase of the radius of the open hole, and the damping ratios of all three groups of cantilever beams with open holes were significantly higher than those of cantilever beams without open holes.

Figure 7(A), (C) show the displacement versus time for the cantilever beam during free decay vibration using the two different finite element simulation methods. Figure 7(A) demonstrates the free vibration response of the cantilever beam with the dynamic explicit method. Since the initial displacement of the free end was set to 20 mm for all the groups, all cantilever beams started to vibrate from the same position, and the peak amplitude decayed with time for each cycle. The vibration period of different groups varied with the radius of the open hole. Figure 7(B) shows the damping ratio calculated from the free vibration response. Similar to the experimental results, the damping ratios predicted by the dynamic explicit method presented the same trend, i.e., cantilever beams with larger open holes had higher damping ratios. The decay factor of cantilever beams with 0, 2.5, 5, and 7.5 mm radii of open holes were 0.4693, 0.4695, 0.4699, and 0.4705, respectively, and the calculated damping ratios were 0.71, 0.74, 0.77, and 0.81%, respectively. Figure 7(B) shows the free vibration response of the cantilever beam obtained with the modal dynamics method. Since in this method, we placed the same weight at the free end of the cantilever beams and released the weight afterwards, the initial vibration displacements of the four groups were different due to the different bending stiffness of each structure. The deflections of each group showed good agreement with the experimental values. As shown in Figure 7(B), the free vibration amplitude of the cantilever beam with larger open holes was greatly reduced. The damping ratios of cantilever beams with 0, 2.5, 5, and 7.5 mm radii of open holes were 0.71, 0.73, 0.77, and 0.81%, respectively. However, the decay factor was the same with a value of 0.47.

Figure 7 shows a comparison of the damped frequencies and damping ratios from the experimental tests and the finite element analyses. The results of the two finite element methods were in good agreement (Table 2). However, the numerical damped frequencies were slightly higher and the damping ratios were slightly lower for most groups compared to the experimental results. This may be related to the lower stiffness values of the experimental samples because of unavoidable defects in the real beams. Nevertheless, the FE predictions of damped frequencies and damping ratios were in good agreement with the experimental observations. Additionally, the damping ratio was found to be more sensitive to the beam structure (i.e., hole geometry) than the damped frequency, since the damping ratio of the cantilever beam with open holes of 7.5 mm radius decreased by 14% compared to the cantilever beam without open holes, while the damped frequency decreased by 12%.

Figure 8 shows the comparison of damped frequencies, damping ratios and decay factor of the cantilever beam with regular holes and random pores by dynamic explicit simulation. The damping frequency of the cantilever beam with regular holes decreased with the increase of the open-hole diameter and the damping ratio increased, while the damped frequency and damping ratio of the cantilever beam with random pores showed the opposite trend. The difference in damped frequency increased dramatically as the porosity increased. As the damped frequency of a structure depends on its mass and stiffness, quasi-static bending simulations were performed for the cantilever beams with random pores, and their deflection and bending modulus are summarized in Table 3. The equivalent bending modulus decreased by only 13.39% as the porosity increased from 0% to 28.09%. The stress distribution diagram shows that the neutral layer was subjected to the least stress (Figure 9). This indicates that the random pores in the neutral layer had a greater effect on the mass of the cantilever beam structure than on the stiffness.

4.   Discussion

Considering the natural frequency, modal analysis is an effective method to calculate and visualize the complex deformation and the dynamic properties of a structure, by means of natural frequencies and their associated modal shapes ^[13,34]. In this paper, the natural frequencies of the structure, including both the original cantilever beam and cantilever beam with open holes, were estimated by discretizing the structure into elements. The n-order natural frequencies were obtained from the stiffness and mass matrices by calculating the eigen vectors and eigen values with commercial finite element software.

The reduction of equivalent mass and the reduction of stiffness are two mechanisms that have opposite effects on the natural frequency. In the mathematical model of the natural frequency of the cantilever beam with open holes, we considered the reduction of stiffness by recalculating the moment of inertia resulting from the open holes, as the total volume of the cantilever beam structure is the combination of the solid and open regions. The mathematical model is intuitive and can help to determine the natural frequencies of such beams efficiently. The proposed analytical model gives similar results compared to the finite element model and the experimental tests, with very low computational cost. However, models based on the continuous beam theories would allow for a more accurate modeling of the structure and realize a more precise estimation of the natural frequency, particularly for higher vibration mode. Finite element simulation was shown to accurately predict the vibrational response of the cantilever beams with open holes, validating its use for analysis of more complex structures. All three methods reflect a decrease in the natural frequency of the first mode as the size of open hole increases, which is in good agreement with the well-accepted perception that lower stiffness leads to a decrease in the natural frequency ^[35]. Therefore, the natural frequency test can also be used to predict the stiffness of a material or structure.

The techniques used to determine the damping properties of materials and structures can be divided into two categories. The first category comprises direct methods that measure energy dissipation directly, such as energy, thermal and hysteresis loop methods. The other category comprises indirect methods measuring the amplitude and frequency associated with energy dissipation, including free damped vibration and resonance curve/half-power bandwidth methods ^[36]. Of these methods, the free damped vibration method is the simplest and most efficient method for accessing the dynamic response of cantilever beams. In the present study, the free vibration responses of all cantilever beams were agreement with the theoretical analyses. Thus, the damped frequency, decay factor, and damping ratio can be extracted from the time histories of the end displacement. This free vibration method does not require a complex setup or vibration generators and is therefore easy to put into practice.

The decay observed in the free vibration response of a structure reflects the damping. This damping phenomenon usually arises from energy dissipation caused by friction, hysteresis and possibly viscoelastic effects in both structural and non-structural members ^[37]. Environmental effects, such as air damping, are an external source of damping for the structure and alter its vibration amplitude as well as shift the resonant frequency ^[38]. In our study, all cantilever beams were tested experimentally under the same external conditions (i.e., lab environment). As for the finite element simulations, air damping was not considered. This may explain why the experimental damping ratios were slightly higher than the simulated damping ratios, especially for cantilever beams without open holes. However, the internal material damping as well as the structure changes are of prime importance.

Material damping includes interfacial friction between the phases of the material itself, as well as damping caused by energy dissipation through molecular and physical bonds ^[39,40]. Typically, materials with large free volumes, loosely packed molecules, weak intermolecular attractions, amorphous nature and flexible molecules have higher damping. Polymers are typically viscoelastic, exhibiting both elastic and viscous properties simultaneously. When a PC cantilever beam vibrates, part of the energy is stored (elastic) and part is dissipated as heat (viscous) ^[36]. In our study, as the stiffness of the cantilever beam decreased with increasing radius of open holes, the beam tended to be more compliant, thus also enhancing its damping properties. Furthermore, literature suggests that damping changes due to the existence of porosity is noteworthy ^{[41,42,43,44]}. The damping due to heat generation and irreversible heat flow in linear, isotropic and homogeneous thermoelastic rectangular plates including uniformly distributed cavities has been described by Stamatopoulos et al. ^[41]. The results showed that damping increased with increasing porosity through a nearly linear relationship. Li et al. ^[42] also presented that a porous magnesium with 3D entangled pore structure exhibited significant damping capacity and a higher porosity and a smaller pore size contributed to larger loss factor, although the physical behavior of metals and polymers are quite different. Golovin et al. ^[44] proposed that damping in porous structures was enhanced by localized stresses in comparison with the corresponding dense materials. A model for mechanical damping in porous materials was suggested for the deformation on the basis of statistical mechanics of micro-heterogeneous materials.

Although we have demonstrated an analytical solution for the beams with regular geometry, many engineering problems are not easily solved by analytical methods. Closed-form solutions of the differential equations describing the physics of the problem are difficult to obtain because of the non-linearity of the materials, the geometrical complexity of the structures and the discontinuity of the structure. To overcome these difficulties, finite element methods are used to obtain approximate solutions of the set of differential equations ^[29]. In order to successfully implement finite element simulations, the material properties, geometry and type of loading configuration should be accurately defined. In our study, we used two finite element methods, dynamic explicit and modal dynamics, to simulate the free damped vibration behavior of a cantilever beam. The explicit dynamic integration method is a mathematical technique for integrating the equations of motion through time and is also known as the central difference or forward Euler algorithm. Often, this method can effectively handle large-scale models with high loading velocities. In addition to the direct integration method, the mode superposition method can also be used to calculate the dynamic response of a structure. Rayleigh damping is a mathematically convenient concept and a linear model, and is still the most popular choice for damping modelling in linear and even non-linear analysis due to its computational efficiency and ease of implementation on commercial software platforms ^[37]. Finite element results demonstrate that the damping of cantilever beams can be modeled effectively and accurately using Rayleigh damping. Furthermore, the mass-proportional damping coefficient, which is part of the Rayleigh damping formulation, has been shown to satisfactorily model the experimentally observed damping response of cantilever beams. As a result, mass-proportional damping effects dominate in low frequencies. Moreover, in dynamic explicit simulations, the damping properties can be seen as material properties, related to the material itself. In modal dynamics, on the other hand, damping is regarded as a purely numerical concept that is implemented upon the structure and intrinsically related to the structural properties of the vibrating system. Therefore, the dynamic explicit method is more suitable for predicting the dynamic response of complex structures if the damping properties of the material are given.

An understanding of the vibrational characteristics of medical implants is essential to guide the selection of appropriate prosthesis configurations to prevent excessive stress. In fact, there are already some examples in the field of prosthodontics that use vibrational features of dental implant prostheses to estimate structural weakness ^[13]. Pedicle screw plate fixation is an effective form of immobilization of the spine used to achieve arthrodesis ^[45]. The pedicle plate is a medically designed implant that is used to provide spinal stability in spinal instrumentation and fusion procedures. In our study, the beam model with open holes has a similar shape to the bone plate. Elucidation of the vibrational characteristics of bone plates can be used to estimate structural weaknesses, in which case stresses accumulate from vibration loading at specific frequencies. In addition, the damping properties of bone plates is essential to reduce the risk of significant vibration and sudden failure of implants. Therefore, engineers can adjust the number and size of open holes to achieve the desired stiffness and damping distribution in the bone plate.

Recently, topology optimization has attracted much attention in the development of medical devices. It is a computational method for optimizing the structure, i.e., by allocating material within a prescribed design domain, according to the given external load and boundary conditions, under the premise that the constraints such as displacement, stress and balance are satisfied ^[46]. Topology optimization is an attractive method for (generally quasi-static) spinal implant design with the goal of best adapting the implant's mechanical response to the requirements of its application. However, topology-optimized designs usually have complex geometric shapes, making it difficult to experimentally determine their dynamic response, especially their mechanical behavior in vivo. In our study, we established a heterogeneous porous beam model as a representative case of an irregular structure. The random pores in the beam were created by Voronoi tessellation which is widely used to generate bone-like structures due to its simple definition, anisotropic properties, controllable by its seeds, etc. ^[47]. The observation that the pores in the model had a greater effect on the mass than on the stiffness of the structure, due to their spatial distribution, suggests that the bending modulus of cantilever beams with substantial porosity can be preserved if the pore distribution is carefully considered. The observation that the natural frequency increased with increasing porosity for the random pore structure, which is the opposite of the response observed for beams with regular hole distribution, further highlights novel opportunities to tailor the static bending stiffness and dynamic response of beam implants through the structure.

5.   Conclusions

In this paper, the dynamic response of PC cantilever beams with open holes was systematically studied, as a model system representing spinal implants experiencing bending loads. Through analytical modelling, analysis and experimental verification, it was shown that the characteristics of the vibration response, such as mode shapes, frequencies and displacement amplitudes, are related to the stiffness of the structure. The analytical model of natural frequency is simplified, yet can accurately predict the natural frequency of beam structures with regular open holes. Dynamic explicit and modal dynamics finite element methods can correctly simulate the damped vibration behavior of beams with open holes by applying the damping properties of the material, and can be applied to analyze the damping behavior of complex, non-regular structures. Natural frequencies and associated modes, and damping properties, were shown to vary with hole or random pore structure, highlighting potential paths to tailor the dynamic response of complex implant structures.

Acknowledgments

This project has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No.812765.

Conflict of interest

There is no conflict of interest.