*2.2. Experimental Setup*

The ISF tests were carried out using a Kondia®HS1000 3-axis milling machine (Kondia, Elgoibar, Spain). The details of the clamping system (Figure 1b) and the setup are described in detail in previous works [10]. The forming forces were measured by a table-type dynamometer Kistler® 9257B (Kistler ibérica SL, Granollers, Spain); surface roughness was determined by means of a Mitutoyo Surftest SV-2000 profilometer (Sariki, Cerdanyola del Vallès, Spain) (Figure 1c) and the maximum depth was recorded by a direct reading off the Kondia milling machine (Kondia, Elgoibar, Spain).

### *2.3. Design of Experiments*

Box-Behnken designs (BBD) [17] are three-level designs that allow second order response surfaces to be fitted efficiently. The design for the four factors consists of 27 experimental runs (Table 2) with 24 unique experimental settings (one replicate) plus three replicates at the central point. The explanatory variables came from the classical approach for metal ISF studies: thickness was kept constant for comparison and spindle speed was added because of its importance in polymer forming. The levels of the chosen variable were as follows:


\* Rotation is considered to be free when any rotation of the tool is due solely to the friction between the sheet and the tool itself.


**Table 2.** Design of experiments and results.

Note: \* In parentheses: 1 = maximum depth accomplished, 0 = sheet fracture.

#### *2.4. Analysis Procedure*

The methodology used to estimate the response surface models is the one described by Myers [18]. The main steps are summarized in Figure 2 and are as follows:


The approach we followed was to retain, in the model, the smallest subset of explanatory variables providing a significant regression test and a non-significant lack-of-fit test as well as good fit statistics (high R2, adj-R2, pred-R2 and low RMSE) together with an appropriate model adequacy. The reason for using a subset of explanatory variables, rather than all of them, is that the estimates of the coefficients will have smaller variance and the predictions will be more precise. The Shapiro Wilk normality test (SWNT) was used to check the normality of the residuals.

Note that coded variables (−1, 0 and 1, indicating low, medium and high level, respectively) were used to compare the size of the coefficients and that the significance level in all cases is α = 0.05. The statistical analysis was conducted using R software [19].

**Figure 2.** Model selection scheme.

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

The effects of varying various forming parameters on the three response variables under study (forming force, roughness and maximum depth) during an ISF process have been explored previously in a number of studies ranging from normal ISF to hot ISF, in experiments or simulations, and with regard to metals as well as polymers. We shall now go on to compare the results of these studies with the results we obtained for the two biocompatible materials in this work.

The factor levels (i.e., the process parameters) as well as response values for all 27 experimental runs are shown in Table 2. Note that the response value ΔRa has been introduced; we shall discuss the reasons for this in the following subsections, in which the statistical models proposed for each material are shown and briefly discussed.

#### *3.1. Maximum Axial Force (Fz Max)*

#### 3.1.1. PCL

The model that best fits Fz max for PCL is as follows:

$$
\circ = 304.83 + 92.12 \cdot \text{Dt} - 31.01 \cdot \text{S} - 22.93 \cdot \text{Dt} \cdot \text{S}, \tag{1}
$$

It is a very good model which describes 94% of the variability of Fz max and has a predicted R<sup>2</sup> of 92%. The residuals of the model are homoscedastic, independent and identically normally distributed (SWNT *p*-value = 0.06).

Two of the first order factors, tool diameter (Dt) and spindle speed (S), are significant, as is the interaction between the two (Table 3) with Dt being the most influential because its coefficient is higher (three times higher than the coefficient of S). In general, the higher Dt (14 mm) produces higher Fz max values. It is clear that the portion of sheet to be formed using a greater tool diameter requires higher forces, as has been pointed out previously for metals (Al7075-O in Li et al. [6], Al3003-O in Bahoul et al. [7] and AISI304 in Centeno et al. [8]) as well as for polymers (PVC in Bagudanch et al. [9]). In contrast, an increase in spindle speed reduces the forces required. The mechanism able to explain

this fact involves the change in the friction conditions, as the heat originating from the increase in friction decreases the various forming force values and, as a result, the mechanical behavior of the polymer material may change. This is also consistent with the work presented by Centeno et al. [8,20] and Baharudin et al. [21] for metal materials and, in the case of polymers, by Davarpanah et al. [11] for PLA and PVC, by Lozano-Sánchez et al. [16] for PCL and UHMWPE and by Bagudanch et al. [22] for PVC, PC, PP, UHMWPE and PCL.

The interpretation of the interaction effect (Dt·S) is as follows: a reduction in S (from 2000 rpm to 0 rpm) slightly increases Fz max when Dt = 6 mm, while the increment is more acute when Dt = 14 mm. As we have said, the increase in spindle speed increases the friction and temperature, causing lower force values, which agrees with results independently observed. Hence, it would be expected that, at higher tool diameters, the heat would concentrate even more in the forming region resulting in a temperature increase [22] which, in turn, would decrease the value of the required forming force; however, this does not occur according to the results we obtained. In the Dt·S interaction, it seems that the effect of the increased contact zone at higher tool diameters—which causes an increase in the forming force [4]—is a more dominant influence than the increase in temperature caused by the higher spindle speed.


**Table 3.** Maximum force (Fz max) model results.

#### **Table 3.** *Cont.*

#### 3.1.2. UHMWPE

The model that best explains Fz max on UHMWPE material is as follows:

$$\begin{aligned} \text{dy} &= 626.33 + 187.88 \cdot \text{Dt} - 98.82 \cdot \text{S} - 15.64 \cdot \text{F} + 6.40 \cdot \text{Az} - 64.90 \cdot \text{Dt} \cdot \text{S} \\ &+ 46.34 \cdot \text{Dt} \cdot \text{Az} + 38.24 \cdot \text{S}^2 \end{aligned} \tag{2}$$

This is a very good model able to explain 98% of the variability of Fz max and with a high capability of predicting new response values (pred R2 = 96%). The residuals of the retained model are i.i.d. distributed (SWNT *p*-value = 0.3303).

Like the PCL model, the UHMWPE model depends highly on Dt which is the most influential first order factor and also appears in the two significant two-way interactions. In general, high values of Dt result in higher values of Fz max. It is already known that an increase in tool diameter entails an increase in force. The second most influential factor is spindle speed, S, and it has a significant quadratic effect and interacts with Dt. The feed rate, F, is the least influential factor because of its lower coefficient; however, it is still significant. The effect of F can be appreciated in the (F, Dt) surface plot in Table 3; when F moves from 1500 to 3000 mm/min, Fz max is reduced on average by 15.64 N.

The interaction effects can be better understood from the surface plots. With regard to Dt·S, it can be appreciated how changing from S = 2000 rpm to S = 0 rpm increases Fmax in all cases, however this increment is more acute when Dt = 14 mm than it is when Dt = 6 mm.

The other significant interaction (Dt·Δz) has the opposite effect: on the one hand, when Dt = 6 mm, a reduction in Δz implies an increase in Fz max; on the other hand, when Dt = 14 mm, a reduction in Δz implies a reduction in Fz max. The reduction in step size decreases the portion of sheet being formed, which should also lead to a reduction in the force: the model shows this trend for high tool diameters. However, it is not the case for low tool diameters for which the forming force slightly increases.

#### *3.2. Surface Roughness (Ra)*

Roughness assessment needed to be conducted in a different way for each material, since the blank sheet was produced in different ways in each case. Furthermore, since it is known that sheet metal mark orientation can affect the results of roughness measures [23], it should be pointed out that the roughness values were obtained, in both cases, perpendicular to the step-down direction.

Tool diameter and step down have been identified as parameters that influence surface roughness, Ra, [13,14] mainly in metals, whereas for PVC and PC, step down and spindle speed with some interactions were significant [23]. The manner in which the tool diameter and step down parameters improve or worsen the surface roughness can be seen in the diagrams in Figure 3, which shows how increased Dt reduces Ra, while increased Δz increases it.

**Figure 3.** (**a**) Effect of different tool diameters, Dt, on surface roughness, Ra (**b**) Effect of different step down, Δz, on surface roughness, Ra.

#### 3.2.1. PCL

The model that best explains Ra for the PCL material includes all four factors of the experiment:

$$\begin{aligned} \mathbf{y} &= 0.5\mathbf{\hat{r}} + 0.43 \cdot \mathbf{D} \mathbf{t} + 0.49 \cdot \mathbf{S} - 0.1\mathbf{\hat{r}} \cdot \mathbf{F} + 0.1 \cdot \mathbf{\hat{A}} \mathbf{z} + 0.42 \cdot \mathbf{D} \mathbf{t} \cdot \mathbf{S} - 0.22 \cdot \mathbf{F} \cdot \mathbf{A} \mathbf{z} \\ &+ 0.51 \cdot \mathbf{S}^2 + 0.16 \cdot \mathbf{F}^2 + 0.21 \cdot \mathbf{A} \mathbf{z}^2 \end{aligned} \tag{3}$$

The model explains 88% of the variability of the response Ra and has a good predicted R2 of 74%. Residuals are independent, homoscedastic and normally distributed (SWNT *p*-value = 0.3082). The model is complex because of the high numbers of terms included in it, but the surface plots shown in Table 4 can help in interpreting it.



6 FRGHG

 

 

) FRGHG

In general, and in contrast to previous studies, a higher Dt is associated with a higher Ra (worse surface quality). However, this effect is more acute when S is at its high level (S = 2000 rpm), due to the friction and the heat that the surface is receiving [9]. Similarly, a higher Δz level is associated with a higher Ra, which is as expected, but again, this effect is more acute when the feed rate is low, F = 1500 mm/min, since at lower feed rates, contact time with the surface is higher, therefore there is more heat, which worsens surface quality. The effect of the quadratic terms S and F show a simple minimum pattern on the response surface. The stationary point of the response surface is Dt = 1.3, S = −1, F = 0.6, Δz = 0.1 in coded units or Dt = 15 mm, S = 0 rpm, F = 2688 mm/min, Δz = 0.36 mm. This stationary point represents a minimum for the Ra value.

#### 3.2.2. UHMWPE

The parameter of roughness for UHMWPE analyzed here will be the difference between the roughness before and after processing (ΔRa = Rabefore processing–Raafter processing). In contrast to PCL in which the sheets proceed direct from a mold, in the case of UHMWPE, the sheets were sliced in layers from a thicker molded sheet, which meant the surface already began with a texture whose starting roughness was high, at around 1.05 μm on average. With one exception, the ΔRa value is always positive, which means that the processing operation improves the roughness of the material. The roughness is reduced on average by 0.464 μm with 95% CI = [0.36, 0.57].

The model selected to explain the differences in roughness in UHMWPE is as follows:

$$\mathbf{y} = 0.4\mathbf{\hat{o}} + 0.08\mathbf{\hat{Dt}} - 0.03\mathbf{\cdot}\mathbf{S} - 0.03\cdot\mathbf{F} - 0.36\cdot\mathbf{Dt}\cdot\mathbf{S} \tag{4}$$

The model, however, explains only 21% of the variability of the response variable and has a bad predicted R2. It was decided to keep F in the model to achieve a non-significant lack-of-fit test. The general regression test is significant. Residuals are normally distributed (SWNT *p*-value = 0.87).

The most important factor explaining the differences in roughness in UHMWPE is the interaction between Dt and S which has a negative coefficient. The surface plots shown in Table 4 can help interpret the model. They show a saddle point (the stationary point) near the center of the plot. From the center point, increasing or decreasing Dt and S at the same time produces a reduction in ΔRa, while increasing one factor and decreasing the other leads to an increase in ΔRa. In other words, Ra after processing is slightly reduced when Dt and S are both at their most positive or most negative values (low ΔRa values), whereas Ra is highly reduced when Dt and S are at opposite levels (very low ΔRa values).

#### *3.3. Maximum Achieved Depth (Zmax)*

The response variable, maximum achieved depth, Zmax, both on UHMWPE and PCL, cannot be analyzed using a response surface model because the data shows a highly right-skewed distribution which is truncated at different depths (42.7 mm for a step down of 0.35 mm or 43.0 mm in the other cases) depending on the experimental settings. For example, for UHMWPE, 20 out of 27 experiments reached the maximum specified depth (Zmax = 100%), five were in the interval [90,100) and two in the interval [80,90). All tested models have a significant lack of fit test and show a dependence on the residuals vs. fitted values, mainly when Zmax = 100%.

The reason for analyzing the Zmax response is to find the point, z, at which the material breaks and to find which factors explain the breaking depth (Zmax). Generally, this can be modeled using survival (or reliability) analysis which aims to analyze the relationship of time to an event. In this context, the event is defined as the breaking (or failure) of the polymer and time is represented by the depth, z.

The experimental data shows many censored observations, that is, experiments in which the material did not break. Each experiment is an observation of the type (zi,ci) where zi is Zmax and ci = 1 if there has been a failure, and ci = 0 (censored) if otherwise. These ci values are shown in parentheses in Table 2 alongside the Zmax values.

The aim of this section is to explain the probability of the material surviving at depth z, that is, the survival, or reliability, function: S(z) = 1 − F(z) = P(Z ≥ z)–being careful not to confuse S(z) with the factor S = Spindle-speed. Another important function in survival analysis is the hazard function, λ(z). It assesses the instantaneous risk of demise at next mm, conditional on survival to that depth. In other words, it is the expected number of events that will occur in the next mm given that there has not been an event to that depth (z).

In this analysis, the Kaplan Meyer method (a non-parametric approach) is used to estimate S(z). Moreover, the Mantel-Haenszel test (log-rank test) is used to check for differences between survival curves (factor levels).

It is finished by fitting a Cox proportional hazard regression model (semi parametric approach) to determine whether the model including the significant factors is significant itself and to look for significant interactions. It is used with the likelihood ratio, Wald and logrank tests.

Note that the number of experiments at each level is not balanced. For example, there are six experiments at S = −1, 15 at S = 0 and 6 at S = 1. This is because data was not collected to carry out a survival analysis, but rather a BBD of experiments.

#### 3.3.1. PCL

Kaplan Meyer curves (Figure 4) can be interpreted as follows: at z = 20 mm, all observations/experiments are without event and the survival (S(z)) is 1, or, equivalently, 100%. The solid lines, colored according to factor levels, show the events with vertical drops (at each z). Colored curves showing different patterns indicate that there are differences in survival between factor levels, such as is the case of factor Dt (Mantel-Haenszel test *p*-value = 0.033): it is the only significant factor in PCL survival.

**Figure 4.** Significant survival curves according to factor levels on PCL and UHMWPE.

Note that the survival curves are not balanced: there are more experiments on Dt = 10 mm than on Dt = 6 or 14 mm (15 vs. 6). However, it is noticeable that survival is lower for Dt = 14 mm.

The Cox proportional hazard regression model has a unique significant factor: Dt. The coefficient of S (b = 0.95) is interpretable in its exponential form (e0.95 = 2.59) as the multiplicative effects of the hazard. That is, increasing Dt one level (e.g., from −1 to 0) increases the danger of breaking by, on average, a factor of 2.59. The overall tests of significance *p*-value (likelihood ratio, Wald and logrank) are significant, indicating that the model is appropriate.

Nevertheless, when S is included in the Cox model together with Dt, the individual regression coefficient of S has a *p*-value of 0.07 and the overall model is significant. In this case, the coefficient associated with S is b = −0.70 meaning that increasing S by one level reduces the danger of breaking by e−0.70 = 0.50. Note, however, that the interaction between Dt and S is not significant. More experiments should be carried out in order to confirm the significance of S in the survival function.

### 3.3.2. UHMWPE

Kaplan Meyer curves for Zmax on UHMWPE show significant differences only across S levels (Mantel-Haenszel test *p*-value < 0.001). From Figure 4, it can be seen that S(z) is lower when S = 2000 rpm. No breaks were observed for S = 1000.

The Cox proportional hazard regression model has S as a unique significant factor. Its coefficient is βˆ = 3.38 meaning that increasing S one level increases, on average, the danger of breaking by e3.38 = 29.49.

The Kaplan-Meyer curves are shown globally without being separated by factor levels in Figure 5. In the graphics comparing PCL versus UHMWPE, it can be observed that PCL sheets break earlier (z around 25) and in greater quantity than UHMWPE sheets in which the first breakage appears at around z = 35.

**Figure 5.** Kaplan-Meyer curves without separating by factor levels for (**a**) PCL and (**b**) UHMWPE.

#### **4. Conclusions**

In this work, SPIF experimental tests using two different polymers biocompatible have been carried out following a Box-Behnken design for four factors and a survival analysis. The maximum forming force, surface roughness and Zmax response achieved in the experiments has been statistically analyzed and empirical models for each material have been obtained. Using the proposed models, it is possible to control the Fz max, Ra and Zmax.

Among the process parameters analyzed and according to the data summarized in Table 5, spindle speed and tool diameter have been found to be the most influential parameters in terms of maximum forming force variation for both materials. From the roughness analysis, one can observe again the importance of spindle speed and tool diameter on these biocompatible polymers, as both involve an increase in temperature due to either the friction at higher spindle speeds or the increase in surface contact and contact time between tool and sheet, both of which worsen the surface because the materials flow, degrade and lose their properties.

Finally, the response Zmax, on UHMWPE and PCL, cannot be analyzed using a response surface model because data is highly right-skewed and truncated at different depths. However, taking into account the objective of the paper, the data can be analyzed using a novel method in this field: survival analysis. The results have shown that the most important factor is spindle speed for UHMWPE and tool diameter tool for PCL.


**Table 5.** Summary of the coefficients of the selected models (bold indicates significance (α = 0.05)).

**Author Contributions:** M.S., M.L.G.-R. and I.B. conceived and designed the experiments. M.S. and M.V.-M. carried out data and the statistical analysis. M.S. and I.B. performed the SPIF experiments. M.L.G.-R. and I.F. directed the research. M.S. wrote the manuscript. The manuscript was finalized through contributions from all authors, and all authors approve the final manuscript.

**Funding:** This research has received funding from the University of Girona (MPCUdG2016/036), the Spanish Ministry of Education (DPI2016-77156-R) and the Catalan Agency for Management of University and Research Grants (2017-SGR-0385).

**Acknowledgments:** The authors would like to thank the Centro de Investigación de Química Aplicada (CIQA) for supplying the PCL sheets. Finally, we are grateful to Marc López for his collaboration during the experimental work.

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

#### **References**


© 2018 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/).
