Numerical assessment of directional energy performance for 3D printed midsole structures

1.   Introduction

Energy is the property to perform work on an object. Although energy is a scalar quantity, the forces exerted by a surface as energy are vector quantities having both magnitude and direction. Consequently, the terminology of directional energy refers to the energy associated with a directional movement ^[1]. When the force towards a direction is multiplied by the amount of deformation, the result is the energy in that direction. The concept of directional energy transfer (DET) is to find a structure that can transfer energy from one direction to another, which was first developed for sport shoes by Fuss ^[1] in 2009. DET in sport shoes maximizes the energy return in the forward direction and thus increases athletic performance in terms of running speed. Nowadays, the state-of-art additive manufacturing has opened the door to fabricate complex structures for many applications (Figure 1), which allows DET to be applied to other products, like airless tires. DET obtained through the tires would help a car move forward, increasing energy efficiency and reducing fuel consumption. It could be possible to generate structures with functional increments to enhance DET. Consequently, a quantitative way to measure the efficiency of transfer is needed for the exploration of different structures. Without loss of generality, the concept of DET is explained using the shoe midsole since most previous works were mainly related to that. However, we do not consider DET as an application-specific concept.

A shoe sole contains three layers: insole, midsole and outsole ^[2] (Figure 1a). Among all these layers, the shoe midsole serves the shoe sole's primary purpose, and traditionally it provides cushioning ^[3]. While an athlete is walking or running, forces applied on the ground by the foot ^[4] can be decomposed into a downward and backward force. According to Newton's third law, reaction forces are exerted by the ground, respectively, to the upward direction helping to stand and to the forward direction pushing the athlete forward. When the athlete contacts the sport surface, there is some work done by the athlete on the sport surface ^[5]. Energy is transferred from the athlete to the surface through the foot and the shoe ^[6]. An athlete's performance can be highly regulated by the interaction between the foot and sports surface ^[7,8]. Athletic footwear has been perceived as a mechanism by which the running economy can be improved ^[9]. Performance enhancement is also a primary motivating reason that runners try new footwear ^[10]. The selection of appropriate footwear is often advocated as an essential requirement for distance running ^[11]. Running shoes with greater shoe cushioning, greater longitudinal shoe stiffness, and greater shoe comfort were associated with improved running economy ^[12].

The sports surface deforms due to work done by the athlete and stores the potential energy. The concept of DET is to utilize the potential energy to move forward more efficiently. In other words, a DET structure in this context is a shoe midsole that can transfer some of the 'vertical' energy into a form of horizontal force and deformation. A quantitative way to measure the efficiency of transfer is thus needed for the exploration of different structures. The purpose of this paper is to develop a general mathematical model for the assessment of DET, which can be used for design optimization and iteration for energy performance. The contributions are:

1). The deficiencies in the current mathematical model of DET using spring elements are analyzed. The concepts and components that are necessary to study DET are identified and remarked.

2). A new formulation is developed to directly calculate the energy in a particular direction, and it is concluded that there is no DET using static equilibrium analyses.

3). The formulation is further extended to dynamic finite element analysis using frame elements, so that the mathematical derivation and calculation are unified and general even for more complicated structures.

Numerical study and analysis are done on different structures and configurations to show that it is possible to realize DET with the new formulation, and the performance of DET is affected by both the geometry and the material of structures.

The paper is organized as follows. Related works are reviewed in Section 2. Section 3 explains the fundamentals of DET, investigates a case, and gives our observations. According to our observations, a new formulation is developed for DET in Section 4, and further generalized in Section 4.2. The developed analyses are done on a few examples, and the results are presented in Section 5. The paper is concluded in Section 6.

2.   Related works

Uses of more complex geometry to modify mechanical properties are becoming more and more common in advanced manufacturing technology ^[13]. In the automotive sector, a kind of application is producing airless tire (Figure 1b). The purpose was to modify the lateral stiffness while maintaining proper contact patch. With the help of three-dimensional (3D) printing technologies, it is possible to design complex structures ^[13,14] to enhance DET performance. Lattice structures have many superior properties, such as a lightweight structure due to its high specific stiffness and strength ^[15]. Mechanical performances of lattice structures depend on various factors such as the cell topology, number of cells, geometric parameters (e.g., strut diameter and cell size), material and manufacturing process, and structural boundary and loading conditions ^[16]. The effect of unit cell size on the elastic modulus, shear modulus, and Poisson's ratio of triangulated lattice structures shows that the elastic modulus and the shear modulus decreased as the cell size increased ^[17]. In a previous study ^[18], three lattice structures corresponding to different loading modes were designed and tested showing some structures have higher energy absorption efficiency than the others. Hence, it is possible to control the energy return by using different lattice structures to have directional energy properties. By modifying density, angle, and anisotropy, it can control both the "softness" and "bending" of a sole ^[19]. A method is presented to precisely control the contact forces and pressure over large contact areas between the foot and a deformable shoe ^[20]. Auxetic materials behave unconventionally under deformation enhancing material properties such as resistance to indentation and energy absorption showing potential to improve sporting protective equipment ^[21]. The functional flexibility of the structures motivates this paper to develop the measurement of DET, so that it can be used as an objective to optimize the structures.

From previous works, the concept of energy return is implemented in sports mainly. Energy return is a similar concept as DET, and it refers to the return of the stored energy ^[22], but it does not consider the direction of the energy. It is believed that sport shoes can play an essential role in a runner's performance ^[23] since running efficiency depends on the interaction between the foot and the sport surface. Wearing a running shoe is changing the geometry of the foot-ground interface and the stiffness of the interface due to the deformation of the midsole ^[24]. Ground reaction forces and kinematic variables vary with shoe hardness and shoe geometry ^[25]. Improving forefoot push-off facilitates the augmentation of forward acceleration and ultimately enhances athletic performances ^[26,27].

A reason behind the scope of energy return being limited is that the returned energy is not being efficiently used. The efficacy of this energy return concept, according to Nigg et al. relies on the energy returned at the right location, at the right time, with the right frequency ^[28]. The concept of energy return led to the development of DET as a critical measure of shoe efficiency. Only a few works have been conducted in the DET area previously. DET can be optimized if the bounce tubes of Adidas shoes are rotated, and the overall stiffness changed (by altering the tube length) ^[29]. Several parameters related to DET energy transferred (as a percentage of vertical input energy), energy returned in the horizontal direction, and the total system energy ^[30]. Dickson's thesis ^[31] showed how energy transfer could be determined for classified systems, and a midsole can be designed capable of transmitting energy from a vertical direction to a horizontal direction. Experimental results were used in the formulation to determine the amount of energy transfer. Although it was a big step in this research area, only the spring systems were studied, and a few assumptions were needed to make the mathematical model work. More studies are needed to enhance DET's fundamental understanding if it is applied to more general cases like 3D-printed lattice structures.

3.   Directional energy transfer (DET)

This section explains the concepts and formulations of DET based on the works of Dickson ^[31] and his supervisor–Fuss ^[1]. Without loss of generality, the concept of DET is illustrated with the midsole structure shown in Figure 2, which is taken from the 1 degree-of-freedom (DOF) prototype designed by Dickson ^[31]. As the boundary conditions are similar along the contact surface, the structure can be generated by repeating a unit cell throughout the midsole. Hence, the analysis can be done on a single unit cell. The study assumes that the foot is in contact with the shoe sole in the whole forefoot strike process. The deformation of a structure happens for the applied load at the contact point between the foot and the sole. The shape and geometry of the unit element are considered during the formulation. It is also assumed that the contact point does not change while the load changes over a short period ($\sim0.1s$). Energy loss due to viscosity is not considered in the mathematical model.

While running, an athlete strikes the ground with one's forefoot and applies a load to the shoe midsole. The midsole deforms due to the load and stores energy as a form of load and deformation (Figure 2c). There are some deformations in the $x$- and $y$-directions, and the structure deforms until it reaches an equilibrium. When the athlete takes off, the deformed structure tends to go back to its initial state and releases the stored energy as a force acting onto the foot. Dickson ^[31] simplified the structure as a torsion spring at the bottom, which is also fixed in translation. The structure is inclined at an angle of $\theta$. When a load ($F_a$) is applied, the spring deforms and stores energy. In the context of DET, the energy in the deformation is separated in $x$ and $y$:

where $E_a^x$ and $E_a^y$ are the energies caused by the applied forces $F_a^x$ and $F_a^y$ in $x$- and $y$-directions, respectively. This is the work done by the applied force causing the deformation. Dickson ^[31] saw that as the energy absorbed by the structure due to the loading.

The next step is to find out how the structure returns the energy, and it is based on how the energy is stored in the structure. In a spring system, energy is stored in the twisted springs. For example, in this system with only one torsion spring, the potential energy stored in the spring $E_s$ is calculated by its spring constant $k$ and deformation $\Delta\theta$, which is also the total energy $E_{tot}$ in the system:

The spring energy ($E_{s}$) also needs to be separated into the $x$- and $y$-directions: $E_s^x$ & $E_s^y$. Dickson ^[31] achieved this separation by decomposing the spring force $F_s$ into the $x$- and $y$-directions – $F_s^x$ and $F_s^y$, and integrated them with the deformations, resulting in

It is assumed that there is no energy converted to other forms besides deformation, and thus the energy is conserved, i.e.,

The amount of transfer can be measured by the gain of energy from the loading to the stored amount in the $x$-direction, or the loss of energy in the $y$-direction:

In this example, since the load was applied only in the $y$-direction, $F_a^x$ and thus $E_a^x$ was zero, i.e., $E_{trans} = E_s^x$. Therefore, since there is deformation in $x$, $E_s^x$ should be non-zero, and Dickson concluded that DET was observed in this 1-DOF structure.

3.1.   Mechanism with two linear springs

The DET principle was also applied to a 2-DOF mechanism with two linear springs ^[31], and it is briefly reviewed here to demonstrate the model further. The mechanism is shown in Figure 3a. The two springs are inclined at $90^\circ$, and the central node is constrained in that inclined plane. The contact point is at the top, and the bottom end (A) has a fixed boundary condition. The mathematical formulas are the same as the ones of 1-DOF, except for the spring forces, which need to be calculated by both springs:

where

$F_{s1}$ and $F_{s2}$ are the forces of the two springs ($1$ & $2$) obtained by their deformations $\Delta l_1$ and $\Delta l_2$. Similar to the 1-DOF case with load in the $y$ direction, the structure deforms in both the $x$- and $y$-directions, and thus it was concluded that DET was observed.

Experimental analysis was also performed on the structure to confirm there is DET, and the corresponding prototype of the mechanism was fabricated as in Figure 3c. This figure is taken from Dickson's thesis ^[31], and it overlaps two pictures of the prototype before and after the load is applied to show the structure is indeed moved in the $x$-direction. The parts in white are made of rigid materials, and the comparatively transparent block between them are made of elastic material. The block allows compression from both sides: parallel and perpendicular to the interface of the two rigid parts, simulating the two-linear-spring mechanism. More details can be found in Dickson's thesis ^[31]. The prototype was compressed to a suitable maximum. The vertical motion caused the prototype to move horizontally, and the structure moved back to the original position after the upper platen returned to its original location. Since the $x$-direction forces were not measured, the energy transfer was calculated based on the $y$-direction. The load profiles for the vertical component against the deflection in $y$ are shown in Figure 3d. Three different loads were tested corresponding to the three curves in the chart. The energy in the $y$-direction was obtained by calculating the area under the curves, where the top part is the loading curve, and the bottom part was the unloading curve. The area enclosed by this loop gave a measurement of the energy loss in the $y$-direction, which could reflect the energy gain in the $x$-direction as in Eq (3.4). As energy loss in $y$ was observed from the curves, it was concluded that the 2-DOF mechanism is capable of transferring energy from the $y$-direction to the $x$-direction, and thus it has DET.

3.2.   Observations and remarks

Although Dickson ^[31] provides a foundation for the study of DET, the study was based on simple spring systems and relied on several assumptions that might not be met. This resulted in some inaccurate conclusions and prohibited the model from being applicable to general cases. Here, we make a few observations and remarks on the formulation, which will help develop a complete numerical analysis.

Firstly, the DET was confirmed mainly because of the case that the applied force in $x$ was zero and thus the energy $E_a^x$ from Eq (3.1) inputted in the system was zero, so any stored energy $E_s^x$ in $x$ was transferred from $y$. It should be noted that the work causing a change in energy is the work done by the net force, rather than by an individual force. Therefore, Eq (3.1) is valid only if the force is applied to a free body (see Figure 4a), i.e., there is no other force. When it is a spring system like the one shown in Figure 4b, the net force should be the difference between the applied force ($F_a$) and the spring force ($F_s$). Work is the multiply of a force and a displacement in a direction, and it can be either positive or negative. If the force is in the same direction as the displacement, the force is doing positive work; otherwise, it does negative work. Even there is no applied force in the $x$-direction, the spring force does a negative work to have a movement in $x$. The positive work in the unloading phase compensates the negative work in the loading phase, instead of coming from $y$. Disregarding the negative work will overestimate the amount of DET in the system, and thus we have the following remark.

Remark 1. To find the energy inputted into a system in a particular direction due to loading, it should integrate the net force (sum of applied and internal forces) in that direction.

Secondly, although the physical experiment showed that there was energy loss in the $y$-direction, the loss could be due to heat or others; and more importantly, it was not sure that even there is energy gain in the $x$-direction, it really helps moving forward. Indeed, the energies in Eq (3.4) are separated in different axes, but they are still not directional. For an illustration, assume an athlete is running in the right direction. When a force $F$ to the right is applied, the work done by it to move a distance $\Delta x$ to the right is $E = F\Delta x$, and this energy is in the same direction as the running direction. However, if there is a force equal in magnitude but towards the left (i.e., $-F$) applied to move the same distance to the left ($-\Delta x$), the work done is the same: $E = F\Delta x$, but it is against the running. Therefore, it is not enough to only separate the energies in different axes, but they should also have a forward or backward classification, and we have the following remark.

Remark 2. The work needs to be analyzed on the axis of interest and also in a desired direction (e.g., with position and negative values on the axis).

One may argue that there is no energy loss in the mathematical model as it is not modeled, and the system is at equilibrium, so the energy gain in one direction must be transferred from another. This brings to our third and last observation.

Thirdly, the conservation of energy was based on the assumption of equilibrium state, and actually, the whole DET study was using static analysis. If a system is at equilibrium, all the forces are balanced, and thus the sum of forces (i.e., the net force) in any direction is zero. In other words, the applied forces should equal to the spring forces at equilibrium, i.e., $F_a^x = F_s^x$ and $F_a^y = F_s^y$. If that is the case, referring to the original formulations, the energies inputted to one direction (Eq (3.1)) should always equal to the energies stored in the same direction (Eq (3.3)), and the energy transfer $E_{trans}$ must be zero all the time (Eq (3.4)). This contradicts the origination of the mathematical model. Even if the net force is considered in the calculation of work, the net force at equilibrium is zero, and thus the work is also zero. Then, $E_{trans}$ is just related to the spring energy stored in the system, which is not meaningful.

Although the static equilibrium analysis is simple and can ease the energy calculation, assuming every instant during loading and unloading is at equilibrium is not realistic. This is because besides deformation, the energy is in multiple forms like acceleration and velocity, and the whole motion is dynamic. To truly compute the DET, we have the following remark.

Remark 3. The study of DET should consider the dynamics, i.e., the instantaneous changes in both applied forces and internal forces, of the whole loading cycle.

Based on the observations and remarks we have made, the next section presents our methodology to develop the formulations for DET.

4.   Methodology

To address the Remarks 1 and 2, we reformulate the mathematical model of DET to properly calculate the directional energies. To address the Remark 3 and to generalize the formulation for more complex problems than spring systems, we apply the dynamic finite element method to the computation of DEP. The details are presented in the following.

4.1.   Reformulation of DET

Although the energies in Eqs (3.1) and (3.3) are calculated separately for different axes, they are still scalar quantities in the axes. This is because when the net force causes a deformation in the same direction, the work done by the net force is always positive, no matter if it is going to the left or to the right. The potential energy stored in a spring is always positive too, since the deformation is squared, e.g., Eq (3.2). To distinguish the forward and backward direction, we borrow the concept of negative work by considering the desired direction. Instead of integrating the net force with the amount of deformation, we assign the integral a direction concerning the rest state of the deformation. In other words, with a given desired direction, the limits on the definite integral are always set by the direction.

It might be easier to understand by a graph shown in Figure 5. Assume the direction of interest is towards the right. There is a linear spring that deforms along the $x$-axis by an applied force $F$ to the right. The cycle starts with applying the load at $x_1$ and stretching the spring. The spring force increases linearly with its elongation, and thus the net force decrease linearly. When the spring force is the same as the applied force, and the net force equals zero, the system reaches its equilibrium at $x_2$, i.e., $F = kx_2$. After that, the external load is removed, and the spring is being restored by the spring force at $x_3$. With the decrease in the displacement, the force also decreases in magnitude until it becomes zero at $x_4$. In this illustration, $x_1 = x_4 = 0$ and $x_2 = x_3 = F/k$. In the common definition, the work done by the force in moving an object from point $a$ to point $b$ is given by

and thus the works done by the net force for the loading and restoring phases are defined as

In this case, both will give the same result:

but the works are done in a different direction.

In fact, it is evident on the graph that the areas under the curves are separated by the $x$-axis. Since we want the energy that helps moving forward (i.e., $+x$), the area above the $x$-axis should be positive, and the area below should be negative. To reflect this in the equation of work, the integration should be in accordance with the forward direction. That is, we set the integral to be always from a smaller limit to a larger limit (because the right direction is the positive side along the $x$-axis). Therefore, for $x_4 < x_3$, the work in the restoring phase should be defined as

This gives $W_1+W_2 = 0$, meaning no energy gain or performance in the $x$-direction. This example is the simplest case with a linear spring and a constant force applied; mass and acceleration are also not considered. In a general situation, the material and the structure can be non-linear, and the applied force can be changing over time, even during restoring phases. Fortunately, this modeling is general as it does not make any assumption on the loading conditions, and it only requires the calculation of net force and deformation.

To sum up, the energy on an axis of interest can be found by the areas between the curves of net force and the axis. To measure the directional performance, instead of comparing the energies between the loading and stored ones (the result is still a scalar), we compute the difference between the areas above and below the axis so that the energy is signed with respect to the desired direction. The mathematical representation can be unified into one form by integrating the net force over different intervals $\{\xi\}$ of deformation – the curve needs to be separated into sub-regions since there should not be a $x$-value with two different $F$-values in an integral. Assume the desired direction is $+x$, the work done by a net force $F_n(x)$ over the whole cycle is defined as

It is worthwhile to remark that Dickson ^[31] introduced DET as the energy loss in the $y$-direction being transferred to the $x$-direction. Our study actually finds it difficult to confirm that the energy gain in the $x$-direction is coming only from the $y$-direction. Therefore, our new formulation here no longer considers this loss-and-gain relationship, but directly the energy-change in the direction of interest for the whole deformation cycle. The re-formulation can have the applied loads in any direction and calculate the energy performance in any direction. We used the $+x$ direction as an illustration as it is the direction of interest in the shoe midsole application, but it can be any direction. In the following, we are going to use a new term 'Directional Energy Performance' (DEP) to make clear of this difference. DEP is defined as the net amount of energy in a specific direction after a motion cycle. Hence, more energy would be obtained in the desired direction from a better performing structure.

4.2.   Dynamic modeling

When a system and its conditions are simple (e.g., linear properties and constant loads), the study of energy does not need to worry about what happens in between, because the initial and final conditions often have all the information required. However, when varying forces and accelerations are considered, it needs to look at the action at every instant in time. Following the Remark 3, DEP is a dynamic process, and thus we extend the model using the dynamic finite element analysis (FEA). The element type is also extended from spring to frame. A frame element can model a straight bar of an arbitrary cross-section, which can deform in the axial and perpendicular direction to the rowspan="i"xis of the bar. A frame is capable of carrying both axial and transverse forces, as well as moments (see Figure 6).

A structural system that considers masses, accelerations, and dynamic forces is modeled by the following expression:

where ${\bf{M}}$ is a mass matrix, ${\bf{C}}$ is a damping matrix, ${\bf{K}}$ is a stiffness matrix, $\dot{u}$ and $\ddot{u}$ are time derivatives of the deformation $u$. For time-history analysis for the derivatives, we used the direct integration method, so that no transformation of the equations into different forms is carried out. A popular implicit method for numerical integration called 'Newmark Method' is applied here. This method is unconditionally stable and has no restriction on the size of the time step. Although damping is highly related to vibration in reality, it requires a set of experiments to determine damping ratios of the structure at two separate frequencies, so the damping effect is not considered here for the sake of simplicity.

Let the frame have a uniform cross-sectional area $A$ with a length $L$ and the material be linear and have an elastic modulus $E$. The mass matrix ${\bf{M}}$ and the stiffness matrix ${\bf{K}}$ are expressed as

For a frame element inclined at an angle $\theta$ with the $x$-axis, a transformation matrix is required to include the orientation of the component, and the equation becomes:

with

Given an applied load $F_a = \{F_{1x}, F_{1y}, M_1\}$, the deformations in the axes of interest (e.g., $\Delta x_1$) at different time instants during the motion can be computed by Eq (4.5). The internal forces ($F_s$) at the loading point need to be calculated to find the net force. Since the internal forces only depend on the structure's deformation state, the static part of the equation can be applied with a given deformation ($u$), i.e., $F_s = {\bf{K}}u$. After that, all the information required to measure the DEP of a structure is obtained, and the results can be analyzed.

Note that, although Eq (4.5) is just for one single element, the matrix can be easily extended when there are more elements, as of the finite element method (FEM). The mechanism with two inclined frames is also tested and reported in the result section. In addition, we applied the frame elements over spring elements to develop a more general method (the spring system is a simplified representation of the actual design considering only the axial loads), but other types of elements can be employed as well to deal with other domains. In fact, the stiffness of a particular structure could be complex, which would make the discrete analysis more complicated. Since the FEM represents the small blocks of materials and connects them to other elements, they can be shaped to almost any geometry subjected to any loading irrespective of complexity.

5.   Analyses and results

With the new mathematical model developed for DEP, this section presents the numerical analyses to test how and when the DEP property is obtained. Following the previous study ^[31], we also focus on the shoe midsole structure. The load applied was related to a person's weight, but it should be distributed uniformly on the contact surface, so a fraction of the load is used in our tests, i.e., $\frac{1}{10}$th. To verify that DEP is a dynamic process and it can only be obtained via dynamic analysis, we first perform a static analysis using the new formulation on the same 2-DOF mechanism in Figure 3. A constant load is applied, and both spring and the frame elements are tested. Secondly, a dynamic load (i.e., vary with time) is used to test and compare a quasi-static analysis and the dynamic analysis. Thirdly, to confirm the energy is transferred from one direction to another, we apply various loading cases to analyze the differences. Lastly, a case study on lattice structures is conducted to show how a structure can be selected for the performance of DEP. In this paper, as we are only studying directional energy on unit cells and there are only a few elements in a cell, the computation time is short for all examples---at most 1.5 s.

5.1.   Static analysis on spring and frame elements

In this section, the structure with two inclined linear springs shown in Figure 3 is re-visited. The two springs are attached to each other at right angles, and the bottom one is inclined at an angle $\theta = 60 °$ to the horizontal. The spring constants are $k_1 = 3333$ N/mm, $k_2 = 6666$ N/mm and their lengths are $l_1 = l_2 = 30$ mm. A downward, vertical load $F = -82.215$ N is applied at the contact point. Linear springs can only carry axial forces, so at equilibrium, the spring forces in Eq (3.5) should be balanced with the corresponding components of the applied force:

which gives $\Delta l_1 = -0.02136$ mm and $\Delta l_2 = -0.00617$ mm, as in compression. The deformations in the $x$-direction can then be calculated:

It results in $\Delta x_1 = -0.01068$ mm, $\Delta x_2 = 0.00534$ mm, and the total $x$-deformation is $\Delta x = \Delta x_1 + \Delta x_2 = -0.00534$ mm. The net spring force in the $x$-direction at equilibrium is

Note that although there is a net $x$-deformation, the net $x$-spring-force is zero, which is not a surprise because the system is at equilibrium. To apply Eq (4.1), the whole cycle is integrated over two intervals: one is from $0$ to $\Delta x$, the other is from $\Delta x$ to $0$. Similar to the analysis in Section 4.1, this gives

This integration should be with respect to $x$, $x_1$, and $x_2$, where $x = x_1 + x_2$. However, since the net spring force ($F_s^x$) is zero, and the applied force ($F_a^x$) is also zero. Therefore, $E^x = 0$ J without the need for further derivation, meaning that no energy performance in $x$ is possible.

Because spring elements can only carry axial forces, the applied force needs to be converted to find the corresponding deformations and then the integration, making the method case-specific. We apply the FEM to generalize the model, and frame elements that can carry arbitrary forces are used here as a demonstration. For the frame element, the spring force can be replaced by internal force extracted using the stiffness matrix. Physical properties of frame elements such as elasticity, cross-sectional area, inclination angle, and length are used in FEM to determine the deformation at the shoe-sole and foot contact point. As in the previous case of the spring element, a similar loading condition is applied. Let the frame elements have a cross-sectional area of $A = 50$ mm² and elastic moduli $E_1 = 2000$ N/mm² and $E_2 = 4000$ N/mm². Assembling the element stiffness matrices (e.g., Eq 4.4), a system of linear equations can be formed as ${\bf{F}} = {\bf{\bar{K}}}u$, where ${\bf{F}}$ is the force vector, ${\bf{\bar{K}}}$ is the global stiffness matrix, and $u$ is the deformation vector. After applying the boundary conditions, the force vector is $F = [ 0 \; 0 \; 0 \; 0 \; \mathrm{-82.215} \; 0 ]^T$, and in turns the deformation vector is computed as $u = [ 0.4133 \; 0.2140 \; \mathrm{-0.0465} \; 2.0184 \; \mathrm{-0.7251} \; \mathrm{-0.0697} ]^T$. Therefore, the directional deformations of interest at the contact point is $\Delta{x_2} = 2.0184$ mm. Similar to the spring system, integration over $x$-intervals gives the energy values in $x$-direction:

The result is similar to the spring system that $E^x = 0$ J, but the calculation here is generic and can be applied to any orientation or combination of forces and structures. The same analysis has been done on other structures too using static analysis, and the results also show no DEP. The static analysis is based on the equilibrium condition, and thus the force matrix shows the internal forces at the very moment when the structure reaches equilibrium. Hence, the applied forces are always the same as the reaction forces. Therefore, we can conclude that they do not show any energy performance for any linear geometry with linear material properties while using static analysis.

5.2.   Quasi-static vs. dynamic analyses

As mentioned in the Remark 3, the DEP is a dynamic process: the applied forces are changing, and the system is not at equilibrium. Therefore, it is crucial to consider the dynamics so that the DEP can be measured. Firstly, the dynamic load used to model the forefoot strike ^[1] is applied here, as shown in Figure 7(a). One cycle of the strike lasts $0.14$ s starting from the initial contact between the foot and the midsole to a point when the foot is off the ground and no longer applies forces on the midsole. Secondly, we compare the quasi-static and dynamic analyses results because both can also take dynamic loads. The 1-DOF structure inclined at an angle of $60 °$ shown in Figure 2 is employed in this comparison, but with a frame element.

In the quasi-static analysis, the load is applied and updated so slowly that the structure also deforms very slowly, meaning that it considers additional static forces to model the dynamic phenomenon. In other words, it separates the dynamic load into many small time steps so that the acceleration of applying the load is small enough to be neglected. Thus a static analysis is done in each time step. Ideally, each increment should be infinitesimal, but then the change in each step would be too small to visualize, so a larger step size is used here for a better illustration of the results. By sampling the dynamic load at different time instants, a sequence of static loads can be obtained: $F_a(t_0), F_a(t_1), F_a(t_2), \dots, F_a(t_n)$. When it is advanced from step $i$ to $i+1$, the net force is increased by $F_a(t_{i+1})-F_a(t_i)$ and then comes back to zero as the internal force increases with the deformation until equilibrium. The force-deformation curve is shown in Figure 7(b), where the triangles on the left of $y$-axis are formed by the net force in $y$ – $F_n^y$, and the triangles on the right is formed by the net force in $x$ – $F_n^x$. That is, the contact point deforms to the right and the down direction, so the deformation is positive in $x$ and negative in $y$.

Each triangle in the curve represents one step and is the result of the static analysis like in Figure 5. The energy $E^x$ computed with this configuration is $1.804$ J. Note that if the step size is infinitely small, the triangles will vanish so that the curve will become a flat line, and the energy $E^x$ will be zero. Therefore, although applying the quasi-static analysis can consider the dynamic load, the core is still the static analysis concerning the structure's equilibrium state, which cannot capture any energy performance.

In the dynamic analysis, the applied force is updated faithfully according to the time, and the system may not reach an equilibrium state. Besides, the inertia forces are also considered so that even the applied force has changed its direction. The structure may still be deforming in the original direction. These result in the internal force different from the applied force, and thus a non-zero net force could be obtained. The result is plotted in Figure 7(c). Similarly, the $F_n^y$-curve is on the left, and the $F_n^x$-curve is on the right of the $y$-axis. The $F_n^x$-curve tells how the energy is in the $x$-direction, which is the direction of interest. The areas above and below the $x$-axis are apparently different, and the net energy in the $x$-direction is found $E^x = 12.246$ J. The net $x$-energy is positive, so there is energy contributing to the forward direction, which is an encouraging result.

One may observe that the area above the $x$-axis being larger than the area below, meaning that the positive net energy comes from when the structure approaches its deformed state. In fact, there are two ways to obtain a positive net energy when the structure deforms and returns: one is that the energy is negative during deformation but positive and larger in return, like bouncing on a compressed spring and using both the reaction and spring force to jump higher; the other is that the return has negative energy, but the deformation has positive and larger energy, like slipping on a banana peel but managed not to fall. This example is a kind of the second situation, i.e., the energy inputted to the forward direction is more than the energy returned to the backward direction. Another observation is that the curve does not reach the starting position at the end of the forefoot strike loading cycle. This means that the structure does not return to its rest state, but it does not imply plastic deformation. Recall that a strike cycle only lasts $0.14$ s, and after that, the foot and shoe leave the ground until the next cycle commences when the foot and shoe regain contact with the ground. There is not sufficient time in a strike cycle for the structure to reach an equilibrium, but there is plenty of time between the cycles for the structure returning to its rest state.

5.3.   Comparing loading cases

A positive net energy in the $x$-direction is observed in the previous test case, but it is not sufficient to conclude that there is DEP. This is because the net change may only be due to the loss in energy generated by the work of $F_a^x$ but not necessarily be transferred from the energy generated by the work of $F_a^y$. Therefore, this section tests the model on some designed loading cases to verify the transfer of energy. Specifically, the ones with only the $y$-loads are applied. This is because if the applied load is working only in the vertical direction and there is energy generated in the horizontal direction, there must be some work in the $x$-direction that is transferred from the $y$-direction. Two additional loading profiles are: variable $y$-loading and constant $y$-loading. The variable $y$-loading is the same as the one of a forefoot strike but without the $x$-component, and the constant $y$-loading is a fixed $y$-load applied in the first half of the cycle and then withdrawn in the second half (Figure 8(a)). They are tested on the same 1-DOF frame structure the same as in the Section 5.2.

The results are plotted in Figure 8(b), (c), each of which includes the curves of applied forces ($F_a^x, F_a^y$), internal forces ($F_s^x, F_s^y$), and net forces ($F_n^x, F_n^y$). For the variable $y$-loading profile, the curves look similar to the one in Figure 7(c), but since $F_a^x$ is zero, the curves of $F_s^x$ and $F_n^x$ are the same. Similarly, the net energy in the $x$-direction can be computed by the areas under that $F_n^x$-curve above and below the $x$-axis, and $E^x = 15.583$ J, which is even higher than the one with non-zero $F_a^x$. From this result, we can confirm that the applied forces in $x$ did not contribute to the net energy positively (in fact, negatively), and thus the energy transferred from the $y$-direction to the $x$-direction in this configuration is actually $E_{trans} = 15.583$ J.

For the constant $y$-loading profile, the curves of net forces look very similar to the one in Figure 5, which is from a static analysis. Indeed, applying a constant load to a dynamic analysis is basically a static analysis. However, due to the acceleration at the instants that the force is applied and removed, the net energy in $x$ is not zero but very small, i.e., $E^x = 2.164$ J. It would be zero if the constant force was applied steadily and slowly. This result verifies the Remark 3 that DEP is a dynamic process, and the loading profile also matters to the DEP performance.

5.4.   Case study: lattice structures

The study here is motivated by the lattice structures used in 3D-printed parts. One of them is shown in Figure 9(a). However, to have DEP, the structure should not be symmetric so that deformation in the $x$-direction is possible for a load applied in the $y$-direction. In this case study, the fore-foot strike loading is applied. Therefore, we modify the lattice structure and obtain three different structures as shown in Figure 9, including a 1-element bar inclined to the right side (inspired by Adidas Airblade Shoe), a 2-element spring (inspired by 3D printed airless tire), and a 4-element structure (stiffen spring). Let it be noted that since a lattice structure can have complex topology resulting the force not uniformly distributed, all the forces throughout the lattice must be considered if we calculate the potential energy in it. However, since we are interested in the athlete who are stepping on the structure rather than the structure itself, we are looking at the foot and calculate how it does work on the structure during the forefoot strike. As work is the energy transferred to or from an object via the application of force along a displacement, the integral of force at the bottom of foot is computed along the trajectory of displacement, and as explained in Section 3.2 and Figure 4, this force should be the net force at the contact point. The force applied by the foot is given, and the FEA done on the lattice is only to obtain the opposite force and the trajectory of the contact point.

To analyze and compare, we used the same material volume for different structures within the same domain, which is a unit cell with both height and width $30$ mm. This considers the average dimension of the shoe sole so that the structure can reach the midsole height and can be repetitively used along the midsole length and width. The objective is to observe any DEP and if the structure is strong enough to stand the loads. Two different materials (Acrylonitrile butadiene styrene – ABS and aluminum) and two separate volumes ($3000$ and $6000$ mm³) are tested to demonstrate the performance. The frame elements are assumed to have a circular cross-section area. The maximum stress was found by comparing the stresses for all elements, and the element subjected to the highest stress was located. The obtained maximum normal stress was compared with the yield strength $\sigma_{max}$, and the maximum shear stress was compared with the shear strength $\tau_{max}$ of the material. Shoe midsole should not undergo plastic deformation in each strike cycle, and hence we considered the yield strength of the material as a benchmark for this analysis.

The net forces against the deformation for ABS with $3000$ mm³ volume are plotted in Figure 10, which shows that the three structures indeed have different performances. The data are summarized in the Table 1. The 2-element structure shows a higher energy performance than the 1-element, whereas the 4-element structure has the lowest energy performance for any material irrespective of volume. In terms of the normal and shear stresses, the 2-element and the 4-element both have higher stress values than the 1-element. Since the 2-element and the 4-element have shallower cross-section area, a higher bending stress is experienced. When a smaller volume ($3000$ mm³) is used, the material with a smaller yield stress (e.g., ABS) causes the material's yield or failure more easily. This eliminates the 2-element structure (although it has a higher DEP), making the 1-element a better choice when less volume is needed. When a larger volume of material ($6000$ mm³) is used, less stress concentration is observed for the 2-element, i.e., normal stress is $-29.790$ N/mm². In this case, the 2-element is more preferable for two reasons. First, the structure is remaining below the yield strength point. Second, the 2-element structure has a higher energy performance. Therefore, if a material has a higher elastic modulus and strength is used, the 2-element structure is preferable over the 1-element.

6.   Conclusions

In this paper, the mathematical model is modified to assess directional energy performance (DEP), proposing significant concepts lacked in preceding research. The model incorporates recognizing net forces uniting spring force in a distinct axis of interest. An analytical method to distinguish the behavior of DEP for any generalized structure is defined. Spring analysis and static FEM analysis confirm their ineffectiveness due to estimating the structure's equilibrium phase. The study verifies DEP is a dynamic process by executing static, quasi-static, and dynamic analysis for various loading circumstances. Energy performance is verified as long as the process is dynamic yet for simple structures with linear materials. This mathematical model serves as a handy computational tool for designing structures to control energy conversion or harvest energy. We compared several predefined case studies and compared results among themselves to find better DEP showing design.

However, our investigation has some limitations. For instance, we did not acknowledge damping, and we used specific orientation of elements that satisfy the boundary conditions. We did not consider the non-linear behavior effect of geometry change on the stiffness matrix. Also, we only considered the circular cross-section and applied the load based on a forefoot strike. Hence, the results may vary depending on the structure's physical properties, cross-section area, or higher load, causing failure of the structure. In the future, more analyses will be conducted considering optimum structure design and material selection in different applications. Further analysis is needed with more complex structure, non-linear materials, and/or dynamic structure to explore DEP's possibility. The strain-state-based plasticity model ^[32] will be used to speed up the computation. In addition, physical experiments will be conducted to validate the DEP concept.

Acknowledgments

This paper acknowledges the support of the Natural Sciences & Engineering Research Council of Canada (NSERC) grant #RGPIN-2017-06707.

Conflict of interest

The authors declare there is no conflict of interest.