ME4508 - Project - Solved

The geometry of a femur bone (data from https://grabcad.com/library/femur-bone-2). The bone is roughly aligned in the x-direction. 

Mechanical analysis of bones is extremely important. For example, many people undergo total knee replacement. According to the American Joint Replacement Registry Annual Report data, in the United States, 0.3–5.5% for over 600,000 primary total knee replacement surgeries may eventually lead to post-surgical bone fractures. One reason is that the interactions between the bone and the prosthetic knee could cause excessive stresses in the bone and lead to its failure. As such, it is important for us to establish a baseline for the stress-strain responses of specific bone samples so that we can model the responses of the bone in the presence of medical devices more accurately.  

Obviously, this assignment is a course project and not a PhD dissertation. So, we will make significant simplifying assumptions. For example, we assume that the bone is homogenous and linear elastic (with small strain assumption). Obviously, bones are not homogenous, but the linear elastic assumption is fairly valid.  

Let us also assume that we have cut a section with the length of 205 mm from the femur bone shown in Fig. 1. This section has a changing cross-sectional area as shown in Fig. 2. If you are curious, the area was calculated using CT-scans and later segmented data in SolidWorks. I have provided the raw data used in this graph in a file called “crossSectionData.csv”, which is available in the course CANVAS.  

•      Part (a) Import this file into MATLAB using readmatrix command and find the value of the averaged cross-sectional area for this sample. Plot area versus distance from the tip and confirm that your graph matches Fig. 2. You can use the command “print” to export your MATLAB figures into high resolution images that you need to include in your report. Do not use print screen. This section has 5 points.

 

 

Figure 2. Changes in the cross-sectional area of the bone sample used in our mechanical testing experiment. 

 

 

Figure 3. Force-displacement relationship obtained from the compression mechanical testing of a section of the femur bone. 

We then used a universal testing device (like an Instron Machine) to find the relationship between the displacement 𝛿 and the force F for this sample, as shown in Fig. 3. We define the apparent stiffness 𝑆! of this sample as the slope of a line that best represents data in Fig. 3. Since displacement should be zero when no force is applied, the apparent stiffness at each step of the experiment should be

                                                                                             𝑆! =              "#                                                                                                                                                   (1)

•      Part (b) Import the data of Fig. 3 into MATLAB using the file “ForceDisplacement.csv” that I have provided for you in CANVAS and using readmatrix command.  Use “polyfit” command in MATLAB to find the slope of the line 𝑆! that best fit the force-displacement data. Use that slope value of 𝑆!   (assume that the intercept is zero) and plot a line over the raw data.  This section has 10 points.

Now let use first assume that the bone is simply a single bar with the average area of A that you calculated in part (a). As such we can simply write:

                                                                                        𝐹 =  $% 𝛿                                                                     (2)

&

or  

                                                                                       𝑆! =  $%&                                                                      (3)

From equation (3), we can simply estimate the elastic modulus of the bone.  

•      Part (c) Use equation (3), L= 205 mm, and A and 𝑆! that calculated in previous sections to find E for this sample. This section has 5 points.

  

However, as you know, the cross-sectional area of the bone is not constant. So, we will use our fancy finite-element knowledge to more accurately estimate the value of the elastic modulus. Here is how we are going to do it. First let us assume that we magically know that E = 10 GPa. We will use finite element method and model the bone as a series of bars connected to each other. For boundary conditions, we assume that the bar is fixed at x=0, and has a compressive load of –1200 N applied on it at x=205 mm. Using the magical value of E = 10 GPa, we would like to calculate the displacement 𝛿 at x=205 mm and then use that for calculation of 𝑆"%$ which is the ratio of the reaction force to the displacement at x=205 mm.  

•      Part (d) Write a MATLAB program and calculate 𝑆"%$, as you increase your element numbers from 1 to 100. For every step that you increase the number of elements, use the average of the area that belongs to that section of the bone as the corresponding area of that element. Plot the 𝑆"%$ versus the number of elements. Find the smaller number of elements, for which 𝑆"%$ does not change more than 0.2% in two consecutive simulations. This section has 50 points.  

•      Part (e) From part d, you should be able to calculate 𝑆"%$ and now you should be able to compare that with 𝑆! that you calculated in part b. Do these two numbers match? Explain why in no more than three sentences. This section has 5 points.  

Now from hereon, it makes sense to only run models with the appropriate number of elements that you identify in part (d). That number might be slightly different from other students depending on how you calculated your errors. So, for consistency, let us choose ten elements from hereon.

•      Part (f) Run your finite element model in a loop with values of 𝐸      𝐺𝑃𝑎, 25           𝐺𝑃𝑎] using 0.1       𝐺𝑃𝑎 intervals. At each step, calculate the error defined as  

 𝑒𝑟𝑟𝑜𝑟 =          5(𝑆"%$ − 𝑆!)’

Plot error versus the values of E. Which value of the E minimizes the error, i.e., the difference between the experiments and the model? Is this value of the E different from what you calculated in part (c)? This section has 15 points.  

•      Part (g) Calculate the element stresses for the –1200 N compressive load using the E that you calculated in part (f). Which element has the largest stress value? Does this make sense? This section has 5 points.   

•      Part (h) Write your reflection about this project. What did you learn? What can I do to improve the projects in future years? This section has 5 points. 

•      Part (i) Instead of using a loop with predefined values of E, use either lsqcurvefit function, fminsearch function, or ga function to find the value of E that minimizes the error. This section has 10 bonus points.