**A 2D Waveguide Method for Lithography Simulation of Thick SU-8 Photoresist**

#### **Zi-Chen Geng, Zai-Fa Zhou \* , Hui Dai and Qing-An Huang**

Key Laboratory of MEMS of the Ministry of Education, Southeast University, Nanjing 210096, China; 220184753@seu.edu.cn (Z.-C.G.); 220184865@seu.edu.cn (H.D.); hqa@seu.edu.cn (Q.-A.H.)

**\*** Correspondence: zfzhou@seu.edu.cn; Tel.: +86-025-8379-2632 (ext. 8817)

Received: 30 September 2020; Accepted: 27 October 2020; Published: 29 October 2020

**Abstract:** Due to the increasing complexity of microelectromechanical system (MEMS) devices, the accuracy and precision of two-dimensional microstructures of SU-8 negative thick photoresist have drawn more attention with the rapid development of UV lithography technology. This paper presents a high-precision lithography simulation model for thick SU-8 photoresist based on waveguide method to calculate light intensity in the photoresist and predict the profiles of developed SU-8 structures in two dimension. This method is based on rigorous electromagnetic field theory. The parameters that have significant influence on profile quality were studied. Using this model, the light intensity distribution was calculated, and the final resist morphology corresponding to the simulation results was examined. A series of simulations and experiments were conducted to verify the validity of the model. The simulation results were found to be in good agreement with the experimental results, and the simulation system demonstrated high accuracy and efficiency, with complex cases being efficiently handled.

**Keywords:** lithography simulation; microelectromechanical system; waveguide method

#### **1. Introduction**

#### *1.1. Research Background*

The lithography process refers to the process of transferring geometric pattern on the mask to the photoresist coated on the surface of the semiconductor wafer. With the rapid development of semiconductor technology, the lithography process has become an important propulsive force to achieve integrated and miniaturized devices. SU-8 photoresist [1] is a chemically amplified negative-tone resist based on epoxy resin, which has excellent chemical stability, resistance against solvents, and mechanical and optical properties [2]. It has become a mainstream micromachining technology for high aspect ratio structures in the field of microelectromechanical systems (MEMS) [3], including micrototal analysis systems [4], microfluidic systems [5], electroplating [6], and so on.

The propagation process of light during exposure is shown in Figure 1. First, we assume that the incident light is parallel, uniform, and perpendicular to the mask surface. The incident light is diffracted when passing through the mask hole. Then, when the light reaches the interface between the air and the photoresist, part of the light will be refracted and will continue to propagate through the photoresist, and the other part will be reflected back into the air. Part of the light that continues to propagate to the inside of the photoresist will be absorbed by the photoresist and generate photoacid at the same time. Finally, at the photoresist/substrate interface, part of the light is refracted into the substrate, and part of the light is reflected back into the photoresist and superimposed with the incident light. In general, the entire propagation process can be regarded as the propagation of diffracted light, taking into account other optical phenomena. Among them, the diffraction effect has the most important influence on light intensity distribution in the photoresist.

**Figure 1.** Light propagation process.

#### *1.2. Research Significance*

Because of expensive equipment and complicated process steps in the lithography process, we need repeated experiments, which is costly and time consuming. It is difficult to understand the inherent principle of lithography technology. Therefore, lithography simulation technology [7] has become a useful tool for predicting and optimizing the manufacturing process, which is also an integral part of lithography technology. Lithography simulation technology avoids the problems of high cost and time consumption caused by repeated plate making and tape production. At the same time, it can optimize the process parameters and promote the development of lithography technology. The flow chart of lithography simulation is as follows (Figure 2).

**Figure 2.** Flow chart of lithography simulation.

The final morphology of thick SU-8 photoresist after development depends largely on the exposure process [8]. The developed photoresist morphology can be predicted by simulating light intensity distribution in the photoresist. Therefore, it has become an inevitable and useful requirement to establish a lithography simulation model [9] of SU-8 thick photoresist to calculate light intensity distribution.

In this paper, a simulation model for UV lithography of SU-8 negative thick photoresist is presented to predict the profiles of developed SU-8 structures based on rigorous electromagnetic field theory. The model schematic is shown in Figure 3. We used the light intensity distribution calculation model to calculate the distribution of light intensity based on mask information, light source information, and photoresist information. The parameters of the air gap, incident angle, and photoresist depth, which have significant influence on the profile quality of developed structures, were studied in the simulation. Simulation results based on profiles of exposure intensity distribution were compared with the experimental results to validate the model.

**Figure 3.** Model schematic.

#### **2. Materials and Methods**

Light intensity simulations are based on two main categories of models. One category is based on scalar diffractive theory [10], and the other is based on rigorous electromagnetic field theory [11], including the finite difference time domain method (FDTD) and the waveguide (WG) method. The limitation of these models comes from simulation accuracy, application capability for special lithography methods such as inclined SU-8 photoresist lithography, and numerical costs for large-scale 3D simulations. Therefore, considering calculation accuracy and calculation time, we finally used the waveguide method to calculate light intensity distribution.

The waveguide [12,13] method has been developed and employed in lithography simulations during the past decade. The idea of the waveguide method is based on rigorous electromagnetic field theory, and it is well adapted to lithography simulation with many typical mask geometries. The first prototype of the waveguide method was presented by Nyyssonen [14] in order to simulate the measurement of optical linewidth. Then, this method was used to simulate diffraction when light passes through a 2D-phase shift mask based on Yuan's theory [15] and later extended to 3D domain [16] by Lucas.

The model is based on the 2D spatial frequency solution to Maxwell's equations, the so-called waveguide method [17–21]. Firstly, the simulated structure is sliced into different rectangular layers with homogeneous optical properties in the *Z*-axis direction, as shown in Figure 4. Then, the material parameters and electromagnetic fields in each layer are expanded by Fourier series. Simultaneous equations formed by these expanded series and Maxwell equations give an eigenvalue problem. After that, a hybrid approach using vector potentials is applied to couple electromagnetic field in each layer. Finally, each layer is coupled according to the proper boundary continuity condition, and the electromagnetic field information in the photoresist is converted into light intensity information.

**Figure 4.** The two-dimensional schematic of the waveguide method.

Firstly, the coupled Maxwell's equation determined by vector potential [22] can be simplified to two differential equations by separating the variable; the equation contains the complex potential *k* 2ε and the eigenvalue α 2 :

$$\begin{cases} \frac{d^2X}{dx^2} + k^2 \varepsilon X = \alpha^2 X\\ \frac{d^2Z}{dz^2} = -\alpha^2 Z \end{cases} \tag{1}$$

Before solving the equation numerically, the dielectric constant ε needs to be expanded by Fourier series. It is assumed that the mask structure is arranged infinitely and repeatedly, with the length *d* as the period in the x-direction:

$$\varepsilon^{j}(\mathbf{x}) = \sum\_{q=-L}^{+L} \varepsilon\_{q}^{j} \exp(i2\pi qb\mathbf{x}) \tag{2}$$

The coefficient of each term can be obtained by inverse Fourier transform:

$$\begin{aligned} \varepsilon\_q^j &= b \int\_0^d \varepsilon^j(\mathbf{x}) \exp(-i2\pi qb\mathbf{x}) d\mathbf{x} \\ X(\mathbf{x}) &= \sum\_l B\_l \exp(i2\pi l\mathbf{x}/d) \end{aligned} \tag{3}$$

Then, the material parameters and electromagnetic fields in each layer are expanded by Fourier series. Simultaneous equations formed by these expanded series and Maxwell equations give an eigenvalue problem:

$$\begin{aligned} DB &= \alpha^2 B \\ D\_{pq} &= k^2 \varepsilon\_{p-q} - (2\pi p/d)^2 \delta\_{p,q} \end{aligned} \tag{4}$$

where α and *B* are the eigenvalues and eigenvectors of a (2L + 1) by (2L + 1) matrix *D*. The computation time depends on the number of orders of the Fourier series expansion and the number of rectangular layers.

The electromagnetic field distribution in layer *j* is as follows. α *j <sup>m</sup>*, *B j <sup>m</sup>* is obtained by solving the eigenvalue problem, and *A j <sup>m</sup>*, *A* ′ *j <sup>m</sup>* is obtained using the boundary continuity condition:

$$\begin{cases} E\_y^i = \sum\_m \left[ A\_m^j \exp(i\,\alpha\_m^j(z - z\_j)) + A\_m^{j\} \exp(-i\,\alpha\_m^j(z - z\_j)) \right] \times \sum\_l B\_{l,m}^j \exp(i2\pi lx/d) \\\ H\_x^j = \frac{i}{k} \frac{\partial E\_y^i}{\partial z} \qquad H\_z^j = -\frac{i}{k} \frac{\partial E\_y^i}{\partial x} \end{cases} \tag{5}$$

After making this transformation and solving the eigenvalue problem, the boundary continuity conditions are matched so that the tangential components of the electric and magnetic fields have to be continuous at *z* = 0 [23]. The boundary continuity condition is as follows:

$$E\_y^j = E\_y^{j+1} \ H\_y^j = H\_y^{j+1} \tag{6}$$

In the air and first layer interface:

$$
\left[\begin{array}{c}
\mathbb{C}^{0}\_{11} \ \mathbb{C}^{0}\_{12} \\
0 \quad 0
\end{array}\right] \begin{bmatrix} A^{1} \\ A^{\prime 1} \end{bmatrix} = \begin{bmatrix} R \\ 0 \end{bmatrix} \tag{7}
$$

$$\begin{aligned} R\_l &= \mathbf{2} \left[ \mathbf{1} - \left( l\_0 b \lambda\_0 \right)^2 \right]^{\frac{1}{2}} \delta\_{l, l\_0} \\ &= \mathbf{2} \left[ \mathbf{1} - \left( \sin \theta \right)^2 \right]^{\frac{1}{2}} \delta\_{l, l\_0} \end{aligned} \tag{8}$$

where θ is the incident angle, and *l*<sup>0</sup> is the order of incident plane wave. It is convenient for us to change the different parameters to obtain the corresponding results.

By matching the boundary conditions at *z* = *T*, we obtain the following:

$$
\left[\begin{array}{c}\mathbb{C}^{n}\_{11}\,\mathrm{C}^{n}\_{12} \\ 0 \quad 0 \end{array}\right] \begin{array}{c} A^{n} \\ A^{\prime n} \end{array} \right] = \left[\begin{array}{c}0 \\ 0 \end{array}\right] \tag{9}
$$

At an intermediate interface *z* = *z<sup>j</sup>* , the waves in layer *j*+1 are related to the waves in layer *j* by

$$
\begin{bmatrix}
\mathbf{C}\_{11}^{j}\mathbf{C}\_{12}^{j} \\
\mathbf{C}\_{21}^{j}\mathbf{C}\_{22}^{j}
\end{bmatrix}
\begin{bmatrix}
A^{j} \\
A^{\prime j}
\end{bmatrix} = \begin{bmatrix}
E\_{11}^{j+1}\mathbf{E}\_{12}^{j+1} \\
E\_{21}^{j+1}\mathbf{E}\_{22}^{j+1}
\end{bmatrix}
\begin{bmatrix}
A^{j+1} \\
A^{\prime j+1}
\end{bmatrix} \tag{10}
$$

where matrix *C* contains only the topography information, such as mask size, material parameters, and so on, and the vector [*A*, *A* ′ ] contains diffraction information.

Finally, the light intensity can be calculated by the electric field value:

$$I(p,q) = n\_r \left| \mathbf{E}\_y^j(p,q) \right|^2 \tag{11}$$

The incident angle is due to refraction at the air/photoresist interface. The structural inclined angle and the structural width of the developed photoresist vary with the incident angle. The relationship between the incident angle δ, the refractive angle θ, and the structural inclined angle α is determined by Snell's law as follows:

$$
\Theta = \sin^{-1}(n\_1 \cdot \sin \delta / n\_2) \tag{12}
$$

$$
\alpha = 90 - \theta \tag{13}
$$

where *n<sup>1</sup>* (= 1) and *n<sup>2</sup>* (= 1.67) are the refractive indices of the air and SU-8, respectively.

More details of the method is described in reference [24–26].

#### **3. Results and Discussion**

#### *3.1. Simulation Results*

The simulation results showed the distribution of normalized intensity calculated under different conditions. In addition to finding light intensity distribution on a certain layer of photoresist, a contour map can also describe the intensity value of a two-dimensional profile of photoresist. This method can help intuitively predict the final development morphology of the photoresist. Observing the influence of different parameters on light intensity distribution provides a guide to set the corresponding process parameters so as to optimize the structure design.

In the exposure process, the surface photolithographic dose *D* (mJ/cm<sup>2</sup> ) is multiplied by the exposure time and the incident light intensity. The surface exposure dose can be considered to be constant over the entire thickness of the photoresist layer. Photolithography of thick-film photoresist implies that the thickness of the photoresist material needs to be taken into account in order to determine the optimal exposure dose for a given photoresist film thickness. The exposure dose is the light intensity at a given thickness and exposure time. According to the material parameters of the photoresist and substrate, the depth of the photoresist can be calculated to obtain the optimal exposure dose under the given photoresist depth. In this calculation, the incident light intensity was 2.6 mW/cm<sup>2</sup> , its wavelength was 365 nm, and the thickness of SU-8 photoresist was 300 µm. According to our simulation conditions, when the thickness of SU-8 photoresist is 300 µm, the optimal exposure dose is 140 (mJ/cm<sup>2</sup> ) [27].

First, we analyzed the simulation of the intensity of light at vertical incidence. The intensity curves for vertical exposure with different thicknesses are shown in Figure 5, and the corresponding contour maps are shown in Figure 6. We can see from the results that the intensity in SU-8 was simultaneously attenuated along the radiation direction with increasing photoresist thickness because of the diffraction and absorption of the incident light in the photoresist.

**Figure 5.** The intensity curve for vertical exposure with different depths.

**Figure 6.** The contour map of intensity distribution for vertical exposure.

The previous simulation was based on the case of a single mask hole. In practice, there are often multiple mask holes, and selection of the spacing between the mask holes is an important parameter to be investigated when designing the layout. Different mask hole spacing will affect the final development morphology of the photoresist. Therefore, Figure 7 shows the intensity distribution of the photoresist in the case of multiple mask holes, and the corresponding contour maps are shown in Figure 8. The mask hole spacing was 40 µm, and the mask hole size was 70 µm. The length of the photoresist was 200 µm, and the abscissa in the figure represents the horizontal position. According to the proportion calculation, it can be seen that the simulation results are correct.

**Figure 7.** The intensity curve for vertical exposure with different depths.

**Figure 8.** The contour map of intensity distribution for vertical exposure.

The above simulation results were obtained without considering the substrate reflection. Part of the light is reflected at the photoresist and substrate interfaces. The refractive index can be calculated according to the refractive index of the photoresist and the substrate, thereby enabling calculation of the distribution of light intensity after refraction. During vertical UV lithography processes of SU-8, reflection enhances light intensity at the same position as the photoresist. During inclined UV lithography processes of SU-8, the incident UV light is reflected at the SU-8/substrate interface. The reflected UV light will be transmitted into other unexposed SU-8, usually creating reflected induced structures. The light intensity distribution under different conditions at the same position is shown in Figures 9 and 10. By comparing the value of light intensity with and without reflection at the same position, we found that the light intensity considering substrate reflection was greater than that without substrate reflection, and this phenomenon was more obvious at the bottom of the photoresist. This is consistent with the previous analysis.

**Figure 9.** The intensity curve for vertical exposure at depth of 150 µm.

**Figure 10.** The intensity curve for vertical exposure at depth of 300 µm.

One of the advantages of substrate reflection is that it can ensure the photoresist bottom is fully exposed, avoiding the disadvantage of the photoresist bottom not being fully exposed due to it being too thick. It should be noted that if the substrate reflection is not well controlled, it will lead to overexposure of the photoresist bottom, which in turn will lead to deviation between the results and the design. The reflection at the photoresist/substrate interface cannot be neglected as the dose of reflected UV light is increased by increasing the exposure time to more than that of a regular exposure process. To eliminate the reflected induced patterns, we can coat a layer of antireflective film on the photoresist or substrate surface, called TARC and BARC, respectively.

Secondly, we analyzed the simulation of light intensity at oblique incidence. The incident UV light was inclined from left to right, and the incident light intensity was 2.6 mW/cm<sup>2</sup> . Other conditions were the same as vertical incidence.

The intensity curves for inclined exposure with different air gaps are shown in Figure 11, and the corresponding contour maps are shown in Figure 12. As can be seen, as the air gap increased from 10 to 30 µm, the exposure area shifted to the right, the diffraction effect became more intense, and the edge diffraction effect also became severe.

**Figure 11.** The intensity curve for inclined exposure with different air gaps (10, 20, and 30 µm) at 23.5◦ incident angle.

**Figure 12.** The contour map of intensity with different air gaps (10, 20, and 30 µm) at 23.5◦ incident angle.

Figure 13 illustrates the intensity curves resulting from different incident angles, and the corresponding contour maps are shown in Figure 14. We can see from the results that the exposure area shifted to the right, and the light intensity began to decrease with increasing inclined angle. Analyzing this phenomenon from a mathematical point of view, according to the expression of matrix *R* in Formula (8), the value of *R* decreased with the increase in incident angle, and therefore the final light intensity also decreased.

**Figure 13.** The intensity curve for inclined exposure with different incident angles in the X–Z cross section.

**Figure 14.** The contour map of intensity with different incident angles (15, 23.5, and 30◦ ) with air gap of 10 µm.

The intensity curves for inclined exposure with different thicknesses are shown in Figure 15, and the corresponding contour maps are shown in Figure 16. It can be seen from the Figure 15 that due to diffraction of light and absorption of the photoresist, the light intensity was different at different depths. On the surface of the photoresist, the light intensity was the highest, while it was the smallest on the bottom of the photoresist. It should be noted that the intensity in SU-8 was simultaneously attenuated along the radiation direction with increasing photoresist thickness.

**Figure 15.** The intensity curve for inclined exposure with different depths at 23.5◦ incident angle.

**Figure 16.** The contour map of intensity distribution for inclined exposure at incident angle 23.5◦ .

Compared with vertical incidence, the intensity curve under oblique incidence decreased from top to bottom and shifted to the right. Further simulation showed that this translation amplitude would be larger with the increase of incident angle. Finally, the computation time depended on the number of orders of the Fourier series expansion. With an increase of the expansion order, the calculation accuracy was higher and the calculation time increased accordingly. Therefore, in actual simulation, it is necessary to select the appropriate expansion order in order to consider the calculation time and calculation accuracy.

#### *3.2. Experimental Results and Analysis*

Experiments for UV lithography were conducted to validate the proposed lithography simulation model, and simulation results based on exposure intensity distribution patterns are illustrated for demonstration. In this calculation, the incident light intensity was 2.6 mW/cm<sup>2</sup> , and its wavelength was 365 nm. The thickness of SU-8 photoresist was 300 µm, and the thickness of the air gap was 10 µm. Besides, in order to efficiently evaluate and compare the simulation and experimental results, the contour value of light intensity patterns was assumed as 0.38 in this study.

SEM photos of SU-8 structures fabricated on glass wafers are illustrated in Figure 17a, and the corresponding simulation results are shown in Figure 17b. From Figure 17a, we can see that the width of the top SU-8 was 31.37 µm, and the width of the bottom SU-8 was 32.06 µm. Therefore, after development, the width of SU-8 pillars was enlarged over the entire height of the photoresist. In particular, the variation of the bottom width was 0.76 µm. Due to the absorption of light by photoresist, for a relatively low exposure dose near the bottom of the SU-8 photoresist layer, the width at the bottom is slightly larger than that at the top. By means of threshold prediction (the contour value was assumed as 0.38) of the light intensity pattern, the top and bottom widths derived from the simulation result in Figure 17b were 33.28 and 34.7 µm, respectively. In comparison, the top width in the intensity pattern was approximately equal to that in the experimental result. Therefore, we drew the conclusion that the simulation results were in good agreement with the experimental results.

**Figure 17.** SEM photo of a SU-8 pillar (**a**) and simulation result (**b**) for mask hole size of 30 µm.

Figures 18 and 19 show a series of SU-8 pillars with different proportions of linewidth and spacewidth and the corresponding simulation results of exposure intensity. SEM photos of the SU-8 structures fabricated on glass wafers are illustrated in Figures 18a and 19a, and the corresponding simulation results are shown in Figures 18a and 19b. In practice, there are often multiple mask holes, and the spacing between the mask holes is an important parameter. Different mask hole spacing will affect the final development morphology of photoresist.

**Figure 18.** SEM photo of a SU-8 pillar (**a**) and simulation result (**b**) for linewidth/spacewidth of 30/30 µm.

Figure 18 shows the simulation and experimental results for linewidth/spacewidth of 30/30 µm. Through measurement and calculation, the top and bottom widths of the photoresist cylinder, shown in Figure 18a, were found to be 32.35 and 35.29 µm, respectively. Therefore, the top and bottom widths of developed SU-8 pillars were enlarged by 2.35 and 5.29 µm, respectively. Figure 19 shows the simulation and experimental results for linewidth/spacewidth of 30/50 µm. The top and bottom widths of the photoresist cylinder, shown in Figure 19a, were 31.78 and 33.08 µm, respectively. Therefore, the top and bottom widths of developed SU-8 pillars were enlarged by 1.78 and 3.08 µm, respectively. The width of the bottom of SU-8 glue is larger than that of the top. This is because the polymer is not very dense when the exposure dose at the bottom of SU-8 photoresist is relatively low. To compare the experimental and simulation results, the dimensions of simulation results were also measured according to the contour value of light intensity patterns. The top and bottom widths for linewidth/spacewidth of 30/30 µm, shown in Figure 18b, were 31.61 and 30.58 µm, respectively, while the top and bottom widths for linewidth/spacewidth of 30/50 µm, shown in Figure 19b, were 31.43 and 30.68 µm, respectively. A comparison between Figure 18a,b and Figure 19a,b showed that the simulation results exhibited a similar trend as the experimental results. However, the disparity between the experimental and simulation results in Figure 18 is notable. This was due to the decrease of mask hole spacing. The results revealed that the interaction between neighboring slits became more significant with diminution of the interval between two slits, and the diffraction effect was enhanced with a decrease in the interval between two slits, leading to the linewidth of the SU-8 structure with linewidth/spacewidth of 30/30 µm becoming a little larger than that with linewidth/spacewidth of 30/50 µm.

**Figure 19.** SEM photo of a SU-8 pillar (**a**) and simulation result (**b**) for linewidth/spacewidth of 30/50 µm.

SEM photos of the oblique SU-8 structures fabricated on glass wafers are illustrated in Figure 20a, and the corresponding simulation results are shown in Figure 20b. Figure 20a,b shows the intensity distribution of photoresist in the case of multiple mask holes with the incident angle of 30◦ and linewidth/spacewidth of 30/60 µm. From Figure 20a, we can see that the width of the top SU-8 was 31.75 µm, and the width of the bottom SU-8 was 33.40 µm. Therefore, the top and bottom widths of developed SU-8 pillars were enlarged by 1.75 and 3.40 µm, respectively. In particular, the variation of the bottom width was 1.75 µm. We also calculated the top and bottom widths in the simulation

results. In comparison, the top width in the intensity pattern was approximately equal to that in the experimental result. Therefore we drew the conclusion that the simulation results were in good agreement with the experimental results. The measured structural inclined angles of developed SU-8 photoresists were 73.76◦ for the incident angle of 30◦ . Correlating Formula (12), the theoretical structural inclined angle α of developed structures was 72.58◦ for the incident angle of 30◦ , which was approximately the same as the measured value and thus demonstrated good agreement with the experimental results.

**Figure 20.** SEM photo of an oblique SU-8 pillar (**a**) and simulation result (**b**) for incident angle of 30◦ and linewidth/spacewidth of 30/60 µm.

#### **4. Conclusions**

In order to simulate UV lithography and predict the profiles of SU-8 structures, a lithography simulation model based on the waveguide method and rigorous electromagnetic field theory is presented in this paper. According to this model, vertical and inclined UV exposure of SU-8 thick photoresist was successfully simulated, and the parameters significantly influencing light intensity distribution were studied. The presented model using the WG method works faster than other electromagnetic field methods and produces precise results. A series of experiments were designed to investigate the performance of the simulation model. By comparing simulation and experiments results, the method was successfully verified. This model is ideally suited for 2D lithography simulations for complex cases.

**Author Contributions:** Z.-C.G. and Z.-F.Z. developed the model and designed the experiments; Z.-C.G. performed the experiments; H.D. analyzed the data; Q.-A.H. proposed the idea for the model and made suggestions for improvement. Z.-C.G. and Z.-F.Z. edited this manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was supported by the National Key R&D Program of China under grant number 2018YFB2002600.

**Conflicts of Interest:** The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; and in the decision to publish the results.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

**Euan Langford <sup>1</sup> , Christian Andrew Griffiths 2,\* , Andrew Rees <sup>2</sup> and Josh Bird <sup>3</sup>**


**\*** Correspondence: c.a.griffiths@swansea.ac.uk

**Abstract:** This paper studies the forces acting upon the Intraosseous Transcutaneous Amputation Prosthesis, ITAP, that has been designed for use in a quarter amputated femur. To design in a failure feature, utilising a safety notch, which would stop excessive stress, σ, permeating the bone causing damage to the user. To achieve this, the topology of the ITAP was studied using MATLAB and ANSYS models with a wide range of component volumes. The topology analysis identified critical materials and local maximum stresses when modelling the applied loads. This together with additive layer manufacture allows for bespoke prosthetics that can improve patient outcomes. Further research is needed to design a fully functional, failure feature that is operational when extreme loads are applied from any direction. Physical testing is needed for validation of this study. Further research is also recommended on the design so that the σ within the ITAP is less than the yield stress, σs, of bone when other loads are applied from running and other activities.

**Keywords:** prostheses; ITAP; micro topology; ANSYS; MATLAB; additive manufacture

#### **1. Introduction**

The failure features of the femur are unique to the human skeleton because the bone is designed to transfer σ throughout the structure when compressive or tensile loads are applied. The σ is also dispersed to the surrounding muscle, cartilage and connected bone in the skeleton. This allows the femur to handle large σ when the user is in motion such as running and jumping. However, the femur is not indestructible as the bone structure is susceptible to fracture from extreme shear loads. It has been theorised by that a sheer force of 4000 N applied perpendicular to the bone would cause a fracture [1]. The properties of the individual's bone structure are also key factors for the propagation of fracturing as the physiology is unique to the individual. This means that some individuals and ages are more susceptible to fractures particularly when bone decreases after the age of 30 [2]. The femur can also fail under different extreme loads from, for example, collision or devices such as landmines, improvised explosive devices, IEDs, and other explosives as identified by Mckay et al. [3]. The direction and position of the applied σ is difficult to define, because the load could be applied from any direction at any time.

The intraosseous transcutaneous amputation prosthesis, ITAP, is a direct skeleton attachment that is used by amputees [4]. This is used rather than a standard socket interface that causes swelling, discomfort and potential infection for the user. The socket interface often prohibits the life span of the use of the prosthetic [5]. As a result, the development of the ITAP and similar direct skeleton attachments have been designed to mitigate these problems, to improve the quality of life for the user [4]. ITAPs were initially designed from studying animals with protruding bone such as deer and rhinos [6]. The mechanics and biomechanics of these animals' bone development have been analysed due to the lack of infection of the tissue surrounding the protruding bone which is similar to an embedded prosthetic. This is vital, as an embedded prosthetic can be susceptible to infection and

**Citation:** Langford, E.; Griffiths, C.A.; Rees, A.; Bird, J. The Micro Topology and Statistical Analysis of the Forces of Walking and Failure of an ITAP in a Femur. *Micromachines* **2021**, *12*, 298. https://doi.org/ 10.3390/mi12030298

Academic Editors: Davide Masato and Giovanni Lucchetta

Received: 17 February 2021 Accepted: 11 March 2021 Published: 12 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

damage of the local tissue and if a prosthetic is not properly manufactured with biomedical procedures with appropriate postproduction methods [4]. Consequently, the ITAP is constructed in a manner that replicates the antlers of deer. An open insertion to the bone with reduced bacterial colonization infection was realised using a diamond-like carbon (DLC) ITAP [6].

The biomechanics of the ITAP and the human body need to be thought of during design, simulation and production as the material must not be harmful to the patient's body. As such, the main material used is that of Titanium (Ti). This is because Ti is biocompatible with the human body, with a high processability through Computer-aided design (CAD) and additive manufacturing (AM) processes. Ti is also strong enough to handle the loads associated with walking and other movements [7,8]. However, it can deteriorate in to smaller fragments in the body over time due to movement. This is most commonly attributed to hip and knee implants because there is constant friction between the implant and the joint [9]. Additional alternative materials can be used for prosthetics such as Polyetheretherketone (Peek). This can be processed by AM methods and inserted into the bone to help repair and supplement the bone after injury. This is commonly used for American football players for bone repair therapy as Peek is resistant to infection [10]. However, Peek is not suitable for use in a load bearing ITAP as there is poor bonding between the bone and the material. Therefore, it could easily become possible for the ITAP to be partially or fully dislodged from the bone, causing internal damage to the user.

The ITAP can be used for both cosmetic and load bearing prosthetics; the former can range from prosthetic eyes, fingers, noses, ears amongst others [11], while the latter can be used for prosthetic limbs with loadbearing forces acting upon extremities. This includes, but is not limited to, both high and low leg amputations, either above or below the knee. However, the work of Newcombe et al., identified that an ITAP implant was not advisable for an amputation greater than one quarter of the original length of the bone. This was due to the residual σ in the femur housing, known as the anchor, being too great so damage to the cortical bone could occur [12]. Further suggested amendments to the ITAP have been developed to reduces the σ with in both the bone and the ITAP. This includes a safety notch that is positioned outside the user's congealed tissue. The proposed modification could prevent an extreme load from damaging the user as the ITAP would failure at the safety notch [13].

Other direct skeletal prosthetics are available for use such as Osseointegrated Prostheses for the Rehabilitation of Amputees, OPRA. This system involves screwing into the bone cavity which increases the residual σ all along the anchor. This system works because it increases the user's hip flexibility and extension in comparison to traditional socket interfaces. However, over time it was recorded that pain was detected due to the build-up of residual σ [14]. The ITAP is not only suited to use in humans but it has had clinical trial in canines. Noel Fitzpatrick et al. studied the procedure for installation and use in four separate canines. The study identified that all subjects' quality of life improved due to the prosthetic attachment. However, canines 1–3 were euthanized due to metastatic disease spreading throughout the body; it was unclear if this was related to the ITAP [15]. As the ITAP can be fitted for different purposes and species, bespoke parts can be produced. This is extremely useful as the properties of bone vary as well as the length of the remaining bone after amputation, see Figures 1 and 2. Bone shape and modelling can be recorded by using Computed Tomography (CT) scanning. This creates a need for design iterations to model and build an ITAP that suits the individual's lifestyle. In this process it is essential to design a build that does not over engineer the forces acting upon the leg. Failure to work in this design constraint could cause near fatal damage to the user resulting in further amputation or loss of life. This is due to the failure of the anchor supporting the ITAP and the potential rupturing of arteries and veins surrounding the localized area of the prosthetic. The work of Sullivan et al. concluded that the use of direct skeleton attachments for prosthetics increases the user's quality of life and the function of the prosthetic limb. Comfort also increased due to the lack of sores developing on the stump. This is because

there is no contact between the stump of the limb and the prosthetic sock fitted over the top. However, the reports concluded that aspects of the development of the ITAP needed further advancement because the "trans-femoral amputee" case needs extensive review before it becomes a common clinical procedure for prosthetic rehabilitation [16]. It has been identified by Bird et al., that the initial area of failure and high σ concentration of the ITAP is located at the root of the ITAP [13]. The results of the study identified that a single safety notch was not enough for the ITAP to be suitable for use. This is because the strength of the Ti used for the ITAP is greater than the bone, so the bone would fail before the ITAP, resulting in damage to the user. Thus, it is identified that a comparable σ between bone and the ITAP at the point of insertion is needed for successful development of the product. This is categorized as designing in failure for success. needed further advancement because the "trans femoral amputee" case needs extensive , that the initial area of failure and high σ concentration σ between bone and the ITAP at the point of insertion is needed fo needed further advancement because the "trans femoral amputee" case needs extensive , that the initial area of failure and high σ concentration σ between bone and the ITAP at the point of insertion is needed fo

**Figure 2.** (**A**) Canine prosthetic ITAP embedded into its front left leg. (**B**) A Canine bespoke ITAP embedded into its front right leg [15].

3). Both walking and extreme σ studies are considered because the operational σ of the design need 3). Both walking and extreme σ studies are considered because the operational σ of the design need The benefit of using AM methods of production is that the topology of the ITAP can be adapted to meet the users bespoke requirement. The manufacturing process allows for modification of the internal structure, of the ITAP and each can be printed to have the same properties and micro topology as the bone being amputated, as seen via CT scanning. The first aim of this study is to identify the optimum topology of the ITAP, locating critical material that is integral to the structure of the design when simulated to both walking and extreme loads. The secondary aim is to optimise the design of the ITAP topology so that a controlled failure can be achieved when extreme loads are applied. Thus, producing an ITAP that prevents a further damage to bone when subjected to excessive loads. For both aims the simulation will consider the design of a quarter amputation of a femur with an ITAP embedded into the bone (Figure 3). Both walking and extreme σ studies are considered because the operational σ of the design needs to be evaluated in order to design a safe prosthetic that can be embedded into the bone. To ensure reliable results a multiple modelling software approach is used to generate the topologies. The paper is organised as

follows: the next section uses simulation methods to identify the maximum stress, σMax, on the standard ITAP due to walking loads. This is followed by topology optimisation of the ITAP model using ANSYS and MATLAB. Then, Section 4. considers the influence of the topology of the ITAP design when under extreme loads. Finally, conclusions are presented on the optimisation of ITAP designs. stress, σ

**Figure 3.** (**A**) Quarter Amputated Femur with accompanying ITAP CAD models, (**B**) Femur and ITAP assembly [17].

#### **2. Experimental Methodology**

focus will be on performing Simulations of σ The purpose of the research is to establish the strength of the ITAP, and to identify the behaviour of a safety mechanism to prevent excessive forces transferring to bone. The focus will be on performing Simulations of σMax on the ITAP due to walking loads. In order to perform this the following section will identify initial design, boundary conditions and mesh used. In 2.2 the Simulated results of the ITAP standard design will be presented and then used as a benchmark for the ITAP designs with topology optimisation in mboxsectsect:sec3-micromachines-1132465.

#### *2.1. The Initial Design, Boundary Conditions and Mesh*

The A quarter amputation of the femur requires an ITAP with an embedded depth of 160 mm, with a diameter of 14 mm [12]. The femur provided in this study was obtained by CT scanning of a deceased male aged 44 weighing at 85 kg and at a height of 185 cm [17]. Based on this femur the standard ITAP design can be seen in Figure 3 The ITAP has been segmented into 4 section as to identify the areas of the recommended change of the topology, when simulated with the loads associated with walking and bone failure. This has been done as it is imperative that the ITAP bonds to the bone and allows the muscle to congeal around the protruding end. As such the changing topology must be correctly identified and implemented, in different section to reduce the risk of failure of insulation. As the designed in failure mechanic must failure externally of the body. Thus, the segments utilized are Section A, B, C and D theses represent the region of the ITAP: embedded in to the femur that must have a constant surface contact with bone; located at the end of the users stump where the collecting skin, fat and muscle fusses; were the desired failure mechanic is needed to be implemented for a safe failure of the ITAP; were the prosthetic leg is attached to the user, respectively. For the simulation of the design the following assumptions are mad.


The mesh was developed by conducting a mesh sensitivity study. Ten simulations were carried out with different mesh refinements to determine the following σ results:


From the study the optimum controls for the simulation where the following;


This process removed errors and many delocalisations within the mesh to ensure each section is uniform and in alignment with each other. However, when reviewing the full geometry and specific section minor desolations where identified. This is due to the simulated model being constructed from separate sections, forming a full geometry, rather than one single section. This purpose of this was to allow an in-depth review of each section of the ITAP and change the failure mechanic within section C of the ITAP.

The simulation boundary locations have been set at the femur joint with the pelvis. Frontal (Fx) Lateral (Fy) and Axial forces (Fz) are acting on the ITAP in the region where the prosthetic would be applied as this is the same area as the patellar surface where the knee would be attached [19], as seen in Figure 4.

The work of Georg et al. identified that the forces acting on the femur are dependent on the location of interest within the bone, the body weight (BW) of the user and finally the position of the motion of walking [21]. As such, the Fx, F<sup>y</sup> and F<sup>z</sup> maximum forces are calculated as 958.9 N, −833.8 N and 3168.6 N, respectively.

along the bone's structure [21], **Figure 4.** (**A**) Forces acting upon the femur along the bone's structure [21], (**B**) FBD of the ITAP, built into Sections A, B, C and D.

#### *2.2. Simulated Results of the ITAP Standard Design*

The σ σ increases. It can also be identified that the location of the σ y's core and only Section C's core, this is partly due to minor delocalisation in the mesh within the be identified that the recorded σ in Section C's core for all simu exceed the σ The σMax recorded when the ITAP is subjected to walking loads is shown in Table 1. It can be identified that as the diameter of the safety notch decreases from 14 mm to 5 mm the σMax increases. It can also be identified that the location of the σMax within the ITAP changes its position as the size of the notch diameter decreases, moving closer to the notch from the point of insertion. Discrepancies arise when analysing the full geometry's core and only Section C's core, this is partly due to minor delocalisation in the mesh within the model, further improvements have been utilised for additional readings. However, it can be identified that the recorded σ in Section C's core for all simulated models does not exceed the σ<sup>s</sup> of Ti. Further conclusions can be drawn upon from this preliminary study, these being:


The results indicate that the ITAP design needs further optimisation to reduce the σ at the joint of bone and ITAP to less than the σ of bone, as well reducing the σ across the whole design so that the σ is less than the σ as the ITAP must fail within the core of Section C when the σ generated by external loads life such as walking, as seen in a 5 mm safety notch were the recorded σ is at 3313.8 MPa. The results indicate that the ITAP design needs further optimisation to reduce the σ at the joint of bone and ITAP to less than the σ<sup>s</sup> of bone, as well reducing the σ across the whole design so that the σMax is less than the σ<sup>s</sup> of Ti. Topology analyses shall be utilised as the ITAP must fail within the core of Section C when the σ generated by external loads are greater than the limits of the system, thereby avoiding irreparable damage to the user. However, the core and ITAP must not fail under standard forces generated in day to day life such as walking, as seen in a 5 mm safety notch were the recorded σ is at 3313.8 MPa.


**Table 1.** Simulated σMax of the ITAP designs.

#### **3. Topology Optimisation of for the ITAP Model**

In the following section the approach to using topology modifications to the ITAP will be show. Firstly, the ANSYS model will be described followed by the 2D and 3D approaches using MATLAB functions.

#### *3.1. ANSYS Model*

With a material reduction (50%) the design of the standard ITAP and notched designs have been simulated in ANSYS. In total 20 different simulations where conducted utilising the topology reduction calculations in order to identify key structural information for the IP when underload. The model mass is reduced in areas that are not critical to the mechanical structure, as seen in Figures 5–7. The same forces identified in Section 2.2 are used and the fixed supports used are as follows. –


**Figure 5.** Fixed contact between the bone-anchor and ITAP, (–**A**–**D**) Standard ITAP to 5 mm safety notch topology reduction models results, (**E**) FBD of the applied walking loads.

– **Figure 6.** Fixed contact for the head of the bone-anchor to ITAP, (**A**–**D**) Standard ITAP to 5 mm safety notch topology reduction models results, (**E**) FBD of the applied walking loads. – –

– – **Figure 7.** Fixed contact for the head of Sections B and C of the ITAP, (**A**–**D**) Standard ITAP to 5 mm safety notch topology reduction models results, (**E**) FBD of the applied walking loads. –

σ distribution in the 5 mm safety notched ITAP σ distribution in the 5 mm safety notched ITAP σ distribution in the 5 mm safety notched ITAP **Figure 8.** (**A**) Initial ITAP designed, (**B**) <sup>σ</sup> distribution in the 5 mm safety notched ITAP.

#### *3.2. 2D Simulations Using MATLAB Functions*

In this section a 2D infill for the entirety of the interior of the ITAP geometry is developed. The coding for designing an optimised interior where material can theoretically

be removed originated as a 2D MATLAB problem written by Otomori et al. [22]. This has since been utilised to identify if the ANSYS topology studies identify the correct locations to remove non-critical material from the cylindrical cross-section of the ITAP's external Sections C and D. The volume of material remaining is denoted by Vmax, and is the only variable in these simulations. Different iterations were completed looking at Vmax (50–95% in intervals of 5%), in total 10 different simulations were conducted to complete this section of research. Using MATLAB code for a level set-based topology optimization the assumptions made in the model are as follows; section of the ITAP's external –


The resulting design structures (Figure 9) show that as the maximum volume of the model increases, the size of the voids forming in the ITAP decrease. This is most prominent in the voids forming at the boundary wall and the voids forming where sections C and D meet. This study has identified where the areas of the topology of the ITAP should be changed depending on the volume of the material used in additive manufactured production processes for complex internal structures. The topology simulation using MATLABs' applied load only consisted of one shear force, rather than the three forces that are applied to the ITAP in motion. Furthermore, the force cannot be set to the applied loads identified in Section 2.1. As such, the results from this study can only be used as a guideline. ' MATLABs' applied load only consisted of one shear force, rather than the three forces that

**Figure 9.** Simulated results of the 2D MATLAB topology simulations identifying the theoretical areas of material reduction of 50 to 95%.

#### *3.3. 3D Simulations Using MATLAB Functions*

To identify the influence of multiple forces on ITAP topology further simulations were needed. The work of Liu et al. allows for simple 3D models to be designed and simulated, where fixed support and applied loads can be controlled with multiple forces being utilised at a single location [23]. Using cubic models three separate simulation models have been studied looking at Sections C and D of the ITAP. These consist of a 85 mm by 14 mm by 14 mm model (Standard), as well as two models with cross-sections through the centre of the model at 85 mm by 14 mm by 1 mm (Vertical) and 85 mm by 1 mm by 14 mm (Horizontal). A 3D multsicale topology optimization code for lattice microscale response is used (top3d, Nelx, Nely, Nelz, Volfrac, Penal, Rmin) and the assumptions made for this study are as follows:

	- # Standard: X1 = 0, X2 = 0, Y1 = 0, Y2 = 14, Z1 = 0, Z2 = 14.
	- # Vertical: X1 = 0, X2 = 0, Y1 = 0, Y2 = 14, Z1 = 0, Z2 = 1.
	- # Horizontal: X1 = 0, X2 = 0, Y1 = 0, Y2 = 1, Z1 = 0, Z2 = 14.
	- # Standard: Fx, Fy and Fz at X1 = 85, X2 = 85, Y1 = 0, Y2 = 14, Z1 = 0, Z2 = 14.
	- # Vertical: Fx and Fy at X1 = 85, X2 = 85, Y1 = 0, Y2 = 14, Z1 = 0, Z2 = 1.
	- # Horizontal: Fx and Fz at X1 = 85, X2 = 85, Y1 = 0, Y2 = 1, Z1 = 0, Z2 = 14.

The MATLAB 3D simulations provide the most detailed results. From the study, the areas where material can be removed to save weight due to the minimal σ present is consistently located at two positions. These being at the centroid of the fixed support at the boundary wall where the ITAP's Sections B and C meet. The next area where material can be removed is located within Section D, 57 mm < x < 85 mm see Figure 10A, where material is being removed on the opposite side of the applied Lateral and Frontal loads. These two positions of material reduction are concordant with both the ANSYS and 2D MATLAB simulations as well as all three models of the 3D MATLAB simulations. The size of the voids forming is dependent on the volume percentage of material being simulated. The results from the Standard models identify how the ITAP could be shaped to meet the individual's body characteristics using the loads calculated. The results identify the areas of material that are less important to the structural integrity of the external sections of the ITAP. The Vertical and Horizontal simulations look at the cross-section of the ITAP with the Axial and either Lateral or Frontal forces applied, Furthermore, the remaining internal structure mimics the internal structure found in bone, Figure 10, forming pathways in-between voids in the ITAP, similar to the findings of both Otomori and Wu with their

respective studies and coding [22,24]. Both the Vertical and Horizontal results show this pattern. The voids are geometrically triangular in design as they are extremely capable of distributing σ within their shape [24]. There is minimal variation between the two sets of results, and this is expected as the forces being applied are similar in size. As the percentage of the material used increases, the definition of the voids and interconnecting pathways decreases. For these study 30 simulations were conducted utilising 10 different volume fraction limits for three different tested geometries.

The visual results of the simulations can be seen in Figure 10 for selected volume percentages of material.

**Figure 10.** (**A**) FBD of the 3D MATLAB topology simulation utilising the external Sections, C and D, of the ITAP, dimensions of 85 mm by 14 mm by 14 mm with applied loads, (**B**) Vertical model using Fx and Fy 95% volume of material, (**C**) Vertical model using Fx and Fy 55% volume of material.

#### **4. Topology of the ITAP Design under Extreme Loads**

areas where material can be removed to save weight due to the minimal σ present is conboundary wall where the ITAP's Sections B and C meet. The next area where material can dividual's body characteristics using the loads calculated. The results identify the areas of In this study the shear force being applied is set to the interconnecting surface between the prosthetic and the ITAP, represented by F<sup>y</sup> as −4000 N. There is a secondary axial force of 416.9 N applied, denoted by Fz. This has been derived from half of the user's weight resting on the ITAP, from the 85 kg test individual. From analysing the results of the simulation, the location of the σMax can be seen to be at the joint between Section A and B, where the ITAP leaves contact with the bone anchor, as seen in Figure 8B. This is the same location as in the walking simulations. The only difference is the position on the circumference of the joint of Sections A and B, as the forces are acting in different directions. The σMax has been recorded above the σ<sup>s</sup> of Ti and both the compressive and tensile σ<sup>s</sup> of bone; therefore, a fracture has a high chance of occurring under this load. From these findings, it is recommended that further study is needed into the design of the ITAP to remove the risk of critical failure of the bone when exposed to extreme loads.

#### *4.1. ANSYS Topology of the ITAP with Extreme Loads Applied*

For this study the assumptions and controls are the same as in Section 3.1. The results identify that voids are forming in the same locations as in the walking simulations, these

distributing σ within their shape [24]. There is minimal variation between the two sets of

being at the boundary walls between sections and fixed supports. Material is also being removed along the sides of the ITAP in the same manner as the walking simulations. However, the main differences in the results is the voids forming within Section D, where more material remains in these locations to disperse the σ throughout the model from the applied loads. As such 4-separate standard ITAP models where simulated utilising ANSYS's high computational topology modelling process. As seen in Figures 11 and 12. more material remains in these locations to disperse the σ throughout the model from the SYS's high computational topology mo

−

ed from half of the user's

the simulation, the location of the σ

σ

tions. The σ above the σ

**Figure 11.** Topology simulations of the full standard ITAP with extreme load F, which is the amalgamation of F<sup>y</sup> and F<sup>z</sup> at 50% volume of material, (**A**) FBD of the ITAP in Sections, (**B**) The resulting topology of the individual Sections, (**C**) FBD of the ITAP constructed in one Section, (**D**) The resulting topology of the full ITAP.

**Figure 12.** Topology simulations of Sections C and D of the standard ITAP with extreme load F, which is the amalgamation of Fy and Fz at 50% volume of material, (**A**) FBD of the ITAP in Sections, (**B**) The resulting topology of the individual Sections, (**C**) FBD of the ITAP constructed in one Section, (**D**) The resulting topology of Sections C and D.

within the ITAP's Sections C and D at different volume percentages of material, by removing material that had the lowest σ concentration within the geometry. Common

distribution of σ within the ITAP and what material can be removed to reduce the mass

#### *4.2. 3D MATLAB Analysis of the ITAP with Extreme Loads Applied*

For this study the assumptions and controls are the same as in Section 3.3 here the geometries used are the Standard and Vertical models. The simulations have been run between 50 and 95% of material in intervals of 5%, resulting in 30 separate geometries being developed for simulations. This study has identified the critical areas of material within the ITAP's Sections C and D at different volume percentages of material, by removing material that had the lowest σ concentration within the geometry. Common trends have been identified in these simulations. This has included: voids forming at the centre of the fixed support; divots manifesting within Section D; triangular interconnecting pathways forming within the voids in the centre of the models; and critical material kept along the boundary walls, most prominently the wall with the force directed towards it. All these observations are identical to the walking simulations 3D MATLAB study. As the percentage of the volume increases, the voids decrease in size as expected, but the pattern built into the structure is constant with the exception of the vertical model at 95%. The results from the latter simulation is sporadic because the material removed is not concordant with previous iterations.

#### **5. Discussion**

#### *5.1. Similarities of the ANSYS, 2D and 3D MATLAB Topology Studies on the ITAP for Walking Simulations*

Each simulation process has identified which regions of material are critical for the distribution of σ within the ITAP and what material can be removed to reduce the mass of the ITAP. As stated before, these results identify the material that can be removed as it is not critical to the structure of the ITAP. As such, the inverse of the results should be studied to correctly design the failure features of the ITAP, so that the device can be correctly used in day-to-day life, but will fail in a designed manner when exposed to extreme loads. Looking at Figure 13 which shows the visual results for each simulation method for sections C and D of the ITAP at 50% by volume of material, a clear pattern can be identified for further study. There are similarities in the reduction of material in each process, most notably the reduction of a large proportion of material at the fixed support of Section C, which is present in all models. Other significant similarities are the reduction of material identified in Section D, where material is removed in the opposite direction of the Frontal and Lateral forces for the divot. This is most prominently observed in both the ANSYS and 3D MATLAB simulations as the magnitude of the forces has been directly controlled and adjusted accordingly. The inverse of the results shows that the σ found within these locations is not substantial and therefore not critical for material reduction. The main areas where reduction could be focused on to achieve a controlled failure are the two quarters of the models that do not have any change in design; this can most prominently be seen in D in Figure 13 as well as A and B. Both the ANSYS and the Standard 3D MATLAB simulations are very similar in the locations of material reduction; the main differences arise due the starting geometry used.

The 2D MATLAB simulated model C is the most diverse, because the algorithm calculating the theoretical structure that could be capable of withstanding shear forces only at a point location rather than a dispersed load that is controlled. As such, the 2D MATLAB simulations do not work out the best structure to an individual load parameter, but rather work out the optimum geometry for shear loads being applied.

The results from this study have identified which internal structures should be utilised to distribute σ while saving weight. There are similarities in the simulated models such as the aforementioned void forming at the top of Section C in the ITAP, as well as the voids in the main body of the section being triangular based; this is similar to the results of E and F in Figure 13. However, from analysing and inverting the results provided from the use of Otomori's code, the areas of critical importance for distributing σ when designing the internal topology for shear σ are the boundary walls of the model. This is concordant with the range of percentages of volume simulated in this study. The resulting design

structures (Figure 9) show that as the maximum volume of the model increases, the size of the voids forming in the ITAP decrease. This is most prominent in the voids forming at the boundary wall and the voids forming where sections C and D meet. This study has identified where the areas of the topology of the ITAP should be changed depending on the volume of the material used in additive manufactured production processes for complex internal structures. The topology simulation using MATLABs' applied load only consisted of one shear force, rather than the three forces that are applied to the ITAP in motion. Furthermore, the force cannot be set to the applied loads identified in Section 2.1. As such, the results from this study can only be used as a guide. Bye analysing the studies it can be seen that only critical material is left, thus it can be inferred that this is the area where design changes are required for a controlled failure. This is concordant with the work of Bird et al. [13], where the location and depth of the safety notches affected the ITAP performances.

The findings from the 3D Vertical and Horizontal MATLAB studies support this claim. However, the findings identify that one section of the material has more importance than the other. In Figure 13 both E and F's cross-section of the ITAP has a thicker boundary wall on one side, located in the direction of F<sup>y</sup> or Fz, respectively. Both these locations are more regular and thicker than the opposite side of the ITAP. Inferring that the thicker side is more important for structural integrity and therefore should be the focus location when designing in the controlled location of failure. In total 60 different simulations were conducted to identify and consolidate the critical material needed for the structural stability of the ITAP when simulated with walking loads.

one solid model using ANSYS's control, ing ANSYS's simulation controls, **Figure 13.** Visual comparison between topology simulation of Section C and D of the ITAP at 50% volume of material, (**A**) one solid model using ANSYS's control, (**B**) 2 sections joined together using ANSYS's simulation controls, (**C**) 2D MATLAB simulations with only the shear force being applied to identify the potential areas of material reduction, (**D**) 3D MATLAB standard model, (**E**) 3D MATLAB vertical model, (**F**) 3D MATLAB horizontal model.

From the results of these studies (Table 2), the following conclusions and recommendations can be made for further development in simulation and physical testing of the topology of the ITAP:

 either side to increase contact with bone and reduce residual σ • Integrate a stress raiser on the ITAP at the end of Section B at the bone anchor, with appropriate fillets on either side to increase contact with bone and reduce residual σ under load at the connection between the ITAP and anchor.

**σ**

σ

σ

study's triangular

Simulated σ

σ AP

σ A



**Table 2.** Simulated σMax and mass reduction of the ITAP when studied at 50% topology reduction.

*5.2. Discussion and Similarities of the ANSYS and 3D MATLAB Topology Studies on the ITAP when Exposed to Extreme Loads*

From analysing both the data from ANSYS and 3D MATLAB topology studies for an ITAP exposed to extreme loads, see Figure 14, that are expected to cause fracture to the femur, the following observations and conclusions can be drawn:


one solid model using ANSYS's control, using ANSYS's simulation controls, **Figure 14.** Visual comparison between topology simulation of Sections C and D of the ITAP at 50% volume of material, (**A**) one solid model using ANSYS's control, (**B**) 2 Sections joined together using ANSYS's simulation controls, (**C**) 3D MATLAB standard model, (**D**) 3D MATLAB vertical model.

ent risk of using an ITAP is that because it is attached to the bone any σ on it receives is The findings from these studies are identical with the walking simulations, indicating that critical material should be removed from the model on the boundary wall, to which the shear the force is directed. However, the difference in the positioning of the forces is subject to change for the extreme loads because impact could occur from any direction. Thus, a failure feature must be designed that is capable of failing from any direction of shear force. This was concordant with all 34 simulations. The similarities of the location and σMax recorded for both the walking and extrema loads can be seen in Table 2. This also identifies the mass reduction calculated from the ANSYS calculations.

#### **6. Conclusions and Recommendations for Future Work**

bespoke to the users' lifestyle. The ITAP is an alternative prosthetic product that has the potential to remove the inflammation and pain from traditional prosthetics while replicating the feeling of having a leg via a direct skeletal implant, thus improving the quality of life of the user. The inherent risk of using an ITAP is that because it is attached to the bone any σ on it receives is transferred to the bone. To prevent further limb damage a failsafe safety notch design is considered. The optimisation challenge is to ensure the ITAP is functional for normal use while allowing for a designed failure if excessive forces are experienced. The following conclusions are made.


the 3D code is better suited to live modelling as there is more control on both the material properties and applied loads. The benefit of this is a visual guide to the critical material that is accurate to the model. However, the problem is the simulation does not yield numerical results and have limited geometry modelling. As such it is best suited for quick accurate visual modelling of cross-sections. In total 10 2D and 60 3D Matlab simulations where carried out within this research, all results yield from these simulations gave clear insight into the topology of the ITAP under load.


It is recommended that further research is undertaken on the design of the ITAP's failure feature utilising the topology study in this paper so that a safety design can be successfully developed for each individual case. Physical testing is also highly recommended as to ascertain the strength compression on different volume fraction of the ITAP, when tested with walking loads and destructive forces. Further simulation on different loads including running and jumping to fully facilitate the user's mobility would gain an insight in the loading capabilities of the ITAP.

**Author Contributions:** E.L.; formal analysis, investigation, software, visualization, methodology. C.A.G.; supervision, conceptualization, methodology, writing—original draft preparation. A.R.; resources, data curation, writing—review and editing. J.B. validation, formal analysis, visualization. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors would like to acknowledge the support of the Future Manufacturing Research Institute, College of Engineering, Swansea University and Advanced Sustainable Manufacturing Technologies (ASTUTE 2022) project, which is partly funded from the EU's European Regional Development Fund through the Welsh European Funding Office, in enabling the research upon which this paper is based. Further information on ASTUTE can be found at www.astutewales.com (accessed on 11 March 2021).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Micromachines* Editorial Office E-mail: micromachines@mdpi.com www.mdpi.com/journal/micromachines

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18

www.mdpi.com ISBN 978-3-0365-3270-7