# Finite Element Analysis of Shock Absorption of Porous Soles Established by Grasshopper and UG Secondary Development

Authors: Xiaoying Liu; Yong Yue; Xuyang Wu; Yanhua Hao; Yong Lu
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

On the basis of computer aided modeling technology, this paper proposes a porous structure modeling method based on Grasshopper visual programming language and Unigraphics NX (UG) secondary development platform. The finite element model of the foot was established, and then models of shoe soles with four basic porous structures of cross, diamond, star, andx were established. Each structure was set with a cylindrical radius of 1, 2, and 3 mm, and a total of 12 porous structure sole models were established. The shock absorption effect of the sole on the foot was evaluated by the deformation of the sole, the peak plantar pressure, and the peak stress of metatarsal bones. It is found that the maximum value of the sole deformation of the diamond porous sole is 4.725 mm, the peak plantar pressure is 105.1 Pa, and the first and second metatarsal peak pressures are 2.230 MPa and 3.407 MPa, which have the best shock absorption effect. It shows that the porous structure plays an important role in the cushioning of the sole. The biomechanical effects of porous soles on feet are studied by computer-aided technology and finite element analysis. This study provides a new research method for the cushioning design of shoe soles and has important reference value for the design of footwear.


---

## Body

## 1. Introduction

The foot is a natural shock-absorbing structure [1–3]. The arch of the foot plays a very important role in cushioning vibration of the foot [4]. In the pressure state of walking, running, jumping, etc., the arch will be deformed. Through the deformation, the impact force of the human body on the ground is reduced, and the human tissue is protected from damage [5]. The arch of the foot and the muscles attached to it, like a shock absorber, protect the human foot. Heel fat pads are natural shock-absorbing tissues of the human body [6, 7]. The plantar skin is dense, with a thicker cuticle, and the subcutaneous tissue is strong and has a fat pad, which can buffer the pressure of the sole on the nerves and blood vessels [8, 9]. At the same time, the fat pad is filled with scattered small and strong connective tissue fiber bundles, which can fix the skin and limit the displacement of fat. Their arrangement direction is compatible with the main pressure direction [10]. The plantar skin has the functions of wear resistance, pressure resistance, and weight bearing, which is conducive to walking and load-bearing stability. Over the last 40 years, the popularity and the practice of running have considerably increased worldwide [11, 12]. There is great interest in developing new products and technologies, such as improvement of running shoes. The materials and structures of the sole could influence the biomechanics, dynamics [13–15], etc.In terms of the prevention of foot diseases, footwear is of great importance. In addition to protecting the foot from external injury and beautifying the body appearance, consumers are increasingly paying attention to the wearing comfort of shoes [16]. As a basic function of footwear, the shock-absorbing performance of sole is one of the important factors affecting the comfort of footwear [17]. Porous materials refer to a type of hole-combination which is composed of solid corrugated edges or walls. It is widely used in our life or engineering applications. For example, honeycomb, foam, cork, bone, and even bread and steamed bread are porous. Porous structure materials have good characteristics of cushioning, sound absorption, heat insulation, low density, and high specific rigidity, and it can be used as heat insulation, packaging, structure, and floating materials. Because of the excellent physical properties of porous structures, their performance and modeling methods have been extensively and deeply studied [18]. At present, there are a variety of methods for preparing porous structures, such as familiar hexagonal honeycomb structure materials. However, only a few of these methods can achieve batch production. Most of the methods have complicated processes, and it is difficult to effectively control the porous structure during the preparation process so that the obtained porous structures still have some performance deficiencies [19]. With the development and application of 3D printing technology, a new approach has been brought to the preparation of porous structures. Computer modeling can be used to accurately model, and the porous framework and macroscopic pores can be designed according to the needs. For example, common research works in related fields include scaffold design in tissue engineering [20–22] and internal structure design in mechanical property optimization [23, 24]. Due to the increasing computer power, the complexity of the model is growing. Generally, the more complex the models are, the larger becomes the uncertainty in the model outputs due to randomness in the input parameters. Vu-Bac et al. [25] provided a sensitivity analysis toolbox consisting of a set of MATLAB functions that offer utilities for quantifying the influence of uncertain input parameters on uncertain model outputs, and it was proved that the method was useful and effective.General sneaker sole is made up of the outsole, midsole, and insole. The outsole is the part that comes into contact with the ground. The middle bottom is sandwiched between the outsole and the inner bottom. The top of the inner bottom comes into contact with the sole of the foot. The outsole is mainly used to start and stop sliding and wear resistance. The inner sole is mainly used to support the foot surface and disperse the plantar pressure. The middle sole has the effect of vibration reduction and buffering. Researchers have long studied the cushioning structure of soles, Nike’s Nike Air, Li-Ning’s Li-Ning Bow, Honeycomb shock absorption, mechanical device, and more. At present, there are few studies on the application of porous structure to shoe soles. In this paper, a method of porous structure modeling based on Grasshopper and UG secondary development is proposed, and it is applied to the midsole of shoes, which is of great convenience and operability. The porous modeling parameters are simple, with only diameter as a parameter, so the probabilistic sensitivity analysis is not carried out.In recent years, the plantar pressure is measured by experiments, and the cushioned effect of the sole is analyzed by comprehensive subjective evaluation. However, it is very difficult to obtain the internal stress of the foot directly by experimental method. The finite element model can effectively simulate the motion state of the foot and obtain the stress distribution rule of the internal bone and soft tissue. In addition, the finite element method is easy to simulate and analyze, and the reasonable finite element model of the foot can well simulate the biomechanical behavior of the foot bone and soft tissue, providing a strong method and theoretical support for further study of the biomechanics of the human foot skeletal and musculoskeletal system. Hu et al. [26–28] built a complete finite element model of the foot. Yu et al. [29] established the finite element model of female foot and high-heel sole structure to explore the effect of high-heel shoes on foot biomechanics. Cheung et al. [30] calculated the influence of plantar fascia stretching and the distribution of stress and strain by finite element simulation. Li et al. [31] simulated the jumping moment with the finite element model and analyzed the stress distribution of the posterior bone of the foot. It was found that the stress was mostly concentrated in the posterior talus joint and the back of the talus. In conclusion, the finite element method can greatly reduce the workload of human foot research and can simulate the foot movement well and provide the foundation for the research of foot biomechanics.In this study, a porous structure modeling method was proposed and applied it to the midsole structure of shoes. Combining with the finite element modeling method, the biomechanical analysis of foot-shoes is carried out, and a reasonable evaluation method of shoe damping performance is formulated, which is of great reference value to the design of shoes and the customization of orthopedic shoes.


## 2. Materials and Methods

This study consists of two phases, each with distinct goals. The first phase is to propose a modeling method of porous structure based on Grasshopper visual programming language and UG secondary development platform. Specifically, a complex parameterized 3D wireframe model was established by Grasshopper, and the 3D wireframe was expanded into a solid model through UG NX open solid modeling technology to complete the modeling of a porous structure. The second phase is the finite element analysis of the sole model with different porous structures. Specifically, to improve the shock absorption and cushioning performance of sports shoes, porous structures with different pores and different lattice types were constructed and applied to the heel area of the middle sole of sports shoes, and then ABAQUS was used to perform finite element analysis on the soles with porous structures of different pores and different lattice types, and the plantar pressure on the surface of these soles was compared to obtain the pore and lattice type of the porous structure with the best vibration absorption effect on the soles.


### 2.1. Porous Modeling

In this section, Grasshopper and UG are used to establish the porous structure. Figure1 is the process flow chart.Figure 1
A flow chart for generating a porous structure.


#### 2.1.1. Generation of 3D Wireframes

A complex parameterized 3D wireframe model was established by Grasshopper. Grasshopper is a node-based visualization language built into Rhinoceros. Designers can connect and adjust the visual operation nodes intuitively and perform parametric design quickly. Unlike the traditional programming environment, Grasshopper does not need to manually tap the code for editing. Instead, it provides another programming method that the internal battery components are connected to form a program according to a certain logic. This modeling method is more acceptable to designers, and it is found that this programming method is more efficient through practice [32].Grasshopper provides the IntraLattice plug-in, which contains many crystal structures and is convenient for the modeling of porous structures. The lattice wireframe battery generation graph is shown in Figure2.Figure 2
Lattice wireframe battery generation graph.Grasshopper provides 10 crystal structures, and the size and number of permutations of each lattice structure can be controlled as needed, which greatly improves the efficiency of drawing. Grasshopper provides a variety of unit cell structures. Taking into account the efficiency of subsequent finite element meshing and analysis, this paper uses four common lattice structures, namely,x-type, star-type, cross-type, and diamond-type. Figure 3 shows the structures of Grasshopper.Figure 3
Porous wireframe baking graph. (a)X. (b) Star. (c) Cross. (d) Diamond.
(a)(b)(c)(d)After the wireframe generated in Grasshopper is baked into Rhinoceros, it can support the conversion of multiple formats, which greatly facilitates the integration with Computer Aided Design (CAD) software. The wireframe can be saved in igs or x_t_ format for subsequent UG import operations.


#### 2.1.2. Expand 3D Wireframe into Solid Model

The 3D wireframe was expanded into a solid model through UG secondary development to complete the modeling of a porous structure. UG is a large-scale modeling and drawing software during the entire life cycle of a product. It has three major modules: CAD, computer aided engineering (CAE), and computer aided manufacturing (CAM). It also has a leading position in the global digital product development system and can significantly improve design efficiency and quality [33]. UG secondary development is divided into three major blocks: menus and toolbars, interface, and application [34]. For the first and second blocks, NXOpen provides NX Toolkits secondary development tools, including Block UI Styler and MenuScrip. According to the provided standard controls, Block UI Styler, as a visual user interface generation tool, can design an interface similar to the NX interactive dialog box, and a user can call the custom function file by editing the dialog box callback function. The Block UI Styler dialog box can be launched through the MenuScript menu bar. MenuScript is a unique scripting language for menu options in the NX system. Developers can find ∗ men files in the NX software installation directory and edit the contents of the files to modify or create new NX menus. At the same time, you can sort the original functions of the system or delete unnecessary function options in the system by modifying the ∗ men file, and you can also load your own designed menu options to customize the NX menus [35]. This article adopts the NX Open C++ application programming interface to program. Open C++ as an object-oriented API is fully compatible with NX Open, and programmers can access class hierarchies and even custom classes by referencing object features to establish permanent objects in NX. In this article, we will use NX Open API and NX Toolkits tools for the development of the program.After the menus and dialog boxes are created, NX Open C and NX Open C++ programming interfaces can be used to write programs in Visual Studio software and then generate dynamic link library files (∗ dll files). Clicking the dialog button will call the corresponding dynamic link library file, which achieves seamless linking of application development programs and the NX platform. In this article, the programming idea is as follows: randomly choose the imported wireframes, generate cylinders with a specified radius according to the length of each wireframe, and finally merge each cylinder into a single body. Its core functions are as follows: first, obtain the data of the selection dialog. One of the dialog boxes designed in this article is to get the selected line. Each line has a unique TAG value. Second, create cylinders. After the tag value of each line is obtained, the cylinder operation is carried out. Finally, merge all cylinders. Due to the irregularity of the wireframe, not every line has an intersection point, so a judgment distance function is needed.Through secondary development, it is possible to determine the diameter as needed to generate a cylinder from the wire frame and merge it into one body. Taking into account the size of the sole, the use of a too small cylindrical diameter will reduce the overall stiffness of the sole, and the overall shock absorption effect of the sole of a too large diameter is not obvious, so the diameters used in this article are 1 mm, 2 mm, and 3 mm, respectively. The porous cylinder structures are shown in Figure4.Figure 4
Different types of porous modeling graphs. (a)X. (b) Star. (c) Cross. (d) Diamond.
(a)(b)(c)(d)This method can quickly establish a porous model and greatly improve modeling efficiency. It can be applied to solid modeling which needs a porous structure. In this article, this program is used to establish the porous structure in the heel area of the middle sole of shoes to improve the shock absorption and cushioning performance of sports shoes.


### 2.2. Finite Element Models

#### 2.2.1. Theoretical Basis

Partial differential equations (PDEs) are fundamental to model different phenomena in science and engineering mathematically. Solving them is a crucial step towards a precise knowledge of the behavior of natural and engineered systems. Probably, the best-known option is the finite element method (FEM) for engineering problems. However, powerful alternatives such as mesh-free methods and isogeometric analysis (IGA) are also available. Samaniego et al. [36] explored deep neural networks (DNNs) as an option for approximation which could deal with such problems without any discretization. They have shown impressive results in areas such as visual recognition. Finite element analysis method is used in this paper. According to the previous research, all the components in the model are regarded as homogeneous and isotropic linear elastic materials.The basic idea of finite element is to discretize the continuum, then seek for the governing equation and satisfy the boundary conditions and continuity conditions, and then solve the numerical solution. Different from the static problem, the inertial force, damping force, and other terms are added in the basic equation of dynamic finite element; that is, the acceleration and the load related to velocity are described by adding discrete and consistent time measures of structure mass distribution, which can be expressed as follows:Equilibrium equation:(1)σij,j+bi−ρui¨−μui˙=0.Geometric equation:(2)εij=12ui,j+uj,i.Physical equation:(3)σij=Dijklεkl.Boundary conditions:(4)ui=ui¯on the boundary ofSu,σijnj=qion the boundary ofSσ.Initial conditions:(5)uix,y,z,0=ui0x,y,z,ui˙x,y,z,0=ui˙0x,y,z,where ρ is the mass density; u is the damping factor; Su is the displacement boundary condition; Sσ is the force boundary condition; σij is the stress components at any point in an elastic body; bi is the component of a volumetric force vector of an elastic body along the rectangular axis; qi is the the component of the surface force vector of an elastic body along the rectangular axis; εij is the component of a strain state at any point in an elastic body; Dijkl is the tensor representation of an elastic matrix; ui is the component of a displacement vector at a point in an elastic body along the rectangular axis;  nj is the direction cosine of the angle between the outer normal of the elastomer boundary Sσ and the three axes.In dynamic analysis, time coordinate is introduced, and space and time discrete method is generally used in finite element analysis. For spatial dispersion, the corresponding difference function is as follows:(6)uix,y,z,t=∑k=1nNkuikt,where n is the number of element nodes, and the node displacement parameter uikt is a function of time. With the help of the variational method, the governing equation of dynamic problem analysis can be obtained:(7)Mu¨t+Cu˙t+Kut=Ft,where M, C, K, and Ft are respectively the mass matrix, damping matrix, stiffness matrix, and node load vector of the system, which are respectively obtained by assembling the matrix and vector of their respective elements and expressed as(8)M=∑e=1nMe,C=∑e=1nCe,K=∑e=1nKe,F=∑e=1nFe,Me=∫VeρNTNdV,Ce=∫VeμNTNdV,Ke=∫VeBTNdV,Fe=∫VeNTfdV+∫SσeNTNdS.where N is the element shape function matrix; B is the element geometric matrix, which represents the relation between strain and node displacement; f is the load vector of the element node.If the effect of damping is ignored, the equation of motion can be abbreviated as,(9)Mu¨t+Kut=Ft.By solving and calculating the equation according to a certain method, the solution of the motion and mechanical parameters of the analyzed system can be obtained.


#### 2.2.2. Finite Element Model of Foot

Philips Brilliance 64-slice spiral CT was used to scan the feet of the volunteers, who weighed 63.5 kg. X-ray examination confirmed that the feet had no bone damage such as fracture, deformity, and tumor. The obtained CT image data was output in DICOM format, and then imported into MIMICS software. The corresponding threshold value was selected for the extraction of bone, and corresponding processing was calculated to generate a rough surface of the foot bone 3D model. The model was further processed through the Geomagic Studio software and UG software, and the final skeleton model is shown in Figure5.Figure 5
Foot bone model.The 3D model was imported into hypermesh for mesh division, material attribute assignment, setting of boundary conditions and loads, and creation of analysis steps. In the transmission of foot force, the skeleton acts as the main support of the human body, which plays a major role in supporting the weight. Therefore, at the lower leg of the foot, the tibia and fibula become the main force points of the human body gravity. In this section, when simulating the neutral condition of human feet, the position of applying gravity load is the upper face of the tibia and fibula. According to previous research, the force on the fibula accounts for about 1/4 of the force on the lower leg, and Ding and Song [37] showed that the force on the fibula accounted for about 1/6 of the force on the lower leg. Although there are some differences in studies, such differences are normal considering that different body positions and body shapes lead to different centers of gravity. The compromise in this section is that the tibia bears 80% of the gravity load and the fibula bears 20%, the resulting foot load is as follows:(i)
Step 1: apply 10 N pressure to the upper end of the tibia and fibula, respectively.(ii)
Step 2: apply pressure of 50 N on the upper tibia and 25 N on the upper fibula.(iii)
Step 3: apply pressure of 235 N on the upper tibia and 50 N on the upper fibula.The reason of three-step loading is that direct loading may cause too much deformation of grid cells and lead to nonconvergence of the calculation results.The complete finite element model of the foot is shown in Figure6.Figure 6
Finite element model of the foot.Table1 shows the material properties of each organizational structure of the foot model [38].Table 1
Material parameters of foot tissues.
NameElasticity modulus (MPa)Poisson’s ratioDensity (kg/m3)Bone73000.31500Joint0.150.42000Ligament2600.48937Plantar ligaments8000.45937Soft tissue0.450.48937Sole40.4123Among them, the ligament adopts beam element, and the bottom plate is used to simulate the ground. The coupling relationship between various structures of the foot model is set as face-to-face contact. The processed model was imported into ABAQUS for analysis and calculation, and the stress analysis results are shown in Figure7. We can see that the maximum stress occurs at the bottom of the calcaneus, with a value of 26.06 MPa. The bottom of the calcaneus is the main stress bearing area. According to the research results of Bayod et al. [39], under a gravity of 300 N, the extreme stress of the foot appears at the contact point between the calcaneal bone and the ground. Moreover, according to the research conclusions of Wang [40], the maximum stress of the calcaneal bone of the foot when the human body is standing normally is 16.608 MPa. In terms of stress distribution, it is similar to the research results of Taha et al. [41], so the complete foot model is true and reliable.Figure 7
Stress nephogram of the foot skeleton.


#### 2.2.3. Finite Element Model of Shoe Midsole

In this paper, porous structure is filled in the midsole heel area of the shoe. Cross,X, star, and diamond crystal structure wireframes are used, respectively. Each crystal wireframe is endowed with a diameter of 1, 2, and 3 mm by the method in Section 2.1 to form a porous solid structure, as is shown in Figures 8 and 9.Figure 8
Overall view of porous soles.Figure 9
Partial view of the perforated sole.The finite element method adopts the numerical discrete method for partial differential equations. Theoretically, the finer the mesh size is divided, the higher the solution precision is. However, the grid cannot be unlimitedly encrypted; there are mainly the following contradictions:(1)
The greater the grid density, the more the cells and nodes, the more the computation, and the longer the computation time.(2)
The rounding error caused by computer floating point operation will increase due to the gradual increase of meshing density.Therefore, in practical engineering problems, users often seek a more appropriate point between the finite element solution accuracy and the calculation cost, and the location of the point is to reach the threshold of grid density independence. In view of this, this section makes a series of meshes for the sole, and the bottom is fixed. A force of 300 N is applied to the upper surface of the sole to simulate the pressure of the foot on the sole. The analysis convergence results are shown in Figure10.Figure 10
Convergence results of soles with different number of elements.As can be seen from Figure10, with the increase of meshing density, the maximum stress of calcaneus gradually increases, and the maximum displacement curve gradually improves. Moreover, when the mesh size starts to be less than 4 mm, the number of the elements is 246274, the convergence result gradually flattens out, and the fluctuation is not large. The final finite element model of foot-shoe is shown in Figure 11.Figure 11
The finite element model of foot-sole.


#### 2.2.4. Dynamic Analysis of Foot-Shoe

ABAQUS/Explicit analysis is suitable for large models with short dynamic response times and extremely discrete events or processes. At the same time, it has a great advantage in dealing with the contact nonlinear problem. It can find the contact among the structures in the model by itself and simulate the complex contact relationship among the structures quickly and effectively.At the beginning of incremental complement (timet), according to the node dynamic equilibrium equation,(10)Mu¨=P−I.Calculated node acceleration:(11)u¨|t=M−1P−I,where M is the mass matrix of the node; u is the acceleration of the node; P is the external force applied; and I is the internal force of the element.The acceleration of any node can be obtained by the mass of the node and the resultant force acting on it. Then, the central difference method is used to determine the velocity of the node in the current increment step:(12)u¨|t+Δt/2=u¨|t−Δt/2+Δt|t+Δt+Δt|t2u¨|t.Thus, the displacement at the end of the increment step can be obtained:(13)u|t+Δt=u|t+Δt|t+Δtu˙|t+Δt/2.Finally, the element stress is calculated according to the strain rate and constitutive relation, and then the internal force of the element node is obtained.Compared with the implicit algorithm, the biggest difference between the two methods lies in the way to solve the node acceleration. For implicit algorithm, it iterates the next result by using the result of the current step and the unknown result of the next step, which leads to high computational cost. The explicit algorithm calculates the result of the previous step and the result of the current step. It is also based on the above characteristics that explicit methods have unique advantages in solving high-speed dynamic events and some extremely discontinuous cases.The element selection, mesh generation, and contact settings in this section are similar to those in 2.2.1. The difference is that under the dynamic analysis, the material of the model must be endowed with the density attribute to consider the influence of inertial forces, and the specific corresponding values are referred to Table1 in Section 3. At the same time, according to the working condition of full foot landing of human body, the boundary conditions are set as follows: Step 1: apply pressure with a total force of 10 N on the two surfaces, respectively, and the time T1 is 0∼10-6 s. Step 2: apply pressure of 50 N on the upper tibia and 25 N on the upper fibula for time T2 of 10-6∼10-5 s. Step 3: apply pressure of 235 N on the upper tibia and 50 N on the upper fibula for time T3 of 10-5∼0.0095 s. As the foot has a certain initial velocity before landing, according to Qian’s research [38], the initial velocity is applied to the whole foot: Vx = 1.65 m/s; vz = 0.23 m/s.


## 2.1. Porous Modeling

In this section, Grasshopper and UG are used to establish the porous structure. Figure1 is the process flow chart.Figure 1
A flow chart for generating a porous structure.


### 2.1.1. Generation of 3D Wireframes

A complex parameterized 3D wireframe model was established by Grasshopper. Grasshopper is a node-based visualization language built into Rhinoceros. Designers can connect and adjust the visual operation nodes intuitively and perform parametric design quickly. Unlike the traditional programming environment, Grasshopper does not need to manually tap the code for editing. Instead, it provides another programming method that the internal battery components are connected to form a program according to a certain logic. This modeling method is more acceptable to designers, and it is found that this programming method is more efficient through practice [32].Grasshopper provides the IntraLattice plug-in, which contains many crystal structures and is convenient for the modeling of porous structures. The lattice wireframe battery generation graph is shown in Figure2.Figure 2
Lattice wireframe battery generation graph.Grasshopper provides 10 crystal structures, and the size and number of permutations of each lattice structure can be controlled as needed, which greatly improves the efficiency of drawing. Grasshopper provides a variety of unit cell structures. Taking into account the efficiency of subsequent finite element meshing and analysis, this paper uses four common lattice structures, namely,x-type, star-type, cross-type, and diamond-type. Figure 3 shows the structures of Grasshopper.Figure 3
Porous wireframe baking graph. (a)X. (b) Star. (c) Cross. (d) Diamond.
(a)(b)(c)(d)After the wireframe generated in Grasshopper is baked into Rhinoceros, it can support the conversion of multiple formats, which greatly facilitates the integration with Computer Aided Design (CAD) software. The wireframe can be saved in igs or x_t_ format for subsequent UG import operations.


### 2.1.2. Expand 3D Wireframe into Solid Model

The 3D wireframe was expanded into a solid model through UG secondary development to complete the modeling of a porous structure. UG is a large-scale modeling and drawing software during the entire life cycle of a product. It has three major modules: CAD, computer aided engineering (CAE), and computer aided manufacturing (CAM). It also has a leading position in the global digital product development system and can significantly improve design efficiency and quality [33]. UG secondary development is divided into three major blocks: menus and toolbars, interface, and application [34]. For the first and second blocks, NXOpen provides NX Toolkits secondary development tools, including Block UI Styler and MenuScrip. According to the provided standard controls, Block UI Styler, as a visual user interface generation tool, can design an interface similar to the NX interactive dialog box, and a user can call the custom function file by editing the dialog box callback function. The Block UI Styler dialog box can be launched through the MenuScript menu bar. MenuScript is a unique scripting language for menu options in the NX system. Developers can find ∗ men files in the NX software installation directory and edit the contents of the files to modify or create new NX menus. At the same time, you can sort the original functions of the system or delete unnecessary function options in the system by modifying the ∗ men file, and you can also load your own designed menu options to customize the NX menus [35]. This article adopts the NX Open C++ application programming interface to program. Open C++ as an object-oriented API is fully compatible with NX Open, and programmers can access class hierarchies and even custom classes by referencing object features to establish permanent objects in NX. In this article, we will use NX Open API and NX Toolkits tools for the development of the program.After the menus and dialog boxes are created, NX Open C and NX Open C++ programming interfaces can be used to write programs in Visual Studio software and then generate dynamic link library files (∗ dll files). Clicking the dialog button will call the corresponding dynamic link library file, which achieves seamless linking of application development programs and the NX platform. In this article, the programming idea is as follows: randomly choose the imported wireframes, generate cylinders with a specified radius according to the length of each wireframe, and finally merge each cylinder into a single body. Its core functions are as follows: first, obtain the data of the selection dialog. One of the dialog boxes designed in this article is to get the selected line. Each line has a unique TAG value. Second, create cylinders. After the tag value of each line is obtained, the cylinder operation is carried out. Finally, merge all cylinders. Due to the irregularity of the wireframe, not every line has an intersection point, so a judgment distance function is needed.Through secondary development, it is possible to determine the diameter as needed to generate a cylinder from the wire frame and merge it into one body. Taking into account the size of the sole, the use of a too small cylindrical diameter will reduce the overall stiffness of the sole, and the overall shock absorption effect of the sole of a too large diameter is not obvious, so the diameters used in this article are 1 mm, 2 mm, and 3 mm, respectively. The porous cylinder structures are shown in Figure4.Figure 4
Different types of porous modeling graphs. (a)X. (b) Star. (c) Cross. (d) Diamond.
(a)(b)(c)(d)This method can quickly establish a porous model and greatly improve modeling efficiency. It can be applied to solid modeling which needs a porous structure. In this article, this program is used to establish the porous structure in the heel area of the middle sole of shoes to improve the shock absorption and cushioning performance of sports shoes.


## 2.1.1. Generation of 3D Wireframes

A complex parameterized 3D wireframe model was established by Grasshopper. Grasshopper is a node-based visualization language built into Rhinoceros. Designers can connect and adjust the visual operation nodes intuitively and perform parametric design quickly. Unlike the traditional programming environment, Grasshopper does not need to manually tap the code for editing. Instead, it provides another programming method that the internal battery components are connected to form a program according to a certain logic. This modeling method is more acceptable to designers, and it is found that this programming method is more efficient through practice [32].Grasshopper provides the IntraLattice plug-in, which contains many crystal structures and is convenient for the modeling of porous structures. The lattice wireframe battery generation graph is shown in Figure2.Figure 2
Lattice wireframe battery generation graph.Grasshopper provides 10 crystal structures, and the size and number of permutations of each lattice structure can be controlled as needed, which greatly improves the efficiency of drawing. Grasshopper provides a variety of unit cell structures. Taking into account the efficiency of subsequent finite element meshing and analysis, this paper uses four common lattice structures, namely,x-type, star-type, cross-type, and diamond-type. Figure 3 shows the structures of Grasshopper.Figure 3
Porous wireframe baking graph. (a)X. (b) Star. (c) Cross. (d) Diamond.
(a)(b)(c)(d)After the wireframe generated in Grasshopper is baked into Rhinoceros, it can support the conversion of multiple formats, which greatly facilitates the integration with Computer Aided Design (CAD) software. The wireframe can be saved in igs or x_t_ format for subsequent UG import operations.


## 2.1.2. Expand 3D Wireframe into Solid Model

The 3D wireframe was expanded into a solid model through UG secondary development to complete the modeling of a porous structure. UG is a large-scale modeling and drawing software during the entire life cycle of a product. It has three major modules: CAD, computer aided engineering (CAE), and computer aided manufacturing (CAM). It also has a leading position in the global digital product development system and can significantly improve design efficiency and quality [33]. UG secondary development is divided into three major blocks: menus and toolbars, interface, and application [34]. For the first and second blocks, NXOpen provides NX Toolkits secondary development tools, including Block UI Styler and MenuScrip. According to the provided standard controls, Block UI Styler, as a visual user interface generation tool, can design an interface similar to the NX interactive dialog box, and a user can call the custom function file by editing the dialog box callback function. The Block UI Styler dialog box can be launched through the MenuScript menu bar. MenuScript is a unique scripting language for menu options in the NX system. Developers can find ∗ men files in the NX software installation directory and edit the contents of the files to modify or create new NX menus. At the same time, you can sort the original functions of the system or delete unnecessary function options in the system by modifying the ∗ men file, and you can also load your own designed menu options to customize the NX menus [35]. This article adopts the NX Open C++ application programming interface to program. Open C++ as an object-oriented API is fully compatible with NX Open, and programmers can access class hierarchies and even custom classes by referencing object features to establish permanent objects in NX. In this article, we will use NX Open API and NX Toolkits tools for the development of the program.After the menus and dialog boxes are created, NX Open C and NX Open C++ programming interfaces can be used to write programs in Visual Studio software and then generate dynamic link library files (∗ dll files). Clicking the dialog button will call the corresponding dynamic link library file, which achieves seamless linking of application development programs and the NX platform. In this article, the programming idea is as follows: randomly choose the imported wireframes, generate cylinders with a specified radius according to the length of each wireframe, and finally merge each cylinder into a single body. Its core functions are as follows: first, obtain the data of the selection dialog. One of the dialog boxes designed in this article is to get the selected line. Each line has a unique TAG value. Second, create cylinders. After the tag value of each line is obtained, the cylinder operation is carried out. Finally, merge all cylinders. Due to the irregularity of the wireframe, not every line has an intersection point, so a judgment distance function is needed.Through secondary development, it is possible to determine the diameter as needed to generate a cylinder from the wire frame and merge it into one body. Taking into account the size of the sole, the use of a too small cylindrical diameter will reduce the overall stiffness of the sole, and the overall shock absorption effect of the sole of a too large diameter is not obvious, so the diameters used in this article are 1 mm, 2 mm, and 3 mm, respectively. The porous cylinder structures are shown in Figure4.Figure 4
Different types of porous modeling graphs. (a)X. (b) Star. (c) Cross. (d) Diamond.
(a)(b)(c)(d)This method can quickly establish a porous model and greatly improve modeling efficiency. It can be applied to solid modeling which needs a porous structure. In this article, this program is used to establish the porous structure in the heel area of the middle sole of shoes to improve the shock absorption and cushioning performance of sports shoes.


## 2.2. Finite Element Models

### 2.2.1. Theoretical Basis

Partial differential equations (PDEs) are fundamental to model different phenomena in science and engineering mathematically. Solving them is a crucial step towards a precise knowledge of the behavior of natural and engineered systems. Probably, the best-known option is the finite element method (FEM) for engineering problems. However, powerful alternatives such as mesh-free methods and isogeometric analysis (IGA) are also available. Samaniego et al. [36] explored deep neural networks (DNNs) as an option for approximation which could deal with such problems without any discretization. They have shown impressive results in areas such as visual recognition. Finite element analysis method is used in this paper. According to the previous research, all the components in the model are regarded as homogeneous and isotropic linear elastic materials.The basic idea of finite element is to discretize the continuum, then seek for the governing equation and satisfy the boundary conditions and continuity conditions, and then solve the numerical solution. Different from the static problem, the inertial force, damping force, and other terms are added in the basic equation of dynamic finite element; that is, the acceleration and the load related to velocity are described by adding discrete and consistent time measures of structure mass distribution, which can be expressed as follows:Equilibrium equation:(1)σij,j+bi−ρui¨−μui˙=0.Geometric equation:(2)εij=12ui,j+uj,i.Physical equation:(3)σij=Dijklεkl.Boundary conditions:(4)ui=ui¯on the boundary ofSu,σijnj=qion the boundary ofSσ.Initial conditions:(5)uix,y,z,0=ui0x,y,z,ui˙x,y,z,0=ui˙0x,y,z,where ρ is the mass density; u is the damping factor; Su is the displacement boundary condition; Sσ is the force boundary condition; σij is the stress components at any point in an elastic body; bi is the component of a volumetric force vector of an elastic body along the rectangular axis; qi is the the component of the surface force vector of an elastic body along the rectangular axis; εij is the component of a strain state at any point in an elastic body; Dijkl is the tensor representation of an elastic matrix; ui is the component of a displacement vector at a point in an elastic body along the rectangular axis;  nj is the direction cosine of the angle between the outer normal of the elastomer boundary Sσ and the three axes.In dynamic analysis, time coordinate is introduced, and space and time discrete method is generally used in finite element analysis. For spatial dispersion, the corresponding difference function is as follows:(6)uix,y,z,t=∑k=1nNkuikt,where n is the number of element nodes, and the node displacement parameter uikt is a function of time. With the help of the variational method, the governing equation of dynamic problem analysis can be obtained:(7)Mu¨t+Cu˙t+Kut=Ft,where M, C, K, and Ft are respectively the mass matrix, damping matrix, stiffness matrix, and node load vector of the system, which are respectively obtained by assembling the matrix and vector of their respective elements and expressed as(8)M=∑e=1nMe,C=∑e=1nCe,K=∑e=1nKe,F=∑e=1nFe,Me=∫VeρNTNdV,Ce=∫VeμNTNdV,Ke=∫VeBTNdV,Fe=∫VeNTfdV+∫SσeNTNdS.where N is the element shape function matrix; B is the element geometric matrix, which represents the relation between strain and node displacement; f is the load vector of the element node.If the effect of damping is ignored, the equation of motion can be abbreviated as,(9)Mu¨t+Kut=Ft.By solving and calculating the equation according to a certain method, the solution of the motion and mechanical parameters of the analyzed system can be obtained.


### 2.2.2. Finite Element Model of Foot

Philips Brilliance 64-slice spiral CT was used to scan the feet of the volunteers, who weighed 63.5 kg. X-ray examination confirmed that the feet had no bone damage such as fracture, deformity, and tumor. The obtained CT image data was output in DICOM format, and then imported into MIMICS software. The corresponding threshold value was selected for the extraction of bone, and corresponding processing was calculated to generate a rough surface of the foot bone 3D model. The model was further processed through the Geomagic Studio software and UG software, and the final skeleton model is shown in Figure5.Figure 5
Foot bone model.The 3D model was imported into hypermesh for mesh division, material attribute assignment, setting of boundary conditions and loads, and creation of analysis steps. In the transmission of foot force, the skeleton acts as the main support of the human body, which plays a major role in supporting the weight. Therefore, at the lower leg of the foot, the tibia and fibula become the main force points of the human body gravity. In this section, when simulating the neutral condition of human feet, the position of applying gravity load is the upper face of the tibia and fibula. According to previous research, the force on the fibula accounts for about 1/4 of the force on the lower leg, and Ding and Song [37] showed that the force on the fibula accounted for about 1/6 of the force on the lower leg. Although there are some differences in studies, such differences are normal considering that different body positions and body shapes lead to different centers of gravity. The compromise in this section is that the tibia bears 80% of the gravity load and the fibula bears 20%, the resulting foot load is as follows:(i)
Step 1: apply 10 N pressure to the upper end of the tibia and fibula, respectively.(ii)
Step 2: apply pressure of 50 N on the upper tibia and 25 N on the upper fibula.(iii)
Step 3: apply pressure of 235 N on the upper tibia and 50 N on the upper fibula.The reason of three-step loading is that direct loading may cause too much deformation of grid cells and lead to nonconvergence of the calculation results.The complete finite element model of the foot is shown in Figure6.Figure 6
Finite element model of the foot.Table1 shows the material properties of each organizational structure of the foot model [38].Table 1
Material parameters of foot tissues.
NameElasticity modulus (MPa)Poisson’s ratioDensity (kg/m3)Bone73000.31500Joint0.150.42000Ligament2600.48937Plantar ligaments8000.45937Soft tissue0.450.48937Sole40.4123Among them, the ligament adopts beam element, and the bottom plate is used to simulate the ground. The coupling relationship between various structures of the foot model is set as face-to-face contact. The processed model was imported into ABAQUS for analysis and calculation, and the stress analysis results are shown in Figure7. We can see that the maximum stress occurs at the bottom of the calcaneus, with a value of 26.06 MPa. The bottom of the calcaneus is the main stress bearing area. According to the research results of Bayod et al. [39], under a gravity of 300 N, the extreme stress of the foot appears at the contact point between the calcaneal bone and the ground. Moreover, according to the research conclusions of Wang [40], the maximum stress of the calcaneal bone of the foot when the human body is standing normally is 16.608 MPa. In terms of stress distribution, it is similar to the research results of Taha et al. [41], so the complete foot model is true and reliable.Figure 7
Stress nephogram of the foot skeleton.


### 2.2.3. Finite Element Model of Shoe Midsole

In this paper, porous structure is filled in the midsole heel area of the shoe. Cross,X, star, and diamond crystal structure wireframes are used, respectively. Each crystal wireframe is endowed with a diameter of 1, 2, and 3 mm by the method in Section 2.1 to form a porous solid structure, as is shown in Figures 8 and 9.Figure 8
Overall view of porous soles.Figure 9
Partial view of the perforated sole.The finite element method adopts the numerical discrete method for partial differential equations. Theoretically, the finer the mesh size is divided, the higher the solution precision is. However, the grid cannot be unlimitedly encrypted; there are mainly the following contradictions:(1)
The greater the grid density, the more the cells and nodes, the more the computation, and the longer the computation time.(2)
The rounding error caused by computer floating point operation will increase due to the gradual increase of meshing density.Therefore, in practical engineering problems, users often seek a more appropriate point between the finite element solution accuracy and the calculation cost, and the location of the point is to reach the threshold of grid density independence. In view of this, this section makes a series of meshes for the sole, and the bottom is fixed. A force of 300 N is applied to the upper surface of the sole to simulate the pressure of the foot on the sole. The analysis convergence results are shown in Figure10.Figure 10
Convergence results of soles with different number of elements.As can be seen from Figure10, with the increase of meshing density, the maximum stress of calcaneus gradually increases, and the maximum displacement curve gradually improves. Moreover, when the mesh size starts to be less than 4 mm, the number of the elements is 246274, the convergence result gradually flattens out, and the fluctuation is not large. The final finite element model of foot-shoe is shown in Figure 11.Figure 11
The finite element model of foot-sole.


### 2.2.4. Dynamic Analysis of Foot-Shoe

ABAQUS/Explicit analysis is suitable for large models with short dynamic response times and extremely discrete events or processes. At the same time, it has a great advantage in dealing with the contact nonlinear problem. It can find the contact among the structures in the model by itself and simulate the complex contact relationship among the structures quickly and effectively.At the beginning of incremental complement (timet), according to the node dynamic equilibrium equation,(10)Mu¨=P−I.Calculated node acceleration:(11)u¨|t=M−1P−I,where M is the mass matrix of the node; u is the acceleration of the node; P is the external force applied; and I is the internal force of the element.The acceleration of any node can be obtained by the mass of the node and the resultant force acting on it. Then, the central difference method is used to determine the velocity of the node in the current increment step:(12)u¨|t+Δt/2=u¨|t−Δt/2+Δt|t+Δt+Δt|t2u¨|t.Thus, the displacement at the end of the increment step can be obtained:(13)u|t+Δt=u|t+Δt|t+Δtu˙|t+Δt/2.Finally, the element stress is calculated according to the strain rate and constitutive relation, and then the internal force of the element node is obtained.Compared with the implicit algorithm, the biggest difference between the two methods lies in the way to solve the node acceleration. For implicit algorithm, it iterates the next result by using the result of the current step and the unknown result of the next step, which leads to high computational cost. The explicit algorithm calculates the result of the previous step and the result of the current step. It is also based on the above characteristics that explicit methods have unique advantages in solving high-speed dynamic events and some extremely discontinuous cases.The element selection, mesh generation, and contact settings in this section are similar to those in 2.2.1. The difference is that under the dynamic analysis, the material of the model must be endowed with the density attribute to consider the influence of inertial forces, and the specific corresponding values are referred to Table1 in Section 3. At the same time, according to the working condition of full foot landing of human body, the boundary conditions are set as follows: Step 1: apply pressure with a total force of 10 N on the two surfaces, respectively, and the time T1 is 0∼10-6 s. Step 2: apply pressure of 50 N on the upper tibia and 25 N on the upper fibula for time T2 of 10-6∼10-5 s. Step 3: apply pressure of 235 N on the upper tibia and 50 N on the upper fibula for time T3 of 10-5∼0.0095 s. As the foot has a certain initial velocity before landing, according to Qian’s research [38], the initial velocity is applied to the whole foot: Vx = 1.65 m/s; vz = 0.23 m/s.


## 2.2.1. Theoretical Basis

Partial differential equations (PDEs) are fundamental to model different phenomena in science and engineering mathematically. Solving them is a crucial step towards a precise knowledge of the behavior of natural and engineered systems. Probably, the best-known option is the finite element method (FEM) for engineering problems. However, powerful alternatives such as mesh-free methods and isogeometric analysis (IGA) are also available. Samaniego et al. [36] explored deep neural networks (DNNs) as an option for approximation which could deal with such problems without any discretization. They have shown impressive results in areas such as visual recognition. Finite element analysis method is used in this paper. According to the previous research, all the components in the model are regarded as homogeneous and isotropic linear elastic materials.The basic idea of finite element is to discretize the continuum, then seek for the governing equation and satisfy the boundary conditions and continuity conditions, and then solve the numerical solution. Different from the static problem, the inertial force, damping force, and other terms are added in the basic equation of dynamic finite element; that is, the acceleration and the load related to velocity are described by adding discrete and consistent time measures of structure mass distribution, which can be expressed as follows:Equilibrium equation:(1)σij,j+bi−ρui¨−μui˙=0.Geometric equation:(2)εij=12ui,j+uj,i.Physical equation:(3)σij=Dijklεkl.Boundary conditions:(4)ui=ui¯on the boundary ofSu,σijnj=qion the boundary ofSσ.Initial conditions:(5)uix,y,z,0=ui0x,y,z,ui˙x,y,z,0=ui˙0x,y,z,where ρ is the mass density; u is the damping factor; Su is the displacement boundary condition; Sσ is the force boundary condition; σij is the stress components at any point in an elastic body; bi is the component of a volumetric force vector of an elastic body along the rectangular axis; qi is the the component of the surface force vector of an elastic body along the rectangular axis; εij is the component of a strain state at any point in an elastic body; Dijkl is the tensor representation of an elastic matrix; ui is the component of a displacement vector at a point in an elastic body along the rectangular axis;  nj is the direction cosine of the angle between the outer normal of the elastomer boundary Sσ and the three axes.In dynamic analysis, time coordinate is introduced, and space and time discrete method is generally used in finite element analysis. For spatial dispersion, the corresponding difference function is as follows:(6)uix,y,z,t=∑k=1nNkuikt,where n is the number of element nodes, and the node displacement parameter uikt is a function of time. With the help of the variational method, the governing equation of dynamic problem analysis can be obtained:(7)Mu¨t+Cu˙t+Kut=Ft,where M, C, K, and Ft are respectively the mass matrix, damping matrix, stiffness matrix, and node load vector of the system, which are respectively obtained by assembling the matrix and vector of their respective elements and expressed as(8)M=∑e=1nMe,C=∑e=1nCe,K=∑e=1nKe,F=∑e=1nFe,Me=∫VeρNTNdV,Ce=∫VeμNTNdV,Ke=∫VeBTNdV,Fe=∫VeNTfdV+∫SσeNTNdS.where N is the element shape function matrix; B is the element geometric matrix, which represents the relation between strain and node displacement; f is the load vector of the element node.If the effect of damping is ignored, the equation of motion can be abbreviated as,(9)Mu¨t+Kut=Ft.By solving and calculating the equation according to a certain method, the solution of the motion and mechanical parameters of the analyzed system can be obtained.


## 2.2.2. Finite Element Model of Foot

Philips Brilliance 64-slice spiral CT was used to scan the feet of the volunteers, who weighed 63.5 kg. X-ray examination confirmed that the feet had no bone damage such as fracture, deformity, and tumor. The obtained CT image data was output in DICOM format, and then imported into MIMICS software. The corresponding threshold value was selected for the extraction of bone, and corresponding processing was calculated to generate a rough surface of the foot bone 3D model. The model was further processed through the Geomagic Studio software and UG software, and the final skeleton model is shown in Figure5.Figure 5
Foot bone model.The 3D model was imported into hypermesh for mesh division, material attribute assignment, setting of boundary conditions and loads, and creation of analysis steps. In the transmission of foot force, the skeleton acts as the main support of the human body, which plays a major role in supporting the weight. Therefore, at the lower leg of the foot, the tibia and fibula become the main force points of the human body gravity. In this section, when simulating the neutral condition of human feet, the position of applying gravity load is the upper face of the tibia and fibula. According to previous research, the force on the fibula accounts for about 1/4 of the force on the lower leg, and Ding and Song [37] showed that the force on the fibula accounted for about 1/6 of the force on the lower leg. Although there are some differences in studies, such differences are normal considering that different body positions and body shapes lead to different centers of gravity. The compromise in this section is that the tibia bears 80% of the gravity load and the fibula bears 20%, the resulting foot load is as follows:(i)
Step 1: apply 10 N pressure to the upper end of the tibia and fibula, respectively.(ii)
Step 2: apply pressure of 50 N on the upper tibia and 25 N on the upper fibula.(iii)
Step 3: apply pressure of 235 N on the upper tibia and 50 N on the upper fibula.The reason of three-step loading is that direct loading may cause too much deformation of grid cells and lead to nonconvergence of the calculation results.The complete finite element model of the foot is shown in Figure6.Figure 6
Finite element model of the foot.Table1 shows the material properties of each organizational structure of the foot model [38].Table 1
Material parameters of foot tissues.
NameElasticity modulus (MPa)Poisson’s ratioDensity (kg/m3)Bone73000.31500Joint0.150.42000Ligament2600.48937Plantar ligaments8000.45937Soft tissue0.450.48937Sole40.4123Among them, the ligament adopts beam element, and the bottom plate is used to simulate the ground. The coupling relationship between various structures of the foot model is set as face-to-face contact. The processed model was imported into ABAQUS for analysis and calculation, and the stress analysis results are shown in Figure7. We can see that the maximum stress occurs at the bottom of the calcaneus, with a value of 26.06 MPa. The bottom of the calcaneus is the main stress bearing area. According to the research results of Bayod et al. [39], under a gravity of 300 N, the extreme stress of the foot appears at the contact point between the calcaneal bone and the ground. Moreover, according to the research conclusions of Wang [40], the maximum stress of the calcaneal bone of the foot when the human body is standing normally is 16.608 MPa. In terms of stress distribution, it is similar to the research results of Taha et al. [41], so the complete foot model is true and reliable.Figure 7
Stress nephogram of the foot skeleton.


## 2.2.3. Finite Element Model of Shoe Midsole

In this paper, porous structure is filled in the midsole heel area of the shoe. Cross,X, star, and diamond crystal structure wireframes are used, respectively. Each crystal wireframe is endowed with a diameter of 1, 2, and 3 mm by the method in Section 2.1 to form a porous solid structure, as is shown in Figures 8 and 9.Figure 8
Overall view of porous soles.Figure 9
Partial view of the perforated sole.The finite element method adopts the numerical discrete method for partial differential equations. Theoretically, the finer the mesh size is divided, the higher the solution precision is. However, the grid cannot be unlimitedly encrypted; there are mainly the following contradictions:(1)
The greater the grid density, the more the cells and nodes, the more the computation, and the longer the computation time.(2)
The rounding error caused by computer floating point operation will increase due to the gradual increase of meshing density.Therefore, in practical engineering problems, users often seek a more appropriate point between the finite element solution accuracy and the calculation cost, and the location of the point is to reach the threshold of grid density independence. In view of this, this section makes a series of meshes for the sole, and the bottom is fixed. A force of 300 N is applied to the upper surface of the sole to simulate the pressure of the foot on the sole. The analysis convergence results are shown in Figure10.Figure 10
Convergence results of soles with different number of elements.As can be seen from Figure10, with the increase of meshing density, the maximum stress of calcaneus gradually increases, and the maximum displacement curve gradually improves. Moreover, when the mesh size starts to be less than 4 mm, the number of the elements is 246274, the convergence result gradually flattens out, and the fluctuation is not large. The final finite element model of foot-shoe is shown in Figure 11.Figure 11
The finite element model of foot-sole.


## 2.2.4. Dynamic Analysis of Foot-Shoe

ABAQUS/Explicit analysis is suitable for large models with short dynamic response times and extremely discrete events or processes. At the same time, it has a great advantage in dealing with the contact nonlinear problem. It can find the contact among the structures in the model by itself and simulate the complex contact relationship among the structures quickly and effectively.At the beginning of incremental complement (timet), according to the node dynamic equilibrium equation,(10)Mu¨=P−I.Calculated node acceleration:(11)u¨|t=M−1P−I,where M is the mass matrix of the node; u is the acceleration of the node; P is the external force applied; and I is the internal force of the element.The acceleration of any node can be obtained by the mass of the node and the resultant force acting on it. Then, the central difference method is used to determine the velocity of the node in the current increment step:(12)u¨|t+Δt/2=u¨|t−Δt/2+Δt|t+Δt+Δt|t2u¨|t.Thus, the displacement at the end of the increment step can be obtained:(13)u|t+Δt=u|t+Δt|t+Δtu˙|t+Δt/2.Finally, the element stress is calculated according to the strain rate and constitutive relation, and then the internal force of the element node is obtained.Compared with the implicit algorithm, the biggest difference between the two methods lies in the way to solve the node acceleration. For implicit algorithm, it iterates the next result by using the result of the current step and the unknown result of the next step, which leads to high computational cost. The explicit algorithm calculates the result of the previous step and the result of the current step. It is also based on the above characteristics that explicit methods have unique advantages in solving high-speed dynamic events and some extremely discontinuous cases.The element selection, mesh generation, and contact settings in this section are similar to those in 2.2.1. The difference is that under the dynamic analysis, the material of the model must be endowed with the density attribute to consider the influence of inertial forces, and the specific corresponding values are referred to Table1 in Section 3. At the same time, according to the working condition of full foot landing of human body, the boundary conditions are set as follows: Step 1: apply pressure with a total force of 10 N on the two surfaces, respectively, and the time T1 is 0∼10-6 s. Step 2: apply pressure of 50 N on the upper tibia and 25 N on the upper fibula for time T2 of 10-6∼10-5 s. Step 3: apply pressure of 235 N on the upper tibia and 50 N on the upper fibula for time T3 of 10-5∼0.0095 s. As the foot has a certain initial velocity before landing, according to Qian’s research [38], the initial velocity is applied to the whole foot: Vx = 1.65 m/s; vz = 0.23 m/s.


## 3. Results and Discussion

In the process of movement, the function of the sole is to attenuate the reaction force of the ground to the human body. As the material or structure of the sole is elastic, the impact energy of the ground facing the sole is converted into the elastic potential energy of the sole. Thus, the more elastic the sole is, the greater the degree of compression can be, and the more energy is absorbed. Therefore, the compression of sole can reflect the cushioning effect of sole.The foot is an important part of human body bearing the load, and the improper distribution of plantar pressure will cause different degrees of damage to the human body. Through analyzing the singularity of various foot parameters caused by different pathologic conditions by the plantar pressure, we can analyze the cause of foot disease and the development trend of disease course and provide important reference for the therapeutic effect and functional evaluation. Therefore, the distribution of sole pressure is one of the indexes to evaluate the cushioning effect of sole.In recent years, biomechanical studies have found that the load of each metatarsal bone of the foot decreases from the first metatarsal bone to the fifth metatarsal bone in order (I 26%, II 20.2%, III 11.8%, IV 7.6%, and V 5.5%). As the first metatarsal plays an important role in weight bearing, excessive stress on the first metatarsal will cause pain in the first metatarsophalangeal joint, and severe cases will cause hallux valgus [42].The load on the second metatarsal is less than that on the first metatarsal and greater than that on the other three metatarsal bones, but the weight bearing area of the second metatarsal is half that of the first metatarsal, so the pressure of the second metatarsal is greater than that of the first metatarsal. It can be seen that the stress of the second metatarsal is the greatest when bearing load. If the second metatarsal takes on too much weight, smaller contact surfaces bearing more concentrated weight can lead to earlier and more severe clinical symptoms [43], and most of the metatarsal fractures found clinically occur in the second metatarsal bone.Therefore, in this study, the displacement of sole heel, the distribution of plantar pressure, and the peak stress of metatarsal bone are proposed as the evaluation indexes of the cushioning effect of sole.


### 3.1. Analysis of the Sole

#### 3.1.1. Displacement of Sole Compression

Figure12 shows the displacement of the sole heel. It can be seen that the deformation of the sole heel is relatively large. By comparing the deformation of the heel, the shock absorption effect of different porous soles is compared.Figure 12
The displacement nephogram of the sole.Figure13 shows the maximum displacement of the sole of different structure under the impact of the foot. It can be seen from the figure that all structures have the maximum deformation when the diameter of the porous cylinder is 1 mm, and diamond type has the maximum displacement, which is 4.725 mm.Figure 13
The maximum displacement of the soles with different porous structures.


#### 3.1.2. Strain Energy Analysis

Strain energy is the energy stored in the structure due to elastic deformation under the action of external forces. With the reduction of load, strain energy is gradually released outwards. The storage and release of strain energy reflect the cushioning effect of structure and material to a certain degree under the action of external load, which is an important embodiment of dynamic energy of structure. Figure14 shows the change of strain energy of different structure soles.Figure 14
Strain energy of soles with different porous structure.It can be seen from the figure that the variation trend of strain energy is consistent, and the strain energy of the sole starts to increase after 0.0015 s, indicating that the force is transferred from the top of the foot to the sole after 0.0015 s. It is found that the strain energy of the sole of the four kinds of porous structure is greater than that of the original sole, indicating that the elastic deformation of the sole of the porous structure is greater and the ability of absorbing foreign energy and releasing energy is more outstanding.


### 3.2. Plantar Pressure

The distribution of plantar pressure was calculated and analyzed by ABAQUS as shown in Figure15.Figure 15
The plantar pressure distribution.It can be seen from the figure that the plantar pressure is mainly distributed at the front palm and the heel, and the maximum stress occurs at the heel.It can be seen from Figure16 that the plantar pressure of the porous sole is significantly lower than that of the original sole model, and the plantar pressure of each porous structure cylinder with a diameter of 1 mm is significantly lower than that of 2 mm and 3 mm. When the radius of the porous cylinder was 1 mm, the peak plantar pressure of the four models reached 105.7 Pa, 105.1 Pa, 108.8 Pa, and 112.1 Pa, respectively. Compared with the other three models, the peak plantar pressure in the diamond model decreased by 0.6%, 3.4%, and 6.2%, respectively. From the above analysis, it can be concluded that the diamond porous model has better cushioning effect and effectively reduces the plantar pressure in the process of motion.Figure 16
The plantar pressure of the soles with different porous structures.


### 3.3. Peak Stress in the First and Second Metatarsal Bones

Figure17 shows the stress nephogram of the metatarsals. It can be seen from Figure 18 that the maximum stress variation rules of the first to second metatarsal bones of all models are consistent, and the metatarsal stress of the original sole is the greatest and the metatarsal stress of the four porous structures is significantly reduced. Four kinds of structure of the first metatarsal maximum stress were, respectively, 2.241, 2.30, 2.561, and 2.466 MPa. Compared with the other three models, the maximum stress in the first metatarsal of the diamond model decreased by 4.9%, 10.2%, and 6.7%, respectively. Four models of the second metatarsal maximum stress were, respectively, 3.416, 3.407, 3.508, and 3.427 MPa. Compared with the other three models, the maximum stress of the second metatarsal in the diamond model decreased by 0.3%, 2.9%, and 0.6%, respectively. By comparison, the diamond porous structure can better reduce the stress of metatarsal bones, indicating that the diamond structure sole has better shock absorption effect.Figure 17
The stress nephogram of the metatarsals.Figure 18
Comparison diagram of maximum stress in the metatarsal bone with different porous structures.Through comprehensive comparison, it is found that the porous structure of the sole can play a cushioning role on the foot. By comparing the deformation of the sole heel, the plantar pressure, and the peak stress of the first and second metatarsal bones, it is found that the diamond porous structure is the optimal damping structure.


## 3.1. Analysis of the Sole

### 3.1.1. Displacement of Sole Compression

Figure12 shows the displacement of the sole heel. It can be seen that the deformation of the sole heel is relatively large. By comparing the deformation of the heel, the shock absorption effect of different porous soles is compared.Figure 12
The displacement nephogram of the sole.Figure13 shows the maximum displacement of the sole of different structure under the impact of the foot. It can be seen from the figure that all structures have the maximum deformation when the diameter of the porous cylinder is 1 mm, and diamond type has the maximum displacement, which is 4.725 mm.Figure 13
The maximum displacement of the soles with different porous structures.


### 3.1.2. Strain Energy Analysis

Strain energy is the energy stored in the structure due to elastic deformation under the action of external forces. With the reduction of load, strain energy is gradually released outwards. The storage and release of strain energy reflect the cushioning effect of structure and material to a certain degree under the action of external load, which is an important embodiment of dynamic energy of structure. Figure14 shows the change of strain energy of different structure soles.Figure 14
Strain energy of soles with different porous structure.It can be seen from the figure that the variation trend of strain energy is consistent, and the strain energy of the sole starts to increase after 0.0015 s, indicating that the force is transferred from the top of the foot to the sole after 0.0015 s. It is found that the strain energy of the sole of the four kinds of porous structure is greater than that of the original sole, indicating that the elastic deformation of the sole of the porous structure is greater and the ability of absorbing foreign energy and releasing energy is more outstanding.


## 3.1.1. Displacement of Sole Compression

Figure12 shows the displacement of the sole heel. It can be seen that the deformation of the sole heel is relatively large. By comparing the deformation of the heel, the shock absorption effect of different porous soles is compared.Figure 12
The displacement nephogram of the sole.Figure13 shows the maximum displacement of the sole of different structure under the impact of the foot. It can be seen from the figure that all structures have the maximum deformation when the diameter of the porous cylinder is 1 mm, and diamond type has the maximum displacement, which is 4.725 mm.Figure 13
The maximum displacement of the soles with different porous structures.


## 3.1.2. Strain Energy Analysis

Strain energy is the energy stored in the structure due to elastic deformation under the action of external forces. With the reduction of load, strain energy is gradually released outwards. The storage and release of strain energy reflect the cushioning effect of structure and material to a certain degree under the action of external load, which is an important embodiment of dynamic energy of structure. Figure14 shows the change of strain energy of different structure soles.Figure 14
Strain energy of soles with different porous structure.It can be seen from the figure that the variation trend of strain energy is consistent, and the strain energy of the sole starts to increase after 0.0015 s, indicating that the force is transferred from the top of the foot to the sole after 0.0015 s. It is found that the strain energy of the sole of the four kinds of porous structure is greater than that of the original sole, indicating that the elastic deformation of the sole of the porous structure is greater and the ability of absorbing foreign energy and releasing energy is more outstanding.


## 3.2. Plantar Pressure

The distribution of plantar pressure was calculated and analyzed by ABAQUS as shown in Figure15.Figure 15
The plantar pressure distribution.It can be seen from the figure that the plantar pressure is mainly distributed at the front palm and the heel, and the maximum stress occurs at the heel.It can be seen from Figure16 that the plantar pressure of the porous sole is significantly lower than that of the original sole model, and the plantar pressure of each porous structure cylinder with a diameter of 1 mm is significantly lower than that of 2 mm and 3 mm. When the radius of the porous cylinder was 1 mm, the peak plantar pressure of the four models reached 105.7 Pa, 105.1 Pa, 108.8 Pa, and 112.1 Pa, respectively. Compared with the other three models, the peak plantar pressure in the diamond model decreased by 0.6%, 3.4%, and 6.2%, respectively. From the above analysis, it can be concluded that the diamond porous model has better cushioning effect and effectively reduces the plantar pressure in the process of motion.Figure 16
The plantar pressure of the soles with different porous structures.


## 3.3. Peak Stress in the First and Second Metatarsal Bones

Figure17 shows the stress nephogram of the metatarsals. It can be seen from Figure 18 that the maximum stress variation rules of the first to second metatarsal bones of all models are consistent, and the metatarsal stress of the original sole is the greatest and the metatarsal stress of the four porous structures is significantly reduced. Four kinds of structure of the first metatarsal maximum stress were, respectively, 2.241, 2.30, 2.561, and 2.466 MPa. Compared with the other three models, the maximum stress in the first metatarsal of the diamond model decreased by 4.9%, 10.2%, and 6.7%, respectively. Four models of the second metatarsal maximum stress were, respectively, 3.416, 3.407, 3.508, and 3.427 MPa. Compared with the other three models, the maximum stress of the second metatarsal in the diamond model decreased by 0.3%, 2.9%, and 0.6%, respectively. By comparison, the diamond porous structure can better reduce the stress of metatarsal bones, indicating that the diamond structure sole has better shock absorption effect.Figure 17
The stress nephogram of the metatarsals.Figure 18
Comparison diagram of maximum stress in the metatarsal bone with different porous structures.Through comprehensive comparison, it is found that the porous structure of the sole can play a cushioning role on the foot. By comparing the deformation of the sole heel, the plantar pressure, and the peak stress of the first and second metatarsal bones, it is found that the diamond porous structure is the optimal damping structure.


## 4. Conclusion

In this paper, a porous modeling method for the sole is proposed. Using Grasshopper and UG secondary development tools, models with porous structure can be easily constructed. By changing and controlling the specific parameters, the characteristic size of the porous structure can be changed and the modeling efficiency of the porous structure can be improved. At the same time, the finite element analysis and optimization of porous structure combined with CAE provide a technical basis for the preparation of porous materials with excellent properties. Twelve kinds of porous sole models were established by this method, and their damping effect was analyzed by finite element method. Three indexes, the deformation of sole, the peak of plantar pressure, and the peak of metatarsal stress were used to evaluate the damping effect of the sole on foot. The study found that the diamond type porous sole had the best shock absorption.This study combines the finite element modeling method for the biomechanical analysis of the foot, provides a new research method for the cushioning design of the sole, has important reference value for the design of footwear and the customization of orthopedic footwear, and will greatly improve the design of functional shoes and the development of the shoe industry.



---

*Source: 2652137-2020-11-19.xml*