Rib Anterior Posterior Bending (Kang 2020)#
Model validation information
Performed by : Johan Iraeus
Reviewed by :
Added to VIVA+ Validation Catalog on : 2022-11-22
Last modified :
Model version (this notebook run for): 0.3.2
© 2019-2023, OpenVT Organization (OVTO)
Available openly under under Creative Commons Attribution 4.0 International License 
Experiment by Kang et al. (2020)#
Summary:#
The simulated outputs are compared to the references from PMHS tests reported by Kang et al. (2020) 1
1 Kang, Yun-Seok, Hyun Jung Kwon, Jason Stammen, Kevin Moorhouse, and Amanda M Agnew. 2020. ‘Biomechanical Response Targets of Adult Human Ribs in Frontal Impacts’, Annals of Biomedical Engineering: 1-12.
Experiment#
Information on the subjects/specimens used for the validation#
261 ribs from 171 PMHS (110 male; 61 female) ranging in age from 22 to 99 tyears (age 59±19) were used to create response corridors for age groups 22-40y, 41-60y, and 61+ years [young, middle, older]
The sample included left and right ribs 4-7
Loading and Boundary Conditions#
The ribs were loaded in the anterior-posterior direction at approximately 2m/s. The bracket at the anterior end was given a prescribed displacement based on the recorded displacement in the physical test. The bracket at the posterior end was fixed. Inside each bracket a potting was mounted on a rotation axis.
Positioning#
The ribs were positioned in the pottings in a similar way as in the physical test. It was assumed that the end of the ribs were at the bottom of the potting, see figure above.
Responses recorded#
The response corridors, representing the middle aged groups was digitized from the paper.
Other Notes for simulation#
Notes on implementation and iterations during the implementation of validation simulations
The force is plotted relative to normalized displacements. Normalization is to the initial length of each rib. Norm disp=100%* displacement/Initial length Scale factor = 100%* 1/Initial length
Rib level |
Sex |
Initial length |
Scale factor |
|---|---|---|---|
Rib 4 |
Female |
159 |
0.629 |
Rib 5 |
Female |
172 |
0.582 |
Rib 6 |
Female |
182 |
0.549 |
Rib 7 |
Female |
193 |
0.518 |
Rib 4 |
Male |
176 |
0.568 |
Rib 5 |
Male |
193 |
0.518 |
Rib 6 |
Male |
207 |
0.483 |
Rib 7 |
Male |
217 |
0.461 |
Model Energies#
The figure below the model energies are plotted for each rib.
Show code cell source
# fig_energy, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2,2)
fig_energy, ((ax1, ax2, ax3, ax4), (ax5, ax6, ax7, ax8)) = plt.subplots(2,4, figsize=(10,10))
fig_energy.suptitle('Model Energies', y=1.02)
#ax.set_title('Simulation #1')
ax1.set_ylabel('Energies (J)')
ax1.set(title='50F Rib4', xlabel='time (ms)', ylim=(0,10))
ax2.set(title='50F Rib5', xlabel='time (ms)', ylim=(0,10))
ax3.set(title='50F Rib6', xlabel='time (ms)', ylim=(0,10))
ax4.set(title='50F Rib7', xlabel='time (ms)', ylim=(0,10))
ax5.set_ylabel('Energies (J)')
ax5.set(title='50M Rib4', xlabel='time (ms)', ylim=(0,10))
ax6.set(title='50M Rib5', xlabel='time (ms)', ylim=(0,10))
ax7.set(title='50M Rib6', xlabel='time (ms)', ylim=(0,10))
ax8.set(title='50M Rib7', xlabel='time (ms)', ylim=(0,10))
ax1.plot(simData_50F_4.MODEL.Hourglass_Energy.time, simData_50F_4.MODEL.Hourglass_Energy.energy, label="Hourglass Energy")
ax1.plot(simData_50F_4.MODEL.Internal_Energy.time, simData_50F_4.MODEL.Internal_Energy.energy, label="Internal Energy")
ax1.plot(simData_50F_4.MODEL.Kinetic_Energy.time, simData_50F_4.MODEL.Kinetic_Energy.energy, label="Kinetic Energy")
ax1.plot(simData_50F_4.MODEL.Total_Energy.time, simData_50F_4.MODEL.Total_Energy.energy, label="Total Energy")
ax1.plot(simData_50F_4.MODEL.Added_Mass.time, simData_50F_4.MODEL.Added_Mass.mass, label="Added Mass")
ax2.plot(simData_50F_5.MODEL.Hourglass_Energy.time, simData_50F_5.MODEL.Hourglass_Energy.energy)#, label="Hourglass Energy")
ax2.plot(simData_50F_5.MODEL.Internal_Energy.time, simData_50F_5.MODEL.Internal_Energy.energy)#, label="Internal Energy")
ax2.plot(simData_50F_5.MODEL.Kinetic_Energy.time, simData_50F_5.MODEL.Kinetic_Energy.energy)#, label="Kinetic Energy")
ax2.plot(simData_50F_5.MODEL.Total_Energy.time, simData_50F_5.MODEL.Total_Energy.energy)#, label="Total Energy")
ax2.plot(simData_50F_5.MODEL.Added_Mass.time, simData_50F_5.MODEL.Added_Mass.mass)#, label="Added Mass")
ax3.plot(simData_50F_6.MODEL.Hourglass_Energy.time, simData_50F_6.MODEL.Hourglass_Energy.energy)#, label="Hourglass Energy")
ax3.plot(simData_50F_6.MODEL.Internal_Energy.time, simData_50F_6.MODEL.Internal_Energy.energy)#, label="Internal Energy")
ax3.plot(simData_50F_6.MODEL.Kinetic_Energy.time, simData_50F_6.MODEL.Kinetic_Energy.energy)#, label="Kinetic Energy")
ax3.plot(simData_50F_6.MODEL.Total_Energy.time, simData_50F_6.MODEL.Total_Energy.energy)#, label="Total Energy")
ax3.plot(simData_50F_6.MODEL.Added_Mass.time, simData_50F_6.MODEL.Added_Mass.mass)#, label="Added Mass")
ax4.plot(simData_50F_7.MODEL.Hourglass_Energy.time, simData_50F_7.MODEL.Hourglass_Energy.energy)#, label="Hourglass Energy")
ax4.plot(simData_50F_7.MODEL.Internal_Energy.time, simData_50F_7.MODEL.Internal_Energy.energy)#, label="Internal Energy")
ax4.plot(simData_50F_7.MODEL.Kinetic_Energy.time, simData_50F_7.MODEL.Kinetic_Energy.energy)#, label="Kinetic Energy")
ax4.plot(simData_50F_7.MODEL.Total_Energy.time, simData_50F_7.MODEL.Total_Energy.energy)#, label="Total Energy")
ax4.plot(simData_50F_7.MODEL.Added_Mass.time, simData_50F_7.MODEL.Added_Mass.mass)#, label="Added Mass")
ax5.plot(simData_50M_4.MODEL.Hourglass_Energy.time, simData_50M_4.MODEL.Hourglass_Energy.energy)#, label="Hourglass Energy")
ax5.plot(simData_50M_4.MODEL.Internal_Energy.time, simData_50M_4.MODEL.Internal_Energy.energy)#, label="Internal Energy")
ax5.plot(simData_50M_4.MODEL.Kinetic_Energy.time, simData_50M_4.MODEL.Kinetic_Energy.energy)#, label="Kinetic Energy")
ax5.plot(simData_50M_4.MODEL.Total_Energy.time, simData_50M_4.MODEL.Total_Energy.energy)#, label="Total Energy")
ax5.plot(simData_50M_4.MODEL.Added_Mass.time, simData_50M_4.MODEL.Added_Mass.mass)#, label="Added Mass")
ax6.plot(simData_50M_5.MODEL.Hourglass_Energy.time, simData_50M_5.MODEL.Hourglass_Energy.energy)#, label="Hourglass Energy")
ax6.plot(simData_50M_5.MODEL.Internal_Energy.time, simData_50M_5.MODEL.Internal_Energy.energy)#, label="Internal Energy")
ax6.plot(simData_50M_5.MODEL.Kinetic_Energy.time, simData_50M_5.MODEL.Kinetic_Energy.energy)#, label="Kinetic Energy")
ax6.plot(simData_50M_5.MODEL.Total_Energy.time, simData_50M_5.MODEL.Total_Energy.energy)#, label="Total Energy")
ax6.plot(simData_50M_5.MODEL.Added_Mass.time, simData_50M_5.MODEL.Added_Mass.mass)#, label="Added Mass")
ax7.plot(simData_50M_6.MODEL.Hourglass_Energy.time, simData_50M_6.MODEL.Hourglass_Energy.energy)#, label="Hourglass Energy")
ax7.plot(simData_50M_6.MODEL.Internal_Energy.time, simData_50M_6.MODEL.Internal_Energy.energy)#, label="Internal Energy")
ax7.plot(simData_50M_6.MODEL.Kinetic_Energy.time, simData_50M_6.MODEL.Kinetic_Energy.energy)#, label="Kinetic Energy")
ax7.plot(simData_50M_6.MODEL.Total_Energy.time, simData_50M_6.MODEL.Total_Energy.energy)#, label="Total Energy")
ax7.plot(simData_50M_6.MODEL.Added_Mass.time, simData_50M_6.MODEL.Added_Mass.mass)#, label="Added Mass")
ax8.plot(simData_50M_7.MODEL.Hourglass_Energy.time, simData_50M_7.MODEL.Hourglass_Energy.energy)#, label="Hourglass Energy")
ax8.plot(simData_50M_7.MODEL.Internal_Energy.time, simData_50M_7.MODEL.Internal_Energy.energy)#, label="Internal Energy")
ax8.plot(simData_50M_7.MODEL.Kinetic_Energy.time, simData_50M_7.MODEL.Kinetic_Energy.energy)#, label="Kinetic Energy")
ax8.plot(simData_50M_7.MODEL.Total_Energy.time, simData_50M_7.MODEL.Total_Energy.energy)#, label="Total Energy")
ax8.plot(simData_50M_7.MODEL.Added_Mass.time, simData_50M_7.MODEL.Added_Mass.mass)#, label="Added Mass")
fig_energy.tight_layout(pad=0.1)
fig_energy.legend(loc='center left', bbox_to_anchor=(1, 0.5));
Force-deflection response#
The figure below compares the predicted force-deflection response to the PMHS results.
Show code cell source
fig_cont, (ax1) = plt.subplots(1,1, figsize=(10,10))
fig_cont.suptitle('Single ribs anterior-posterior bending', y=1.02)
ax1.set_ylabel('Force (N)')
ax1.set(xlabel='Normalized displacement [-]')
F50rib4x=np.array(simData_50F_4.Boundary.conditions_Posterior_force.time) *ureg.parse_expression('milliseconds')
F50rib4y=np.array(simData_50F_4.Boundary.conditions_Posterior_force.force) *ureg.parse_expression('kilonewton')
F50rib4_filtered = SF.cfc(F50rib4x.reshape(-1,1), F50rib4y.reshape(-1,1)* 1, 60, False)
F50rib5x=np.array(simData_50F_5.Boundary.conditions_Posterior_force.time) *ureg.parse_expression('milliseconds')
F50rib5y=np.array(simData_50F_5.Boundary.conditions_Posterior_force.force) *ureg.parse_expression('kilonewton')
F50rib5_filtered = SF.cfc(F50rib5x.reshape(-1,1), F50rib5y.reshape(-1,1)* 1, 60, False)
F50rib6x=np.array(simData_50F_6.Boundary.conditions_Posterior_force.time) *ureg.parse_expression('milliseconds')
F50rib6y=np.array(simData_50F_6.Boundary.conditions_Posterior_force.force) *ureg.parse_expression('kilonewton')
F50rib6_filtered = SF.cfc(F50rib6x.reshape(-1,1), F50rib6y.reshape(-1,1)* 1, 60, False)
F50rib7x=np.array(simData_50F_7.Boundary.conditions_Posterior_force.time) *ureg.parse_expression('milliseconds')
F50rib7y=np.array(simData_50F_7.Boundary.conditions_Posterior_force.force) *ureg.parse_expression('kilonewton')
F50rib7_filtered = SF.cfc(F50rib7x.reshape(-1,1), F50rib7y.reshape(-1,1)* 1, 60, False)
M50rib4x=np.array(simData_50M_4.Boundary.conditions_Posterior_force.time) *ureg.parse_expression('milliseconds')
M50rib4y=np.array(simData_50M_4.Boundary.conditions_Posterior_force.force) *ureg.parse_expression('kilonewton')
M50rib4_filtered = SF.cfc(M50rib4x.reshape(-1,1), M50rib4y.reshape(-1,1)* 1, 60, False)
M50rib5x=np.array(simData_50M_5.Boundary.conditions_Posterior_force.time) *ureg.parse_expression('milliseconds')
M50rib5y=np.array(simData_50M_5.Boundary.conditions_Posterior_force.force) *ureg.parse_expression('kilonewton')
M50rib5_filtered = SF.cfc(M50rib5x.reshape(-1,1), M50rib5y.reshape(-1,1)* 1, 60, False)
M50rib6x=np.array(simData_50M_6.Boundary.conditions_Posterior_force.time) *ureg.parse_expression('milliseconds')
M50rib6y=np.array(simData_50M_6.Boundary.conditions_Posterior_force.force) *ureg.parse_expression('kilonewton')
M50rib6_filtered = SF.cfc(M50rib6x.reshape(-1,1), M50rib6y.reshape(-1,1)* 1, 60, False)
M50rib7x=np.array(simData_50M_7.Boundary.conditions_Posterior_force.time) *ureg.parse_expression('milliseconds')
M50rib7y=np.array(simData_50M_7.Boundary.conditions_Posterior_force.force) *ureg.parse_expression('kilonewton')
M50rib7_filtered = SF.cfc(M50rib7x.reshape(-1,1), M50rib7y.reshape(-1,1)* 1, 60, False)
ax1.plot(experiment['Norm-disp-50M'], experiment['Force-50M'], label = "Experiment female middle aged", **plotExperimentMale)
ax1.plot(experiment['Norm-disp-50F'], experiment['Force-50F'], label = "Experiment male middle aged", **plotExperimentFemale)
ax1.plot(0.629*simData_50F_4.Boundary.conditions_Anterior_displacement.displacement, -1000*F50rib4_filtered, label = "VIVA+ 50F Rib4", linestyle='-', linewidth=3, **plot50F)
ax1.plot(0.582*simData_50F_5.Boundary.conditions_Anterior_displacement.displacement, -1000*F50rib5_filtered, label = "VIVA+ 50F Rib5", linestyle='--', linewidth=3, **plot50F)
ax1.plot(0.549*simData_50F_6.Boundary.conditions_Anterior_displacement.displacement, -1000*F50rib6_filtered, label = "VIVA+ 50F Rib6", linestyle=':', linewidth=3, **plot50F)
ax1.plot(0.518*simData_50F_7.Boundary.conditions_Anterior_displacement.displacement, -1000*F50rib7_filtered, label = "VIVA+ 50F Rib7", linestyle='-.', linewidth=3, **plot50F)
ax1.plot(0.568*simData_50M_4.Boundary.conditions_Anterior_displacement.displacement, -1000*M50rib4_filtered, label = "VIVA+ 50M Rib4", linestyle='-', linewidth=3, **plot50M)
ax1.plot(0.518*simData_50M_5.Boundary.conditions_Anterior_displacement.displacement, -1000*M50rib5_filtered, label = "VIVA+ 50M Rib5", linestyle='--', linewidth=3, **plot50M)
ax1.plot(0.483*simData_50M_6.Boundary.conditions_Anterior_displacement.displacement, -1000*M50rib6_filtered, label = "VIVA+ 50M Rib6", linestyle=':', linewidth=3, **plot50M)
ax1.plot(0.461*simData_50M_7.Boundary.conditions_Anterior_displacement.displacement, -1000*M50rib7_filtered, label = "VIVA+ 50M Rib7", linestyle='-.', linewidth=3, **plot50M)
ax1.set_ylim(0,180)
ax1.set_xlim(0,50)
fig_cont.tight_layout(pad=0.5)
fig_cont.legend(loc='center left', bbox_to_anchor=(1, 0.5));
plt.show();
c:\Users\iraeus\Anaconda3\envs\viva\lib\site-packages\dynasaur\calc\cfc.py:66: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.
sampled_array = np.array(sampled_array)
c:\Users\iraeus\Anaconda3\envs\viva\lib\site-packages\numpy\core\_asarray.py:83: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.
return array(a, dtype, copy=False, order=order)
c:\Users\iraeus\Anaconda3\envs\viva\lib\site-packages\numpy\core\_asarray.py:83: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.
return array(a, dtype, copy=False, order=order)
c:\Users\iraeus\Anaconda3\envs\viva\lib\site-packages\numpy\core\_asarray.py:83: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.
return array(a, dtype, copy=False, order=order)
c:\Users\iraeus\Anaconda3\envs\viva\lib\site-packages\numpy\core\_asarray.py:83: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.
return array(a, dtype, copy=False, order=order)
c:\Users\iraeus\Anaconda3\envs\viva\lib\site-packages\numpy\core\_asarray.py:83: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.
return array(a, dtype, copy=False, order=order)
c:\Users\iraeus\Anaconda3\envs\viva\lib\site-packages\numpy\core\_asarray.py:83: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.
return array(a, dtype, copy=False, order=order)
c:\Users\iraeus\Anaconda3\envs\viva\lib\site-packages\numpy\core\_asarray.py:83: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.
return array(a, dtype, copy=False, order=order)
c:\Users\iraeus\Anaconda3\envs\viva\lib\site-packages\numpy\core\_asarray.py:83: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.
return array(a, dtype, copy=False, order=order)
