Use of Dynamic FEA for Design Modiﬁcation and Energy Analysis of a Variable Sti ﬀ ness Prosthetic Foot

Featured Application: Variable sti ﬀ ness in the force response of a prosthetic foot. Abstract: Di ﬀ erent tasks and conditions in gait call for di ﬀ erent sti ﬀ ness of prosthetic foot devices. The following work presents a case study on design modiﬁcations of a prosthetic foot, aimed at variable sti ﬀ ness of the device. The objective is a proof-of-concept, achieved by simulating the modiﬁcations using ﬁnite element modeling. Design changes include the addition of a controlled damping element, connected both in parallel and series to a system of springs. The aim is to change the sti ﬀ ness of the device under dynamic loading, by applying a high damping constant, approaching force coupling for the given boundary conditions. The dynamic modelling simulates mechanical test methods used to measure load response in full roll-over of prosthetic feet. Activation of the element during loading of the foot justiﬁes the damped e ﬀ ect. As damping is in contrast to the main design objectives of energy return in prosthetic feet, it is considered important to quantify the dissipated energy in such an element. Our design case shows that the introduction of a damping element, with a high damping constant, can increase the overall rotational sti ﬀ ness of the device by 50%. Given a large enough damping coe ﬃ cient, the energy dissipation in the active element is about 20% of maximum strain energy.


Introduction
Energy storage and return (ESR) prosthetic feet serve as designated spring systems, storing energy in mid-stance of gait, that again is released for propulsion of the foot in late stance [1,2]. Increased energy return has been reported to benefit users with increased body propulsion and decreased sound limb loading [3]. Deformation of the foot allows roll-over throughout stance, imitating the response of a biological ankle. Stiffness of the foot contributes to this roll-over, affecting the forces and moments exerted on the residual limb, and users' whole body kinematics [4]. The flexible, energy-storing parts of ESR feet are generally leaf springs, fabricated from fiber reinforced polymer materials. Stiffness curve and damping response is fixed and depends on material properties, shape, and thickness [5,6]. Feet are prescribed to patients in stiffness categories, primarily based on weight and activity levels [7,8].
It has been indicated that ESR foot performance during diverse ambulatory tasks and conditions could be improved further by varying the stiffness of the prosthetic device [4,9,10]. Hence, adjustable or adaptable stiffness could be beneficial to users, especially in tasks that differ from the walking motion for which the foot is optimized. Recent publications demonstrate an increased research interest in prosthetic feet that have adjustable stiffness and the effect of prosthetic stiffness on user in different tasks. Prosthetic feet with modular ankle stiffness have been proposed by Shepard et al. [11] and Glanzer et al. [12]. These devices alter stiffness by shifting a support in a beam deflection system, allowing a change in forefoot stiffness between steps. Shepard et al. [13] reported that users could, on average, feel a stiffness change of approximately 8% and were consistent when choosing preferred stiffness. Adamczyk et al. [4] showed the effect of forefoot and hindfoot stiffness on walking mechanics at different speeds on a prosthetic with adjustable component stiffness.
The application considered in this work is the prosthetic foot, Pro-Flex®(Össur, Reykjavík, Iceland). The current design has a fixed stiffness curve and the modifications considered are aimed at enabling controlled change of stiffness in the device. The objective is twofold. First, to construct a model to simulate the dynamic force reaction of a prosthetic foot throughout the gait cycle. The purpose here is to track shifts in function as design is modified. Secondly, to induce change in stiffness with a controlled damping element, using the model to verify the change in function of the modified design. Abrupt changes in the force response of a prosthetic device during loading will affect the foot roll-over shape. We stipulate that change in force response of a prosthetic device under load, needs to be gradual or damped, to ensure continuity in the roll-over. Therefore, our methodology is to induce stiffness change, by introducing high damping coefficients in a system of springs connected in parallel and series.
Damping effects in the early stance of prosthetic feet have been credited with increased self-selected walking speed of amputees [14,15]. Damping characteristics of prosthetic feet have also been shown increased loading rates on uninvolved side in unilateral amputees [16]. In general, the concept of damping contradicts the design objectives of energy return. Therefore, it is considered important to quantify the energy dissipation that damped force coupling introduces. The force transfer by elements with high damping coefficients are used to simulate force transfer that could be described as a combination of viscose and frictional damping in a functional force coupling. The assumption is made that the functional element has a low off-state damping coefficient and a high on-state damping coefficient.
The modeling is done using finite element (FE) analysis (FEA). Previously published work on the application of FEA in prosthetic foot design has focused on investigating structural characteristics of components, rather than function. These works show FEA of individual mechanical components for stress analysis to models of the whole prosthetic, focusing on the force response in the anisotropic carbon fiber material. Omnasta et al. [17] studied the load bearing response of a prosthetic foot, from spring blade and foot cover deflection, to stresses in foot-pylon connection. Several studies also utilize force values from gait analysis data [17][18][19][20]. The use of FEA for direct design purposes is also common, seen in [21], where different curvatures of a spring blade prosthetic foot are investigated, and [22], where shapes of a solid plastic foot (Niagara foot) are evaluated. In addition, there has been interest in FE modeling of the interface between residual leg and prosthetic socket since the early 90s [23] and is ongoing [24].
FEA that focuses on the dynamic modeling and damping response of ESR prosthetic feet are most notably modal and frequency analysis of sprint feet. Vinney et al. [25] evaluated the suitability of FEA as a method to study dynamic characteristics of prosthetic feet. Noroozi et al. [26] studied the damped elastic response to impulse synchronization of running feet, where it is suggested that user performance can be enhanced by optimizing the prosthetics design to this parameter. A recent study from Rigney et al. [27] uses FEA for mechanical characterization of several types of sprinting feet, including damping under dynamic loading. The use of FEA has also proven useful in more complex studies of the prosthetic gait for design purposes. Mahmoodi et al. [19] used FEA to optimize design parameters of three types of prosthetic feet to estimate optimal stiffness based on gait analysis data, such as ground reaction forces and roll-over shape. Bonnet et al. [18] used kinetic data from gait analysis combined with FEA to calculate stress, strain, and energy storage in a prosthetic foot. The following work presents a FE model, intended to serve as a tool to predict shifts in dynamic function, in reaction to change in mechanical design.
In this study, we build a dynamic FE model of a prosthetic foot and validate it against mechanical test results to support the application described above. The model is used to simulate a standardized mechanical test [28], where the foot performs a full roll-over under a dynamic load profile. The construction of the model is meant to serve as a platform for iteratively modifying designs. The analysis is transient to enable damping in the model. A new joint, restrained with a spring and damper element, is defined in the modified model. This will result in a more compliant force response of the foot, as shown in previous work [29]. The resulting rotational stiffness of the modified design, for varying damping coefficient in the element is then compared to the original model. Furthermore, the strain energy is analyzed, in order to quantify the dissipation of using damped force coupling for change in stiffness.

Preparation of FEA
The prosthetic foot used for the study is a Pro-Flex®(Össur, Reykjavík, Iceland) (size 27, Category 5), a proven energy efficient ESR foot [3] with a high ankle range of motion (RoM) [30]. The model is made using ANSYS Workbench®(Canonsburg, PA, USA). The model is constructed in three steps and each step is validated against mechanical testing to confirm assumptions on material properties and contact behavior. Transient structural analysis is carried out for the 3D geometry model of the foot, where a mechanical testing standard [28] is simulated to allow for real time function of the foot in a continuous step cycle.
The elastic parts of the foot are three carbon fiber (CF) composite leaf springs; top blade, middle blade, and sole blade; see Figure 1 These are defined as flexible bodies, while other components are assumed fully rigid. The blade thickness is retrieved from the layup recipes of the CF composite material and implemented in the geometry model. The blades are meshed with a sweep method and modeled with solid/shell type elements (ANSYS; SOLSH190). These are solid type elements that are free from locking in bending dominated conditions and specially designed to model thin to moderately thick plates and shells. The element type allows for definition of layers with different material properties through the thickness of the body [31].
The CF composite blades are unidirectional, with thin layers of bidirectional woven material on top and bottom. The layer thickness, material properties, and fiber angle are assigned to each layer of the blades. The tensile and compressive moduli in the fiber direction of the unidirectional CF composite dominates the force response of the blades in bending, although the shear modulus is also important in a high deformation condition [32]. To get an accurate estimate of the moduli, a three-point bend test was performed on 5 mm thick, unidirectional CF composite samples. The flexural modulus was found to be 97 GPa and set as the Young's modulus for the material in fiber direction. Other values on material properties were taken from ANSYS pre-defined "Epoxy Carbon (230 GPa) Prepreg". Rigid parts are connected by joint elements and assigned a rotational degree of freedom, imitating the actual hinges and ball joints of the prosthetic foot; see Figure 1. Bonded contact is defined for the connection between the three blades and the connection of blades to the clamps. Frictionless contacts are defined at the wedging of the blades in bending and frictional contact between sole blade and tilt table; see Figure 2. Friction factor is adjusted so that no slip occurs. Contact formulation and normal stiffness settings are adjusted to imitate the physical conditions of each case.

Model Validation
The model is validated in three steps. First, the in-line production tests of mechanical uniaxial loading of the blades, as cantilevered beams, are simulated and compared to measured values. This comparison confirms assumptions made on material properties. Second, simulations are done of static toe and heel deflection testing of prosthetic feet and compared to measured values. Herein, the sample foot is rigidly attached to a uniaxial test machine and aligned at a fixed angle, −15 • for the heel and 20 • for the toe, against a free rolling plate. The machine compresses the foot for a set displacement and the resulting force is measured. In this process, the model assumptions on contact and joint element settings between different parts are adjusted. Lastly, the simulation of the roll-over task of the technical specification ISO/TS 16955 [28] is done and compared to measurements for validation of the model. Figure 2 shows the schematic structure of the test equipment and procedure. One investigated cycle consists of the dynamic stance phase of a single step. The foot sample is attached to a rod and force is applied to a ball joint at the other end. The rod and attached foot sample are thus free to rotate around the ball joint. The M-shaped force curve represents the ground reaction force of heel-to-toe walking; see Figure 2b). The maximum force peak is set to 824 N, representing the typical ground reaction force of a 70 kg person walking at normal speed. The sample is pushed vertically down on a rotating table (tilt table). The rotation is constant displacement controlled according to a synchronized profile; see Figure 2c). The heel contact at the start of the cycle is at a −20 • decline, rolling over to 40 • incline and pushing the foot up at the toe. Control functions for force and rotation are defined as per the technical specification [28]. A load cell positioned at 500 mm above the foot (about the height of the anatomical knee joint) on the rod measures the resulting forces and moments. A sequence of the roll-over is shown in Figure 3b). The stance phase time is 0.600 s, corresponding to a 1 s normalized step time.
Every-day use ESR prosthetic feet are used with a cosmetic foot shell, made of soft elastomers. The foot shell was excluded in the FE modeling as the visco-elastic member would increase the complexity and reduce reliability of the model. In contrast to the standard [28], the comparison mechanical measurements were also executed without a foot cover to allow for direct comparison.
The heel-to-toe profile of the standard does not fully represent normal walking, as it exaggerates the angle at initial contact and loading response of stance phase, resulting in higher moments and anterior-posterior forces than would be expected from clinical gait analysis results. Results of the roll-over simulations are compared to the mechanical testing and presented in a rotational stiffness curve; see Figure 4. The angular moment in the sagittal plane is calculated for the ankle system in the rotational joint between the main body and top blade. This moment reaction is plotted against the ankle angle, characterizing the rotational stiffness of the ankle system throughout stance phase. In the mechanical testing, the force and moment reaction in the load cell are used to calculate the moment in the ankle joint. The ankle angle is measured between the sole blade and the pylon with a 2D video analyzing software, TEMA (Image Systems, Linköping, Sweden). The curve starts at a 0 • angle with no moment reaction, then through plantarflexion, back through the zero point at mid-stance, through positive values in dorsiflexion, and then returns to the starting position.
The difference between the model and measured values is most noticeable at maximum plantarflexion, where the measured rotational stiffness is 5.7 Nm/ • but the model results are 6.5 Nm/ • (14% higher). The results are sensitive to sample alignment and this difference is believed to be due to the shift of foot sample in the dynamic measurements. The coefficient of determination of measured values against the model results was calculated as 0.987. As the measurement data is retrieved at 1000 Hz, the data points are not evenly distributed throughout the curve in Figure 4. The primary objective of the modeling is to track the effect of design changes. Overall, this difference between measured values and simulation is considered within requirements for a valid comparison and valuation of shifts in results as the model is modified. The main source of convergence difficulties in the FE model is related to the instability of the solution in periods of the cycle where abrupt contact changes occur. This is apparent as the heel rolls over the tilt table in the beginning of load response and as the forefoot comes in contact with the tilt table just before mid-stance. A three-point moving average filter was applied to smooth the results before further calculations are made, where abnormal fluctuations in the simulation results occur.

Strain Energy
From the FE analysis, we can retrieve the reacting forces and moments in all parts of the prosthetic foot. Figure 5 shows the total strain energy for the entire foot, as well as, for each individual blade. The results show how each blade is strained and unloaded over stance phase. During early stance, the energy is stored in the middle blade and sole blade. However, during late stance when only the toe is loaded, the strain energy is primarily in the top blade and middle blade. Over the stance phase, the middle blade has the greatest effect on the overall stiffness of the device. The work at the boundaries is calculated from the force and displacement of the pylon, summed with product of torque and rotation on the tilt table at the sole. Total work at the boundary will equal the system energy. If we assume that frictional losses in moving parts are negligible, the boundary work will equal the strain in the spring blades and losses due to material damping. These damping losses are seen as hysteresis in the angular stiffness curve from mechanical testing; see Figure 4. The model does not account for material damping, causing the boundary work to be equal to the strain energy curve over the cycle time.

Simulating Design Modification
The first objective of design modification in the model is to produce a change in the characteristic stiffness curve. The modelled foot is stiffer than recommended for the weight that the load profile assumes. The intention is to make changes in the model that make the foot more compliant and then increase stiffness with a controlled element, preferably to the original stiffness. The aim is to change the stiffness with an element that is purely damping and evaluate the damping factor needed for an efficient force transfer. This could be frictional damping or otherwise damped force coupling. The assumption made is that abrupt force coupling would have adverse effects on the user as discontinuity in the foot roll-over. The simulation assumes an active element that is actuated from mid-stance to late stance in order to affect the forefoot stiffness. The timing of the actuation is also varied to study the effect of actuation under load. The second objective is to quantify the effect of the energy dissipation as a result from the damping on the energy return.
The model is initially made more compliant by introducing a new joint in the foot design. Instead of a rigid mechanical link connecting the middle blade to the main body (see Figure 1), this component is modeled as two parts, connected by a translational joint. Allowing movement over the new joint makes the system unsupported, so a linear spring element is defined over the joint. The spring stiffness is set at 500 N/mm after iterative simulations, aiming for a 30% overall increase in RoM for the given boundary conditions. The introduction of a new spring element that acts both in series and parallel with the existing spring system lowers the overall spring stiffness of the system.
The resulting rotational stiffness of the compliant model of the foot is shown in Figure 6, and is compared to simulation of the original prosthetic foot design. The rotational stiffness curve shows a slight decrease in moment and 25.9% increase in RoM on the hindfoot (plantarflexion). The effect is more on the forefoot, as RoM is increased from 18.9 to 26.3 • and the resulting ankle stiffness is decreased by 42% (6.9 to 4.0 Nm/ • ). A damping element is defined with the assumption that the controlled element has an off-state damping effect, represented with a low damping coefficient of 1 Ns/mm. The effect of this damping is negligible on the hindfoot, but a small hysteresis effect appears in the forefoot reaction curve.

Active Functional Joint
For the actuated element, a damping element is included in the simulation and is defined over the translational joint, parallel to the spring, previously defined in the mechanical link; see Figure 7. The element is inactive at the start of the cycle and activated at mid-stance (0.3 s). It is then turned off shortly before the end of the cycle, at 0.56 s (93% of stance), allowing for the spring energy to be released. The resulting rotational stiffness curves for four values of damping coefficients for the active element are shown in Figure 6, compared to the FE simulation of the original prosthetic foot and the simulation of the compliant model with no active element. The damping coefficients are set to 25, 50, 100, and 300 N s/mm. The rotational stiffness approaches the original rigid link model as the damping coefficient is increased. Rotational stiffness evaluated at maximum moment in dorsiflexion increases up to 51% when applying 300 Ns/mm damping compared to the compliant setup (from 4.0 to 6.1 Nm/ • ). In Figure 6, the curves for actuated damping are cut off at peak ankle angle to clarify the plot, excluding the return curve showing the hysteresis due to damping. The energy dissipation is shown and quantified in more detail in Figure 8. To track the effect of the energy dissipation in the actuated element, the energy and work is calculated for model components. Figure 8 shows the stored and dissipated energy in the model for the simulation of the actuated damping element for three damping coefficients. The lowest damping coefficient (25 Ns/mm) is disregarded, as it induces the least change in stiffness. The plot shows the recoverable energy and damping energy for the three cases. The recoverable energy is the combined strain energy from the strain in the CF blades and the spring energy from the spring element defined over the new joint. The strain in the CF blades increases with higher damping coefficients. The difference between the combined strain energy curves is primarily the higher spring energy, as a result of increased displacement of spring element in the transient joint for lower damping coefficients. Figure 8. Dissipated damping energy (grey) and recoverable strain energy (black) development in the foot over the stance phase for different damping coefficients used in the actuated element, strain energy being from CF blade strain, and spring element strain. The element is actuated at mid-stance (0.3 s) and turned off at the end of the stance, as the element starts to contract after elongation. The damping energy from the actuated damping model rises at mid-stance, as the element is activated, showing the total dissipated energy at the end of the cycle. Dissipated energy for the cycle increases when increasing the damping coefficient from 50 to 100 Ns/mm but decreases when damping coefficient is further increased to 300 Ns/mm. This confirms a maximum dissipation for the given system and boundary conditions at a damping coefficient between 100 and 300 Ns/mm. The dissipated energy reduces the strain energy available for foot propulsion in late stance phase, decreasing the efficiency of the foot as an energy returning spring system. The sum of the strain energy and damping dissipation equals the work done on the system at the boundaries. As the simulations are investigated under constant displacement conditions, the dissipated energy adds to the boundary work rather than decrease the recoverable strain energy. The efficiency can be represented by a ratio between recoverable energy and work at the boundaries. This ratio between peak values at roughly 80% of the stance phase show that efficiency is in line with the absolute value of dissipated energy, or 77.6%, 75.6%, and 79.6% for 50, 100, and 300 Ns/mm damping, respectively.
To further study the effect of the actuated damping, simulations are made where the actuation time for the damping is delayed into the terminal stance. Figure 9 shows the results for the 300 Ns/mm damping coefficient for two additional actuation times. As before, the curves are compared to the simulation of the Pro-Flex and the compliant model. In addition to actuation at mid-stance (0.3 s) the element is actuated at 0.35 and 0.4 s (58% and 67% of stance phase, respectively). These results show that even as the element is activated during considerable loading of the foot, at 67% of stance phase, the moment versus angle curve is continuous. This supports the hypothesis that dampened actuation of change in load support will result in a gradual stiffness change that will not interrupt the roll-over of the foot. As focus in prosthetic design is increasingly shifted towards actuated elements and microprocessor-controlled function, the need arises for a robust method to assist with identifying the best parameters for both sensing and control. These results also demonstrate the effectiveness of using FEA in such tasks.

Conclusions
An FE model of a prosthetic foot was presented and validated against actual measured values on a physical prosthetic foot. The model can be used to predict mechanical reactions to forces, moments, and displacements. It also proved to be a useful tool to analyze the functional characteristics of the prosthetic foot. Furthermore, the model can be used to simulate the response of the device when modifications are implemented in comparison to the original design.
The work showed how the stiffness of a system of series and parallel springs can be altered by introducing a pure damping element. Stiffness change of 51% was demonstrated. Although the change achieved is less than demonstrated in previous work [11,12] it is considerably higher than the reported user perception of 8% [13]. Furthermore, for the given system and boundary conditions, it was presented that the dissipation of a damping element will reach a maximum at a given value of damping coefficient. For higher values, the energy dissipation decreases due to decreased velocity, causing the element to be better described as force coupling, becoming more efficient at transferring the forces as the coefficient is higher. The dissipation in the damping element has been shown to reduce the recoverable energy by over 20% for the utilized load profile. The benefit of being able to change or adapt the foot stiffness during loading, needs to outweigh the effect of this lost energy.
For the load condition of constant displacement, change in the system function such as the increased damping in the case presented, will have more impact on work on the boundary than how the foot deforms. The strain energy does therefore not vary as much as might be expected in gait analysis of a user walking on the foot, where it is expected that the reaction would be more complex. It is anticipated that users will adjust their gait to the force response of the foot, affecting both deformations, as well as, reactive forces and moments.
As stated previously, the force and rotation input functions for mechanical testing [28] differ from what is seen in clinical gait analysis when investigating amputee gait at normal walking speed. This is apparent when the strain energy curves are compared to earlier work. The strain energy results presented here, compares quite well with results for fast walking speed presented by Bonnet et al. [18], while it differs for normal walking speed in the same work. The load in early stance phase is higher than for normal walking speed and the mid-stance minimum occurs later in the stance phase. Further work will include simulations with changed boundary conditions. Force vector data from gait analysis will be used to dictate vertical force curves and tilt table rotations in order to better imitate normal walking. Furthermore, data will be gathered to produce different gait patterns, such as ascent and descent of ramps and stairs.