*Article* **Force Prediction for Incremental Forming of Polymer Sheets**

#### **Gustavo Medina-Sanchez, Alberto Garcia-Collado, Diego Carou \* and Rubén Dorado-Vicente**

Department of Mechanical and Mining Engineering, University of Jaén, EPS de Jaén, Campus LasLagunillas, 23071 Jaén, Spain; gmedina@ujaen.es (G.M.-S.); acollado@ujaen.es (A.G.-C.); rdorado@ujaen.es (R.D.-V.)

**\*** Correspondence: dcarou@ujaen.es; Tel.: +34-953-211-712; Fax: +34-953-212-420

Received: 29 June 2018; Accepted: 30 August 2018; Published: 3 September 2018

**Abstract:** Incremental sheet forming (ISF) is gaining attention as a low cost prototyping and small batch production solution to obtain 3D components. In ISF, the forming force is key to define an adequate setup, avoiding damage and reducing wear, as well as to determine the energy consumption and the final shape of the part. Although there are several analytical, experimental and numerical approaches to estimate the axial forming force for metal sheets, further efforts must be done to extend the study to polymers. This work presents two procedures for predicting axial force in Single Point Incremental Forming (SPIF) of polymer sheets. Particularly, a numerical model based on the Finite Element Model (FEM), which considers a hyperelastic-plastic constitutive equation, and a simple semi-analytical model that extends the known specific energy concept used in machining. A set of experimental tests was used to validate the numerical model, and to determine the specific energy for two polymer sheets of polycarbonate (PC) and polyvinyl chloride (PVC). The approaches provide results in good agreement with additional real examples. Moreover, the numerical model is useful for accurately predicting temperature and thickness.

**Keywords:** incremental forming; FEM; force prediction; numerical model; semi-analytical model; specific energy

#### **1. Introduction**

Sheet forming by means of local deformations or incremental sheet forming (ISF) is a prototyping and small batch production solution to obtain 3D components developed in the late 20th [1]. During ISF, a tool follows a 3D path deforming a flat sheet clamped to a rigid frame. This technique is noted for their low cost and good forming capabilities, beyond the forming limit curve, compared to conventional forming processes [2]. The process has applications in automotive industry [3], aeronautical industry [4] and biomedicine [5].

The main ISF solutions are punctual incremental forming and doubled sided incremental forming assisted by a partial or full die (Figure 1). Single Point Incremental Forming (SPIF) is flexible, easy to implement low cost solution, and has been in the research spotlight in the last years. However, this process induces complex deformation mechanisms (shear, stretch and bending), which cause poor dimensional results [6], and limit its use in industrial applications [7].

**Figure 1.** Main incremental sheet forming (ISF) techniques. (**a**) Single point incremental forming; (**b**) Two point incremental forming; (**c**) Double-sided incremental forming with partial die; (**d**) Double-sided incremental forming with full die.

SPIF is usually accomplished by a conventional Computer Numerical Control (CNC) machine equipped with a clamping frame and a punch that follows a G-code program of instructions. The deformation mechanism is related to the forming forces and, therefore, a good knowledge about forming forces allows tackling several key SPIF problems: How to improve dimensional accuracy; how to avoid machine damage and wear; how to extend the technique to different materials, etc.

Despite the industrial interest of low cost small batch production of polymer parts at room temperature [8], the extension of SPIF to these materials is currently a challenge. In addition to the SPIF drawbacks mentioned above, it is difficult to model the deformation behavior of polymer materials. Most research efforts are focused on developing prediction models for metals of known mechanical properties. When coming to glassy polymers, they suffer relaxation after yielding and, therefore, they have a complex constitutive equation with strain rate and temperature material dependence.

The study of the forming force can be an adequate strategy to advance in forming of polymer sheets. Thus, progress can be made from the extension of forming force approaches for metal sheets to polymer sheets or by developing new forming force models. To the best knowledge of the authors, existing analytical and experimental forming force models were specifically developed for metals, and the scarce numerical models for polymers are not comparable to those for metals.

Regarding the force approaches for metals, analytical models define relations between process parameters, material mechanical properties and force. For example, Bansal et al. [9] relates the forming forces with thickness, meridional and circumferential stresses, contact surface and thickness. Besides experimental models, such as the regression equation as defined by Aerens et al. [10], estimate the forming force by means of experimental data. When considering numerical simulations, they are accurate for simple geometries but computationally expensive because of the nonlinearities produced by contact area changes. Behera et al. [11] shows an extensive classification of numerical simulations works focused on metal sheets.

The existing studies about incremental forming of different polymeric materials are limited and force prediction is not their main goal. In this regard, experimental works check new applications and asses how process parameters influence response variables such as force, roughness and incremental depth [12]. On the other hand, analytical and numerical approaches try to determine deformation mechanics and failure modes [13]. Concerning these solutions, it is worth mention the theoretical model based on membrane analysis developed by Silva et al. [14], the constitutive equation, based on overstress proposed by Alkas et al. [15,16], to model viscoplastic materials with only seven material parameters, and the numerical study of Nguyen et al. [17], who considers a viscoelastic material to estimate sheet thickness and spring-back.

The present work proposes a semi-analytical model based on fitting forming force measurements when forming truncated cones for different values of the deformed volume. This deformed volume is determined as a function of the contact surface and thickness. The novelty of this model is to extend the concept of specific energy, used in the orthogonal cutting model to SPIF. In this regard, a new concept called specific forming energy is introduced. This specific forming energy relates the forming forces with the geometric parameters characteristic of the process. The solution is simpler than analytical approaches.

Moreover, a coupled thermo-mechanical numerical model, with a hyperelastic-plastic constitutive equation suitable for polymer sheets, is presented. Unlike other aforementioned models with viscoplastic and viscoelastic constitutive equations, the proposed approach studies a hyperelastic-plastic material, only defined by six material parameters, which considers mechanical and thermal response. The material model parameters are determined minimizing the differences between experimental and simulated material stress-strain curves at different temperatures.

The structure of the paper comprises: Section 2 describes the considered assumptions, and the developed numerical and semi-analytical models; Section 3 explains the setup and tests carried out; Section 4 presents the main results, and discusses the goodness of the proposed approaches with respect to real measurements; and, finally, Section 5 presents the main conclusions.

#### **2. Forming Force Models**

#### *2.1. Semi-Analytical Force Model*

A simple method that relates the experimental axial forming forces, with a locally deformed volume, is proposed. Note that the main forming force acts along the *z* direction, so that is the one analyzed.

Following a similar reasoning such as that used in orthogonal cutting, the local plastic deformation induced during SPIF along the *z* direction requires a power *P*, so that it should be held that:

$$P = F\_{\overline{z}} \cdot \Delta z,\tag{1}$$

where *Fz* is the axial forming force and Δ*z* is the step down.

It stands to reason that the power needed to deform a sheet volume *V* should be constant for a specific material and temperature. This constant *U*, similarly to the machining case in Reference [18], can be called specific forming energy:

$$
\Delta U = \frac{F\_z \cdot \Delta z}{V}, V = \mathbf{S} \cdot t,\tag{2}
$$

where *S* is the affected punch-sheet area, and *t* is the mean thickness beneath the punch.

When knowing *S* and *t* at a specific position and *U* for a specific material, it is possible to solve the Equation (2) for obtaining *F*z.

The affected area can be calculated from the expression proposed by Bansal et al. [9] as:

$$S = \frac{\pi D}{4} (D/2 + \Delta z) + \frac{\pi D^2}{8(\gamma - a)} (\sin \alpha - \sin \gamma). \tag{3}$$

This area depends of the incremental depth Δ*z*, the forming tool diameter *D* and the angle of the sheet *α*, since *γ* can be calculated as:

$$\gamma = \arccos\left(1 - \frac{2 \cdot \Delta z}{D}\right). \tag{4}$$

*Materials* **2018**, *11*, 1597

The mean thickness of the deformed volume can be calculated from the sheet area deformed in a step down (Figure 2) and the arc of the punch perimeter in contact with the sheet as:

$$t = \frac{A}{(a' + \gamma') \cdot D/2} \tag{5}$$

$$\alpha' = \arccos(\cos(\alpha) - 2 \cdot \Delta z / D),\tag{6}$$

$$
\gamma' = \arccos(1 - 4 \cdot \Delta z / D). \tag{7}
$$

**Figure 2.** Deformed area defined.

Based on all of the above, the deformed volume calculated in Equation (2) depends exclusively on the geometric parameters: Δ*z*, *D*, and *α*. In summary, once the material characteristic *U* is experimentally determined, the axial forming forces can be estimated from the geometric parameters using Equation (2).

#### *2.2. Numerical Estimation of the Forming Force for Polymer Sheets*

The numerical model was carried out in a fully coupled-stress dynamic analysis using the ABAQUS commercial software. The analysis considers the inertia effect, the temperature-dependent of the material response, and the transient thermal response. It also includes the integration of the momentum and the heat flow equations coupling the material energy dissipation during plastic flow, rising the local temperature [19].

The Finite Element (FE) model used for the numerical simulation of the process can be seen in Figure 3. The forming tool and the backing plate were simulated as analytical rigid surfaces in order to perform an efficient numerical contact analysis. The polymer sheet is discretized employing 2052 S4RT elements with reduced integration to avoid the hourglassing, with an active degree of freedom to capture the variation of the temperature at nodes throughout the thickness by bilinear interpolation. The number of integration points used to capture the temperature variation throughout the thickness was 5. The mesh was generated taking into account the tool path that defines the final shape of the part. It was generated uniformly through angle and radius coordinates in order to maintain elements aligned and to generate minimal distortion during the polymer sheet deformation and thickness reduction. The minimum element length in the model was 3 mm, obtaining a minimum stable time increment of 1.2 ×10−<sup>6</sup> s that is small enough to avoid instabilities during the entire simulation. Note that,

mesh size influences the approach accuracy, nevertheless the authors selected the element size in order to limit the computational time to 10 h (for an Intel Core i8 CPU) maintaining an adequate concordance with the experimental results. The contact between the plastic shell and the forming tool was simulated by surface-to-surface interaction with penalty tangential behavior and hard contact behavior for normal direction.

**Figure 3.** FE model for the SPIF analysis.

The thermal contact conductance Ψ between the metal forming tool and the plastic shell was simulated by gap conductance independent of the contact pressure for polycarbonate (PC) and polyvinyl chloride (PVC) with value of 0.183 W/mK [20]. Other heat transfer mechanisms such as convection or radiation were not taken into consideration due to the low temperature of the plastic shell reached during the entire process. The role of the friction in the formability of the thermoplastics and metals during the SPIF process has been widely discussed by several authors [21–25], concluding that the effect of the friction in the thermoplastic has an important role due to the increase of the temperature of the material. In this work, the values of the dynamic coefficient of friction for the lubricated PC-Shell interface *μ<sup>k</sup>* for PC and PVC are taken from Ludena and Bayer [26]. Thermal properties like conductivity *kt* and specific heat *Cp* are taken from the manufacturer's data sheet. Table 1 summarized all thermal properties employed in the numerical model.


**Table 1.** Thermal properties for polycarbonate (PC) and polyvinyl chloride (PVC) employed in the numerical model.

The material model employed for the two thermoplastics is a non-linear hyperelastic model combined with J2-plasticity theory based on isotropic hardening: A simple hardening law that obtains good results with glassy polymers and low material parameters [27]. The hyperelastic component is based on the Arruda-Boyce eight chain model [28] that takes a non-linear Langevin chain statistics into account when deriving the strain energy density function. The predicted stress response of the eight-chain model can be written as follows:

$$
\sigma = \frac{\mu}{f \overline{\lambda^\*}} \frac{\mathcal{L}^{-1} \left( \overline{\lambda^\*} / \lambda\_L \right)}{\mathcal{L}^{-1} (1 / \lambda\_L)} \text{dev} [b^\*] + \kappa (f - 1) \mathbf{I}, \tag{8}
$$

where *μ* is the shear modulus, *k* the bulk modulus, and *λ* is the limiting chain stretch. The variable *b*<sup>∗</sup> = *J*−2/3*b* is the distortional left Cauchy-Green tensor, and *λ*<sup>∗</sup> is the applied chain stretch which can be calculated from:

$$
\overline{\lambda^\*} = \sqrt{\frac{\text{tr}[b^\*]}{3}}.\tag{9}
$$

The J2-plasticity component is based on isotropic hardening, that describes the size change of the yield surface *σ*<sup>0</sup> as a function of the equivalent plastic strain *ε*´*pl*. The model incorporates this effect by an exponential law defined by Equation (10):

$$
\sigma^0 = \sigma|\_0 + Q\_\infty \left(1 - e^{-bt^{pl}}\right),
\tag{10}
$$

with, *σ*|<sup>0</sup> the yield stress at zero plastic strain, *Q*<sup>∞</sup> the maximum change of the size in the yield surface and *b* defines the rate at which the size of the yield surface changes as plastic strain develops. If the material is rate independent, the yield condition is:

$$
\sigma^0 = q\_\prime \tag{11}
$$

with

$$q = \sqrt{\frac{3}{2}S : S} \, \Big|\, \tag{12}$$

and *S* is the deviatoric stress. The yield function is only dependent on the temperature and the equivalent plastic strain (*ε*´*pl*).

The procedure to determinate the values of the needed model parameters (*μ*, *k* and *λ*) is described in Section 3.1.2.

#### **3. Materials and Methods**
