Femur shaft (Ivarsson 2009)¶
Validation model information¶
Postprocessing of femur shaft validations based on Ivarsson et al. 2009
- Performed by: Nico Erlinger
- Reviewed by: Corina Klug Added to VIVA+ Validation Catalog on: 2022-10-21
| Version | Date | Performed by | LS-Dyna |
|---|---|---|---|
| 0.3.2 | 2022-10-19 | Nico Erlinger | 9.3.1 |
| 1.1.0 | 2024-05-22 | Corina Klug | 9.3.1 |
© 2019-2024, OpenVT Organization (OVTO)
The Jupyter notebooks are licensed under Creative Commons Attribution 4.0 International License 
Reference¶
A. Schubert, N. Erlinger, C. Leo, J. Iraeus, J. John, C. Klug (2021): "Development of a 50th Percentile Female Femur Model", IRCOBI conference proceedings, http://www.ircobi.org/wordpress/downloads/irc21/pdf-files/2138.pdf
Experiments by Ivarsson et al. (2009)¶
Ivarsson, B. J.; Genovese, D.; Crandall, J. R.; Bolton, J. R.; Untaroiu, C. D.; Bose, D. The Tolerance of the Femoral Shaft in Combined Axial Compression and Bending Loading in Stapp Car Crash Journal, Vol. 53, 2009.
Summary¶
In total 42 tests of femur shaft were replicated by prescribing the vertical motion to the impactor (approx. 1.5 m/s). In addition, an axial force (approx. 4, 8, 12 or 16 kN) was applied for combined tests (axial compression and three-point bending).
{admonition}
Female specimens (n=16) are listed on the left side and male specimens (n=23) on the right side.
| Specimen | Loading type | Loading direction | Side | Donor age |
|----------|--------------|-------------------|------|-----------|
| 1.18 | co | AP | L | 64 |
| 1.19 | co | PA | R | 64 |
| 1.20 | co | PA | L | 40 |
| 1.21 | co | AP | R | 40 |
| 1.22 | co | AP | L | 45 |
| 1.23 | co | AP | R | 45 |
| 1.24 | co | AP | L | 60 |
| 1.26 | co | AP | L | 50 |
| 1.27 | co | AP | R | 50 |
| 1.28 | co | AP | L | 56 |
| 1.29* | co | PA | R | 56 |
| 1.32 | be | PA | L | 52 |
| 1.33 | be | AP | R | 52 |
| 1.35 | be | PA | L | 62 |
| 1.36 | be | PA | R | 62 |
| 2.10 | be | AP | L | 57 |
| 1.01 | co | PA | L | 51 |
| 1.02 | co | PA | R | 51 |
| 1.03 | co | PA | L | 62 |
| 1.04 | co | PA | R | 62 |
| 1.05 | co | AP | L | 62 |
| 1.06 | co | PA | R | 62 |
| 1.07 | co | PA | L | 49 |
| 1.08 | co | AP | R | 49 |
| 1.09 | co | AP | R | 62 |
| 1.10 | co | PA | L | 44 |
| 1.11 | co | PA | R | 44 |
| 1.12 | co | AP | L | 58 |
| 1.13 | co | AP | R | 58 |
| 1.14 | co | AP | L | 65 |
| 1.15 | co | PA | R | 65 |
| 1.16 | co | AP | L | 53 |
| 1.17 | co | AP | R | 53 |
| 1.30 | co | PA | L | 62 |
| 1.34 | be | AP | R | 63 |
| 1.37 | be | AP | L | 45 |
| 1.38 | be | PA | R | 45 |
| 2.06 | be | AP | L | 39 |
| 2.07b | be | PA | R | 51 |
Abbreviations:
- co : combined (axial compression and three-point bending
- be : three-point bending (no axial compression)
- L : left side
- R : right side
- AP : anterior-posterior
- PA : posterior-anterior
Loading and Boundary Conditions¶
The displacement-time histories from the diagrams published by Ivarsson et al. (2009) were digitised with WebPlotDigitizer v4.4 (https://automeris.io/WebPlotDigitizer).
Experimental responses¶
The force-time curves of the impactor and the axial loadcell were downloaded from the NHTSA Biomchanics Test Database (www.nhtsa.gov/research-data/research-testing-databases).
In [9]:
Copied!
simulation_list = ['1_18', '1_19', '1_20', '1_21', '1_22', '1_23', '1_24', '1_26', '1_27', '1_28',
'1_29', '1_32', '1_33', '1_35', '1_36', '2_10', '1_01', '1_02', '1_03', '1_04',
'1_05', '1_06', '1_07', '1_08', '1_09', '1_10', '1_11', '1_12', '1_13', '1_14',
'1_15', '1_16', '1_17', '1_30', '1_34', '1_37', '1_38', '2_06', '2_07b']
date = datetime.date.today().strftime('%Y-%m-%d')
dynasaur_output_file_name = 'Dynasaur_output.csv'
processed_data_dir = f'data/processed/{date}'
#processed_data_dir=f'data/processed/240522'
simulation_list = ['1_18', '1_19', '1_20', '1_21', '1_22', '1_23', '1_24', '1_26', '1_27', '1_28',
'1_29', '1_32', '1_33', '1_35', '1_36', '2_10', '1_01', '1_02', '1_03', '1_04',
'1_05', '1_06', '1_07', '1_08', '1_09', '1_10', '1_11', '1_12', '1_13', '1_14',
'1_15', '1_16', '1_17', '1_30', '1_34', '1_37', '1_38', '2_06', '2_07b']
date = datetime.date.today().strftime('%Y-%m-%d')
dynasaur_output_file_name = 'Dynasaur_output.csv'
processed_data_dir = f'data/processed/{date}'
#processed_data_dir=f'data/processed/240522'
Results¶
<>:103: SyntaxWarning: "is not" with a literal. Did you mean "!="?
<>:103: SyntaxWarning: "is not" with a literal. Did you mean "!="? C:\Users\klugcor\AppData\Local\Temp\ipykernel_9508\106710733.py:103: SyntaxWarning: "is not" with a literal. Did you mean "!="? if i is not '1_29':
Find values at time of evaluation¶
In [12]:
Copied!
def find_value_at_time_of_evaluation(value_x, value_y):
list = []
counter = 0
for i in simulation_list:
processed_data_path = os.path.join(processed_data_dir, i).replace('\\', '/')
simData = pd.read_csv(os.path.join(processed_data_path, dynasaur_output_file_name),
delimiter=';', na_values='-', header = [0,1,2,3])
x = simData[value_x]
y = simData[value_y]
x = np.array(x).flatten()
y = np.array(y).flatten()
res = np.interp(time_for_strain_eval[counter], x, y)
res = np.round(res, decimals=5)
list.append(res)
counter += 1
return list
value_x = 'BONES', 'Femur_Cortical_L_MPS_history', 'time'
value_y = 'BONES', 'Femur_Cortical_L_MPS_history', 'strain'
MPS_list = find_value_at_time_of_evaluation(value_x, value_y)
value_x = 'BONES', 'Femur_Cortical_L_PS99_history', 'time'
value_y = 'BONES', 'Femur_Cortical_L_PS99_history', 'strain'
PS99_list = find_value_at_time_of_evaluation(value_x, value_y)
value_x = 'FEMUR', 'Impactor_rcforc_z', 'time'
value_y = 'FEMUR', 'Impactor_rcforc_z', 'force'
force_list = find_value_at_time_of_evaluation(value_x, value_y)
value_x = 'FEMUR', 'Impactor_nodout_z', 'time'
value_y = 'FEMUR', 'Impactor_nodout_z', 'displacement'
displ_list = find_value_at_time_of_evaluation(value_x, value_y)
displ_list = [i * (-1) for i in displ_list]
sheet1 = ipysheet.sheet(rows=len(simulation_list), columns=6, column_headers=False, row_headers=False)
column0 = column(0, simulation_list)
column1 = column(1, time_for_strain_eval, numeric_format = '0.0', read_only=True)
column2 = column(2, force_list, numeric_format = '0.00', read_only=True)
column2 = column(3, displ_list, numeric_format = '0.00', read_only=True)
column4 = column(4, MPS_list, numeric_format = '0.00000', read_only=True)
column5 = column(5, PS99_list, numeric_format = '0.00000', read_only=True)
sheet1.column_headers=['Simulation', 'Time of evaluation [ms]', 'Force at evaluation [kN]',
'Displacement at evaluation [mm]', 'MPS [-]', 'PS99 [-]']
display(sheet1)
sheet2=ipysheet.sheet(rows=len(simulation_list), colums=2, column_headers=False, row_headers=False)
column0 = column(0, simulation_list)
column1 = column(1, MPS_list, numeric_format='0.00000')
column2 = column(2, PS99_list, numeric_format='0.00000')
sheet2.column_headers=['Simulation', 'MPS', 'PS99']
df = ipysheet.to_dataframe(sheet2)
# include male and female, exclude outliner
exclude_simulation_list = ["1_10", "1_29", "1_30","1_32", "1_33"]
# for female IRC: exclude male from exported strain sheet
#exclude_simulation_list = ["1_01", "1_02", "1_03", "1_04", "1_05", "1_06", "1_07", "1_08", "1_09", "1_10", "1_11", "1_12", "1_13", "1_14", "1_15", "1_16", "1_17", "1_30", "1_34", "1_37", "1_38", "2_06", "2_07b"]
#for male IRC: exclude female from exported strain sheet
#exclude_simulation_list = ["1_18", "1_19", "1_20", "1_21", "1_22", "1_23", "1_24", "1_26", "1_27", "1_28", "1_29", "1_32", "1_33", "1_35", "1_36", "2_10"]
for i in range(len(exclude_simulation_list)):
df.drop(df.loc[df['Simulation'] == exclude_simulation_list[i]].index, inplace=True)
df.to_csv('data/processed/strain_sheet_for_IRC.csv', index=None, sep=';', mode='w')
def find_value_at_time_of_evaluation(value_x, value_y):
list = []
counter = 0
for i in simulation_list:
processed_data_path = os.path.join(processed_data_dir, i).replace('\\', '/')
simData = pd.read_csv(os.path.join(processed_data_path, dynasaur_output_file_name),
delimiter=';', na_values='-', header = [0,1,2,3])
x = simData[value_x]
y = simData[value_y]
x = np.array(x).flatten()
y = np.array(y).flatten()
res = np.interp(time_for_strain_eval[counter], x, y)
res = np.round(res, decimals=5)
list.append(res)
counter += 1
return list
value_x = 'BONES', 'Femur_Cortical_L_MPS_history', 'time'
value_y = 'BONES', 'Femur_Cortical_L_MPS_history', 'strain'
MPS_list = find_value_at_time_of_evaluation(value_x, value_y)
value_x = 'BONES', 'Femur_Cortical_L_PS99_history', 'time'
value_y = 'BONES', 'Femur_Cortical_L_PS99_history', 'strain'
PS99_list = find_value_at_time_of_evaluation(value_x, value_y)
value_x = 'FEMUR', 'Impactor_rcforc_z', 'time'
value_y = 'FEMUR', 'Impactor_rcforc_z', 'force'
force_list = find_value_at_time_of_evaluation(value_x, value_y)
value_x = 'FEMUR', 'Impactor_nodout_z', 'time'
value_y = 'FEMUR', 'Impactor_nodout_z', 'displacement'
displ_list = find_value_at_time_of_evaluation(value_x, value_y)
displ_list = [i * (-1) for i in displ_list]
sheet1 = ipysheet.sheet(rows=len(simulation_list), columns=6, column_headers=False, row_headers=False)
column0 = column(0, simulation_list)
column1 = column(1, time_for_strain_eval, numeric_format = '0.0', read_only=True)
column2 = column(2, force_list, numeric_format = '0.00', read_only=True)
column2 = column(3, displ_list, numeric_format = '0.00', read_only=True)
column4 = column(4, MPS_list, numeric_format = '0.00000', read_only=True)
column5 = column(5, PS99_list, numeric_format = '0.00000', read_only=True)
sheet1.column_headers=['Simulation', 'Time of evaluation [ms]', 'Force at evaluation [kN]',
'Displacement at evaluation [mm]', 'MPS [-]', 'PS99 [-]']
display(sheet1)
sheet2=ipysheet.sheet(rows=len(simulation_list), colums=2, column_headers=False, row_headers=False)
column0 = column(0, simulation_list)
column1 = column(1, MPS_list, numeric_format='0.00000')
column2 = column(2, PS99_list, numeric_format='0.00000')
sheet2.column_headers=['Simulation', 'MPS', 'PS99']
df = ipysheet.to_dataframe(sheet2)
# include male and female, exclude outliner
exclude_simulation_list = ["1_10", "1_29", "1_30","1_32", "1_33"]
# for female IRC: exclude male from exported strain sheet
#exclude_simulation_list = ["1_01", "1_02", "1_03", "1_04", "1_05", "1_06", "1_07", "1_08", "1_09", "1_10", "1_11", "1_12", "1_13", "1_14", "1_15", "1_16", "1_17", "1_30", "1_34", "1_37", "1_38", "2_06", "2_07b"]
#for male IRC: exclude female from exported strain sheet
#exclude_simulation_list = ["1_18", "1_19", "1_20", "1_21", "1_22", "1_23", "1_24", "1_26", "1_27", "1_28", "1_29", "1_32", "1_33", "1_35", "1_36", "2_10"]
for i in range(len(exclude_simulation_list)):
df.drop(df.loc[df['Simulation'] == exclude_simulation_list[i]].index, inplace=True)
df.to_csv('data/processed/strain_sheet_for_IRC.csv', index=None, sep=';', mode='w')
Sheet(cells=(Cell(column_end=0, column_start=0, row_end=38, row_start=0, squeeze_row=False, type='text', value…
--------------------------------------------------------------------------------------
MPS
--------------------------------------------------------------------------------------
Results from Fit_Weibull_2P (95% CI):
Analysis method: Maximum Likelihood Estimation (MLE)
Optimizer: TNC
Failures / Right censored: 34/0 (0% right censored)
Parameter Point Estimate Standard Error Lower CI Upper CI
Alpha 0.0317518 0.00203919 0.0279963 0.036011
Beta 2.82072 0.379688 2.16662 3.67229
Goodness of fit Value
Log-likelihood 106.119
AICc -207.851
BIC -205.185
AD 0.671476
------Lower---------- Results from Fit_Weibull_2P (95% CI): Analysis method: Maximum Likelihood Estimation (MLE) Optimizer: L-BFGS-B Failures / Right censored: 9998/0 (0% right censored) Parameter Point Estimate Standard Error Lower CI Upper CI Alpha 0.027331 0.000110258 0.0271157 0.0275479 Beta 2.6031 0.0208747 2.56251 2.64434 Goodness of fit Value Log-likelihood 31859.5 AICc -63714.9 BIC -63700.5 AD 13.3172 ------Upper---------- Results from Fit_Weibull_2P (95% CI): Analysis method: Maximum Likelihood Estimation (MLE) Optimizer: L-BFGS-B Failures / Right censored: 9998/0 (0% right censored) Parameter Point Estimate Standard Error Lower CI Upper CI Alpha 0.0369252 0.000127483 0.0366762 0.0371759 Beta 3.06022 0.0230068 3.01546 3.10565 Goodness of fit Value Log-likelihood 30485.8 AICc -60967.7 BIC -60953.2 AD 20.7649
--------------------------------------------------------------------------------------
PS99
--------------------------------------------------------------------------------------
Results from Fit_Weibull_2P (95% CI):
Analysis method: Maximum Likelihood Estimation (MLE)
Optimizer: TNC
Failures / Right censored: 34/0 (0% right censored)
Parameter Point Estimate Standard Error Lower CI Upper CI
Alpha 0.0274591 0.00178048 0.0241821 0.0311803
Beta 2.79418 0.376437 2.14574 3.63856
Goodness of fit Value
Log-likelihood 110.804
AICc -217.222
BIC -214.556
AD 0.673409
------Lower---------- Results from Fit_Weibull_2P (95% CI): Analysis method: Maximum Likelihood Estimation (MLE) Optimizer: L-BFGS-B Failures / Right censored: 9998/0 (0% right censored) Parameter Point Estimate Standard Error Lower CI Upper CI Alpha 0.0236013 9.61276e-05 0.0234136 0.0237905 Beta 2.57831 0.0206764 2.5381 2.61915 Goodness of fit Value Log-likelihood 33252.7 AICc -66501.5 BIC -66487 AD 13.3444 ------Upper---------- Results from Fit_Weibull_2P (95% CI): Analysis method: Maximum Likelihood Estimation (MLE) Optimizer: TNC Failures / Right censored: 9998/0 (0% right censored) Parameter Point Estimate Standard Error Lower CI Upper CI Alpha 0.0319805 0.000111448 0.0317628 0.0321996 Beta 3.03176 0.0227917 2.98742 3.07677 Goodness of fit Value Log-likelihood 31846.9 AICc -63689.7 BIC -63675.3 AD 20.8199
