1 Introduction

Friction stir extrusion (FSE) is a solid-phase processing technique that is used to consolidate and extrude high-performance wires, rods, and tubes from metal powder, chips, flakes, or a solid billet in a single step. The process was invented at The Welding Institute (Cambridge, UK) and patented by Thomas et al. [1] in the early 1990s. Subsequent development at Pacific Northwest National Laboratory (PNNL) has led to patents for using FSE to extrude tubes [2]. FSE process is nominally solid state, and the maximum temperature is self-limited by the loss of flow stress as the extruded material approaches its melting point. Generally, no auxiliary heating is required and only the material near the interface between the rotating die and the charge is heated due to material plastic deformation and friction, which reduces the energy and time required for extrusion. The enhanced microstructure and material mixing inherent to FSE yield products with much better ductility and energy absorption capacity than those made by conventional extrusion [3]. These advantages have spurred a rapid development in FSE in recent years, and FSE has shown tremendous potential in industries such as automotive, aerospace, advanced manufacturing, and metal waste recycling.

FSE is a thermomechanically coupled process involving severe shear plastic deformation, material mixing, and heat generation from plastic deformation and friction. Experimentally measuring the thermomechanical behaviors (e.g., material flow, strains, strain rates, and temperatures) of the processed material is challenging, expensive, or impossible. Thus, numerical simulation is necessary to obtain a better understanding of the complex physics and material behaviors during FSE. Due to the associated computational challenges and cost, few validated models for FSE are available in the literature. Behnagh et al. [4] proposed a multistep 2-D coupled Eulerian–Lagrangian (CEL) model in Abaqus to simulate the processing temperatures in FSE of magnesium chips and consequent distributions of grain size and microhardness in the material with different process parameters. Baffari et al. [5] developed a 3-D, Lagrangian, thermomechanical model in Design Environment for FORMing (DEFORM)® software to simulate FSE of AZ31 milling chips. They showed that the distributions of material strain, strain rate, microstructure, and temperature can be predicted well by the model. Zhang et al. [6] proposed a 3-D CFD model in the ANSYS Fluent for FSE of AA6061 wires with consideration of both heat transfer and material flow. Recently, Karwa [7] developed a finite element method (FEM) model in FORGE NxT® for FSE, and comprehensive studies were carried out on the simulation results to analyze all aspects of the FSE process.

As a solid-phase processing similar to FSE, friction stir welding (FSW) [8] has been intensively investigated in the past few decades, and the relevant simulation approaches for FSW may be applicable to FSE. The existing FSW models can be categorized into fluid-based and solid-based approaches, and their detailed numerical comparison was investigated by Bussetta et al. [9, 10]. In the fluid-based approach, the workpiece is considered as a non-Newtonian fluid due to the large material deformation that occurs in FSW, and the problem is solved using the principles of computational fluid dynamics (CFD). After Ulysse [11] first developed a CFD model for FSW using the fluid dynamics analysis program FIDAP, many studies were carried out based on fluid-based approach. Seidel and Reynolds [12] proposed a 2-D thermal model based on non-Newtonian flow around a circular cylinder. By assuming full sticking contact condition between the tool pin and welded zone, they observed significant vertical material mixing occurred during FSW. Colegrove and Shercliff [13] developed a 3-D thermal and material flow model for FSW of AA7449 thick plate in ANSYS Fluent, and the developed model was used to investigate different pin profiles and rotational speeds on the pin traversing force. Bastier et al. [14] conducted a two-step simulation of the steady state phase of FSW, in which the thermomechanical behavior of FSW was first studied using a CFD model, and the residual stresses induced by the process were studied in the second step. Albakri et al. [15] proposed a 3-D, thermomechanical CFD model for FSW of AZ31. The flow stress constitutive relation was experimentally corrected which yields more consistent temperature, material flow, and strain rate predictions. Despite its success in predicting experimentally validated quantities during FSW, one issue about the fluid-based approach is that it cannot include material hardening, because it assumes the workpiece as a non-Newtonian fluid with an effective viscosity. Moreover, the material elastic strains are assumed to be negligible and are ignored. This is valid for FSW but may lead to inaccurate elastic stage, extrusion force, and pressure predictions in FSE. On the other hand, the solid-based approach is based on the principles of computational solid mechanics (CSM), which allows the inclusion of material strength, elasticity, hardening, and strain rate dependence in the formulation. The solid-based approach also allows the modeling of complex contact conditions that exist at the interface of the workpiece and tool. The existing solid-based models for FSW include Lagrangian model and Eulerian model. Jain et al. [16] proposed an implicit Lagrangian model to simulate the plunging, dwelling, and welding stages of a FSW process. Zhang and Zhang [17] established a Eulerian model for FSW which addressed the large material deformation issue without remeshing. To balance the drawbacks of the Lagrangian and Eulerian finite element simulations, arbitrary Lagrangian Eulerian (ALE) model and CEL model have been applied to FSW problems to avoid the need for remeshing to handle large deformations in the stirred zone and describe the history of field variables more accurately. For example, Assidi et al. [18] implemented a thermomechanical ALE model for FSW process in the Forge3® finite element software. The Norton’s friction model was used for the contact between the tool and work piece, which provides good predictions of welding forces. Dialami and coworkers proposed an apropos kinematic simulation framework (ALE/Eulerian/Lagrangian) for FSW [19,20,21] as well as many enhancements to the proposed framework [22], such as a two-stage speed up strategy [23], subgrid scale stabilized mixed \({\mathbf {v}}/p\) formulations to avoid numerical locking and oscillations [24], and material flow tracing method [25]. Dialami et al. also presented a local-global strategy to predict residual stresses during FSW process [26]. Gupta et al. [27] developed a CEL model in Abaqus to predict weld interface and temperature in the friction stir scribe joining of dissimilar materials. More recently, Robe et al. [28] developed a R-ALE FEM approach to simulate heat transfer during FSW of AA2xxx and AA7xxx joint.

Fig. 1
figure 1

a Cross-sectional schematic of wire FSE apparatus; b flat-face die; c AA1100-O solid billet with AA2050-O marker wires; d extruded wire from the billet

Full size image

All the above-mentioned models for FSW and FSE were belonged to grid-based methods, where the local approximation of functions and their derivatives for solving the governing partial differential equations (PDEs) is based on a predefined grid, mesh, or stencil. Unlike the grid-based methods, SPH is an Lagrangian particle method first invented by Lucy [29] and Gingold and Monaghan [30] for astrophysics simulations. The general idea of SPH is to approximate a continuous field using a set of kernel functions centered at discrete particles, which carry the physical properties of the system [31]. Because no grid or mesh is needed to construct the set of algebraic equations with nodal unknowns for the field variables, SPH is capable of handling extremely large material deformations and mixing, which would cause severe mesh distortion issues in grid-based methods. The particle representation in SPH also enables a natural description of material interfaces, free surfaces, and moving boundaries during the simulation process, which is a challenging task for grid-based numerical methods. Furthermore, the Lagrangian nature of SPH allows for explicit tracking of material field variables. Over the past three decades, SPH has been mainly developed and applied to CFD simulations and thermal problems [32]. Libersky and Petschek [33] first used SPH to solve a solid mechanics problem based on the elastic–perfectly-plastic material constitutive model. Since then, interest has grown in applying SPH to a variety of solid mechanics problems. For instance, Gray et al. [34] addressed the tensile instability issue in SPH by introducing an artificial stress to the principal stresses. The improved SPH model was applied to elastic dynamics problems such as oscillating beams, colliding rings, and brittle solids. Gray and Monaghan [35] introduced a SPH formulation for elasticity problem to predict the stress fields and brittle fracture caused by stress concentration. It was found that SPH is a suitable tool for modeling damage and fracture due to its meshless nature. Cleary et al. [36] applied SPH to a 2-D metal forging problem with AA6061, and it was shown that the very large plastic material deformation during forging process can be well handled by SPH. Ba and Gakwaya [37] proposed a total Lagrangian SPH formulation for thermomechanical problems under large deformation. The effectiveness of the developed model was tested on a high velocity Taylor bar impact problem and hot forging problem.

As an ideal method to simulate solid-phase processing problems, SPH was first used to predict the material flow and temperature in FSW by Tartakovsky et al. [38] using the fluid-based approach. Later, Bhojwani [39] implemented solid-based SPH in the LS-DYNA program for FSW. In this work, only stress and strain fields were represented, using the Johnson–Cook constitutive law for the stress flow. The heat generation and material thermal softening were not considered. Pan et al. [40] predicted the temperature, grain size, microhardness, and texture evolution during FSW using a 3-D fluid-based SPH model. Fraser et al. [41] developed a robust and efficient SPH code, SPHriction-3D, for finding the optimal FSW process parameters. A coupled, thermomechanical, solid-based SPH model was implemented to predict the temperature, stress, strain, damage, and defect formation during the process. Ansari and Behnagh [42] developed a thermomechanical SPH model in Abaqus to simulate the plunging phase of FSW. More recently, Stubblefield et al. [43] created a fully thermomechanical coupled SPH approach for additive friction stir-deposition (AFS-D) process. As to work using SPH to simulate FSE, no studies have apparently been published yet.

In this study, a 3-D, thermomechanically coupled, solid-based SPH model is proposed to simulate FSE of aluminum wire from a billet. Johnson–Cook model is utilized to describe the thermo-viscoplastic behavior of the billet material during FSE, and Coulomb–Tresca friction law is used to describe the contact condition between the die tool and billet. Heat is generated from both the plastic deformation of the bulk material and the frictional sliding at the interface between the tool and billet. The evolutions of the material flow, strain, temperature, and extrusion force were predicted to gain a better understanding of the extreme thermomechanical conditions during the FSE process. The proposed model is thoroughly validated against the experimental observations for different sets of processing parameters. The remainder of this paper is outlined as follows: Sect. 2 describes the experimental procedure of FSE. The thermomechanical SPH model, together with the governing equations and constitutive model, is presented in Sect. 3. The model setup is detailed in Sect. 4. The numerical results are presented, experimentally validated, and discussed in Sect. 5. Finally, conclusions and remarks are drawn in Sect. 6.

2 Experimental procedure

In this study, the FSE process was performed by a one-of-kind Shear Assisted Processing and Extrusion (ShAPE) machine designed by PNNL. Figure 1a shows the schematic representation of the FSE apparatus, which includes a cylindrical die tool with a flat face and a central orifice (Fig. 1b), a ring-shaped container for the billet, and a backing block at the billet end. The billet ring was placed in a larger water-cooled container. As the die tool rotates and advances, high shear and compressive forces are applied to the billet, plasticizing and generating heat from both friction and plastic deformation to soften the billet material, which is then extruded through the orifice of the die to form a wire, as shown in Fig. 1d.

A flat-face die fabricated from H13 tool steel with a 100:1 extrusion ratio (defined as the ratio of the cross-sectional area of the billet to that of the extruded wire) was used to extrude 2.54 mm diameter wire from a 25.4-mm-diameter cylindrical billet of AA1100-O. A marker insert technique was used to visualize the material flow in the wire. Specifically, 2.54 mm diameter AA2050-O wires were inserted into the center and at 1/3 of the radius (1/3r) from the center of the billet, as shown in Fig. 1c. The marker wires extended completely through the billet longitudinally. Constant die rotational and advance speeds were applied in this process, resulting in a constant extrusion speed with varied extrusion force. The die rotational and advance speeds used in the experiment are shown in Table 1, along with the final extruded wire lengths for each case. The extrusion force, die position, torque, and power were recorded by the ShAPE machine. A K-type thermocouple was placed 0.5 mm below the die face and 5 mm from the die central axis through a predrilled hole to measure the processing temperature of die face. After extrusion, the produced wires were sectioned transversely at different locations. The cross sections of the wire were ground, polished, etched, and imaged by optical microscopy to show the distribution of the marker wire at different locations along the wire. The strains at 1/3r of the wire cross section at the sectioned locations were estimated based on the spiral morphology of the 1/3r marker wires at those locations. The methodologies used for marker wire insertion and strain calculation can be found in greater detail in Li et al. [44].

Table 1 FSE parameters and lengths of extruded wires
Full size table

3 Model description

3.1 Smoothed particle hydrodynamics (SPH) method

In SPH, a continuous function \(f(\varvec{x})\) and its spatial derivative \(\nabla \cdot f(\varvec{x})\) at point \(\varvec{x}\) in a continuous domain \({\Omega }\) can be approximated as \(\langle f(\varvec{x}) \rangle \) and \(\langle \nabla \cdot f(\varvec{x}) \rangle \), respectively, by a convolution integral based on the information of its neighboring particles as follows:

$$\begin{aligned}&\langle f(\varvec{x}) \rangle =\int _{{\Omega }}f(\varvec{x}')W(\varvec{x}-\varvec{x}',h)d\varvec{x}' \end{aligned}$$
(1)
$$\begin{aligned}&\langle \nabla \cdot f(\varvec{x}) \rangle =-\int _{{\Omega }}f(\varvec{x}')\cdot \nabla W(\varvec{x}-\varvec{x}',h)d\varvec{x}' \end{aligned}$$
(2)

where \(W(\varvec{x}-\varvec{x}',h)\) is the kernel function in terms of the spatial distance between the function evaluation point \(\varvec{x}\) and the interpolation point \(\varvec{x}'\), and the smoothing length h defines the support domain of the kernel function. A h value between 1.05 and 1.3 is typically suggested, \(h=\) 1.2 is used in this study. The kernel function W should satisfy unity, support compactness, positivity, the Dirac delta function property, and smoothness conditions [32]. Among various kernel functions, the cubic B-spline kernel [45] is used in this study, which can be expressed as:

$$\begin{aligned} {W}(R,h)\,=\alpha _{d}\, \left\{ \begin{array}{lll} \displaystyle \frac{2}{3}-R^2+\frac{1}{2}R^3 &{}&{} 0\le {R}<1\\ \displaystyle \frac{1}{6}(2-R)^3&{}&{} 1\le {R}<2\\ 0&{}&{} R\ge {2} \end{array}\right. \end{aligned}$$
(3)

where R is the relative distance between two particles \(\varvec{x}\) and \(\varvec{x}'\), \(R=r/h=|\varvec{x}-\varvec{x}'|/h\), with r denoting the distance between the two points. The coefficient \(\alpha _{d}={3}/(2\pi h^3 )\) is defined in 3-D space. The definition of \(W=0\) when \(R\ge 2\) in Eq. (3) is to satisfy the support compactness condition. Therefore, the support domain with nonzero kernel function can be defined as \(\kappa h\), where \(\kappa \) is a constant related to the smoothing function for a particle at \(\varvec{x}\). The gradient of the cubic B-spline function can be expressed as:

$$\begin{aligned} \nabla {W}(R,h)\,=\nabla R \alpha _{d}\, \left\{ \begin{array}{lll} \displaystyle \frac{3}{2}R^2-2R &{}&{} 0\le {R}<1\\ \displaystyle -\frac{1}{2}(2-R)^2&{}&{} 1\le {R}<2\\ 0&{}&{} R\ge {2}\\ \end{array}\right. \end{aligned}$$
(4)

In SPH, the continuous domain \({\Omega }\) of interest is discretized by a finite number of particles located at \(\varvec{x}_i(t), i\in {1,2,..,N}\). Each particle i carries unique physical properties such as density, velocity, stress, and temperature. These particles can be used not only to represent the field properties of the domain, but also for numerical integration, differencing, or interpolation. Now the continuous form of Eqs. (1) and (2) can be written in quadrature formulas as:

$$\begin{aligned}&\langle f(\varvec{x}_i) \rangle =\sum _{j=1}^{N_i}\frac{m_j}{\rho _j}f(\varvec{x}_j)W(\varvec{x}_i-\varvec{x}_j,h) \end{aligned}$$
(5)
$$\begin{aligned}&\langle \nabla \cdot f(\varvec{x}_i) \rangle =-\sum _{j=1}^{N_i}\frac{m_j}{\rho _j}f(\varvec{x}_j)\cdot \nabla W(\varvec{x}_i-\varvec{x}_j,h) \end{aligned}$$
(6)

where \(N_i\) is the number of particles within the support domain of W at particle i and \(m_j\) and \(\rho _j\) are the mass and density of particle j. \({m_j}/{\rho _j}\) defines the “weight” of particle j in the quadrature. A schematic concept of a 2-D discrete particle approximation in SPH is shown in Fig. 2.

Fig. 2
figure 2

Schematic view of a 2-D SPH particle approximation for particle i using neighboring particles within the support domain of the kernel function W

Full size image

3.2 Governing equations

In this study, the billet is considered as a thermo-elastoplastic solid. Therefore, the billet material deformation is governed by the equations of mass, momentum, and energy conservation:

$$\begin{aligned}&\frac{\mathrm{d}\rho }{\mathrm{d}t}=-\rho \nabla \cdot \varvec{v} \end{aligned}$$
(7)
$$\begin{aligned}&\frac{\mathrm{d}\varvec{v}}{\mathrm{d}t}=\frac{1}{\rho }\nabla \cdot \varvec{\sigma }+\varvec{b} \end{aligned}$$
(8)
$$\begin{aligned}&\rho C_p\frac{\mathrm{d}T}{\mathrm{d}t}=\varvec{\sigma }:(\nabla \cdot \varvec{v})-\nabla \cdot (k\nabla T)+\dot{q} \end{aligned}$$
(9)

where \(\rho \) is the solid density, t denotes the time, \(\varvec{\sigma }\) is the Cauchy stress tensor, \(\varvec{v}\) and \(\varvec{b}\) are the velocity and external force vectors, respectively, \(C_p\) denotes the specific heat capacity, T is the temperature, \(\dot{q}\) is a heat source and sink term that describes the heat generation and loss in the system, \(\nabla \) is the spatial gradient operator, and \(\nabla \cdot \) denotes the spatial divergence operator. With the SPH kernel and particle approximation shown in Sect. 3.1, the three conservation equations can be written in discrete form and then solved for the unknowns at each discrete particle.

In the FSE process, heat is generated from two sources: (1) frictional sliding energy at the billet and the die tool interface and (2) mechanical work caused by the plastic deformation of the billet material. The primary heat sinks are the mass of the apparatus pieces and the active water-cooling system outside the billet container. The total heat source and sink term in Eq. (9) is then comprised of plastic work \(\dot{q}_\mathrm{p}\), frictional sliding energy \(\dot{q}_\mathrm{f}\), and surface convection \(\dot{q}_\mathrm{c}\) as follows:

$$\begin{aligned} \dot{q} =\dot{q}_\mathrm{p}+\dot{q}_\mathrm{f}+\dot{q}_\mathrm{c} \end{aligned}$$
(10)

Note that all the quantities, coordinates, and derivatives in Eqs. (7)–(10) are represented in the Lagrangian reference frame, making the SPH method a truly Lagrangian particle method.

3.3 Material constitutive model

The thermo-elastoplastic flow behavior of the AA1100-O billet is modeled using a Johnson–Cook (J–C) constitutive relation [46], which is expressed as follows:

$$\begin{aligned} \sigma _\mathrm{f}=(A+B\varepsilon _{\mathrm {eff}}^n)(1+C\mathrm {ln} {\dot{\varepsilon }})(1-T^{*^m}) \end{aligned}$$
(11)

where \(\sigma _\mathrm{f}\) is the flow stress, \(\varepsilon _{\mathrm {eff}}\) is the effective plastic strain, \({\dot{\varepsilon }}={\dot{\varepsilon }}_{\mathrm {eff}}/{\dot{\varepsilon }}_0\) is the normalized effective plastic strain rate, and \({\dot{\varepsilon }}_{\mathrm {eff}}\) and \({\dot{\varepsilon }}_0\) denote the effective plastic strain rate and reference total strain rate (typically 1.0 \(\mathrm{s}^{-1}\)), respectively. A, B, C, n, and m are material parameters: A is the initial yield stress, B and n are the strain hardening constants, while C denotes the material sensitivity to plastic strain rate. \(T^*=(T-T_\mathrm{r})/(T_\mathrm{m}-T_\mathrm{r})\) is the homologous temperature: \(T_\mathrm{r}\) is the reference temperature at which the material parameters are evaluated, \(T_\mathrm{m}\) is the material solidus temperature, and T is the instantaneous temperature. The terms in the first and second sets of parentheses in Eq. (11) represent the strain and strain rate hardening effects on the flow stress, whereas the term in the third set of parentheses describes the thermal softening effect.

Equation (11) describes the instantaneous stress required to sustain material plastic deformation or flow. Once the von Mises stress exceeds the flow stress, a radial return algorithm [47] is used to calculate the stress state resulting from the strain increment and the associated plastic strain. After that, the heat source term from plastic deformation \(\dot{q}_\mathrm{p}\) can be calculated as

$$\begin{aligned} \dot{q}_\mathrm{p}=\eta \varvec{\sigma }:\dot{\varvec{\varepsilon }}^p \end{aligned}$$
(12)

where \(\eta \) is the coefficient characterizing the fraction of plastic work dissipated into thermal energy and \(\dot{\varvec{\varepsilon }}^p\) is the plastic strain rate tensor.

Fig. 3
figure 3

Active cooling system on the billet container

Full size image
Fig. 4
figure 4

SPH model for FSE: a Geometries and dimensions of the parts; b assembled model (Units: mm)

Full size image

3.4 Contact condition and heat generation

The FSE process relies on the shear stress from the die tool face to generate heat to soften the billet material and deform it. Thus, a correct die–workpiece contact condition is of critical importance for more accurate material deformation and heat generation. Among the existing contact conditions for solid-phase processing such as FSW, the full sticking condition assumes no slippage between the materials on either side of the contact such that both sides have the same interfacial velocities [11]. However, the sticking condition can only account for the heat generation from material plastic deformation and it usually overestimates the maximum temperature [48]. It was indicated by Gerlich et al. [49] and Schmidt et al. [50] that the contact condition for in solid-phase processing is actually partial sliding/sticking, and frictional sliding energy and material plastic deformation are both responsible for temperature elevation. Using the frictional coefficient and slip rate parameter, many studies have numerically shown that the partial sliding/sticking can yield results that are closer to experimental observations [51,52,53,54].

In this study, a partial sliding/sticking contact condition is realized by Coulomb friction law with Tresca limit, as suggested by Schmidt and Hattel [55]. The Coulomb–Tresca friction law can be naturally implemented in the Lagrangian particle-based model, and it can consider heat generation from both plastic deformation and frictional sliding between workpiece and die. Adopting the Coulomb–Tresca friction law, the friction shear stress \(\varvec{\tau }_\mathrm{s}\) can be written as:

$$\begin{aligned} \varvec{\tau }_\mathrm{s}=\left\{ \begin{array}{l} \displaystyle {\mu _\mathrm{f}{p}\frac{\varDelta \varvec{v}_\mathrm{s}}{\Vert \varDelta \varvec{v}_\mathrm{s}\Vert }} \qquad \mathrm {for} \quad \varvec{\tau }_\mathrm{s}<\varvec{\tau }_{\max }\\ \displaystyle {\varvec{\tau }}_{\max }\qquad \qquad \,\,\,\, \mathrm {for} \quad \varvec{\tau }_\mathrm{s}\ge \varvec{\tau }_{\max }\\ \end{array}\right. \end{aligned}$$
(13)

where \(\mu _\mathrm{f}\) is the friction coefficient, p is the contact pressure determined by the normal slave node force and the slave node area, \(\varDelta \varvec{v}_\mathrm{s}\) is the relative die–workpiece sliding velocity, \(\varvec{\tau }_{\max }=\frac{\sigma _y}{\sqrt{3}} \frac{\varDelta \varvec{v}_\mathrm{s}}{\Vert \varDelta \varvec{v}_\mathrm{s}\Vert }\) is the limiting shear force with \(\sigma _y\) being the equivalent material flow stress. This mechanical contact is modeled using the soft constraint penalty formulation that is suitable for treating contact between bodies with dissimilar material properties (e.g., steel-aluminum) [56]. In this method, each particle on the slave side is first checked for penetration through the master surface. Once penetration occurs, an artificial force proportional to the penetration distance and penalty stiffness is introduced to prevent penetration. The penalty stiffness is the minimum of the master segment and slave node stiffness. To avoid the excessive penetration on the contact side with weaker material, a stability contact stiffness is further added. The details of the soft constraint penalty formulation can be found in [56] and are omitted here.

The energy converted into heat resulting from frictional contact is expressed as

$$\begin{aligned} \dot{q}_\mathrm{f}=\varvec{\tau }_\mathrm{s}\cdot \varDelta \varvec{v}_\mathrm{s} \end{aligned}$$
(14)

which is partitioned to each side of the contact. The frictional sliding energy partition fraction \(f_\mathrm{part}\) is dependent on the thermal diffusivity of the two materials in contact, which is discussed and determined in Sect. 4.2. The thermal contact between SPH particles and solid elements needs further elaboration and is discussed in Sect. 4.1

Note that although many studies showed satisfactory results with the Coulomb–Tresca friction law for both FSE and FSW [42, 54, 57, 58], the friction coefficient \(\mu _\mathrm{f}\) is purely empirical or hypothetical, and its value may affect the accuracy of the simulation result. The Coulomb–Tresca friction law is used in the current SPH model because it is a natural option when the Lagrangian framework is used. To obtain more accurate results, a more physics-based contact model should be developed and incorporated into the SPH model, which is out of the scope of this study.

3.5 Thermal boundary condition

To avoid overheating and preserve the die tool during the FSE process, pumped water provides active cooling to the billet container, as shown in Fig. 3. This thermal dissipation boundary condition is realized in the model by considering a heat flux (\(\dot{q}_\mathrm{conv}\)) on the container surface, which is expressed in a convection equation as:

$$\begin{aligned} \dot{q}_\mathrm{conv}=h_\mathrm{conv}(T_\mathrm{ref}-T)A_\mathrm{s} \end{aligned}$$
(15)

where \(h_\mathrm{conv}\) is the convection coefficient, T denotes the container surface temperature, \(T_\mathrm{ref}\) is the reference temperature in the cooling system which is set as the room temperature in this study, and \(A_\mathrm{s}\) is the container surface area exposed to the cooling system. In accord with the experiment, the convection heat sink term in Eq. (15) is only applied on the free surface of the billet container in the model. The radiation source term is not considered in our model because it is negligible when the active cooling is present.

Fig. 5
figure 5

Finite element meshes and SPH particles for simulation (The numbers of particles or elements are given in parentheses for each part)

Full size image
Table 2 Johnson–Cook model parameters [59,60,61]
Full size table
Table 3 Elastic and thermal properties [62, 63]
Full size table

4 Model setup

The SPH model was implemented in the commercial software LS-DYNA [56]. Explicit dynamic analysis was used to solve this thermomechanically coupled problem with complex contact conditions.

4.1 Model geometry, assembly, and discretization

The model geometry is adapted from the experimental setup in Fig. 1a. The billet ring and backing plate are combined as a single billet container part with 1 mm container wall thickness and 2 mm backing plate thickness. The length of the die tool is 10 mm. The billet container and die tool were assumed to be rigid and discretized into solid elements for contact calculation. The AA1100-O billet part is deformable and will experience large deformation during the process, so it is discretized into SPH particles.

Accurate temperature field prediction is essential for FSE as it directly affects the temperature-dependent material and contact properties. Unlike solid elements that discretize the predefined boundaries and material interfaces of the solid parts, SPH uses particles to represent the changing free surfaces and moving material interfaces at different times. This introduces challenges to identify the SPH boundary at each time step to enable thermal contact between the SPH parts and solid parts. To this end, Fraser [41] recently developed detection methods for particles located at both free surfaces and particle–solid interfaces, which enabled the application of thermal boundary conditions to SPH parts and thermal contact between SPH parts and solid parts. Later, Xu and Wang [64] proposed a virtual SPH method to allow heat transfer between SPH particles and solid parts. In this method, the solid elements are adaptively transformed into virtual SPH particles, which move together with their associated solid elements. Since the virtual SPH particles share the same kernel function with the original SPH particles for the heat transfer equilibrium equation, an authentic thermal coupling between the original SPH parts and new virtual SPH parts can be achieved. The temperature fluctuations in the virtual SPH particles are then inherited by the associated solid elements, be transferred to other solid parts, and interact with the thermal boundary conditions defined on solid parts. In this study, the virtual SPH method proposed by Xu and Wang [64] was used to enable thermal contact between billet SPH particles and rigid parts solid elements.

To reduce the total number of SPH particles, the billet container was split into an inner one and an outer one, and the die was split into a top one and a bottom one. Only the solid elements in inner container and the top die were transformed into virtual SPH particles. As a result, there are five parts in the proposed model with detailed dimensions shown in Fig. 4a. Two longitudinal (z-direction) AA2050-O marker wires with a \(2\times 2\)-particle cross section are inserted: one at the billet center and the other at 1/3 r (radius) away from the billet center. The parts shown in Fig. 4a were assembled in the z-direction tightly with no gaps between them, as shown in Fig. 4b. Both dies rotate counterclockwise around the central axis (viewed in positive z-direction) at a constant rotational speed, and advance in negative z-direction at a constant speed.

For the discretization, an 8-node, 3-D, bilinear element was used for the rigid solid parts. To strike a good balance between the computational cost and numerical accuracy, the billet was discretized into 0.4 mm diameter particles, which leads to 96,768 particles. An element sized in 0.25 mm was used to mesh the solid parts, resulting in 486,984 solid elements, as shown in Fig. 5. Each solid element in the inner billet container and top die were transformed into one virtual SPH particle. The final model featured in this study consists of 198,704 SPH particles, as indicated in Fig. 5.

Fig. 6
figure 6

Sliding energy partition fraction \(f_\mathrm{part}\) between AA1100-O and H13

Full size image
Fig. 7
figure 7

Histories of kinetic energy, internal energy, and their ratio of the billet versus simulation time, with die rotational speed of 300 rpm and die advance speed of 8 mm/min

Full size image
Fig. 8
figure 8

Snapshots of billet material (only) and effective plastic strain during FSE at: a t = 3.0 s; b t = 4.5 s; c t = 7.5 s; and d t = 9.0 s, with die rotational speed of 300 rpm and die advance speed of 8 mm/min

Full size image

4.2 Material properties and model parameters

Because J–C model material parameters for AA1100-O were not available in the literature, the corresponding values for a similar alloy such as aluminum 2024-T351 are often used instead [65]. In our work, the strain rate hardening parameters C and \(\dot{\varepsilon _0}\) were obtained from [59] for AA1100-H12. The thermal softening exponent m was obtained from [60]. A major limitation of the J–C model is that the strain hardening is weakly correlated with the thermal softening. As a result, an unsaturated flow stress will be predicted even at high plastic strain and temperature [43]. To address this issue, the yield stress for AA1100-O in undeformed condition was assigned to initial yield stress A at room temperature, and the strain hardening parameters B and n were set at zero. J–C model parameters for AA2050-O were obtained from [61] also with B and n set to zero. The detailed Johnson–Cook material parameters used in this study are detailed in Table 2.

The rigid billet container and die tool were assigned the properties of H13 tool steel. The elastic and thermal properties of AA1100-O, AA2050-O, and H13 tool steel were obtained from a commonly available material database [62], which are tabulated in Table 3. Because the thermal properties of aluminum alloys are typically sensitive to temperature variation, same temperature-dependent thermal material properties taken from [63] were used for AA1100-O and AA2050-O. The thermal properties of H13 tool steel were kept constant throughout the simulation because they are insensitive to the temperature variation [63]. The density, Young’s modulus, and Poisson’s ratios are assumed to be temperature-independent since they were tested to have small impacts on the simulation results.

An effective friction coefficient \(\mu _\mathrm{f}\) in Eq. (13) was assumed for the contact between the billet and the container or die, the value for which was chosen as constant 0.57 between aluminum and AISI H13 tool steel [66]. All the frictional sliding energy is assumed to be converted into heat. The partition fraction \(f_\mathrm{part}\) of the sliding energy transferred to heat the billet can be calculated as [56]

$$\begin{aligned} f_{\mathrm {part}}=\frac{1}{\displaystyle {1+\sqrt{\frac{(\rho C_p \kappa )_{\mathrm {tool\;material}}}{(\rho C_p \kappa )_{\mathrm {billet\;material}}}}}} \end{aligned}$$
(16)

Based on the material properties in Table 3, \(f_\mathrm{part}\) as a function of temperature is shown in Fig. 6. Since the steady-state temperature at the contact interface is about 400–600 \(^\circ \)C for FSE, \(f_\mathrm{part}\) is assumed a constant value of 0.7 in this study for simplicity. A dissipation coefficient of \(\eta =\) 0.9 was used for the plastic dissipation, implying that 90\(\%\) of the plastic work is converted into heat energy. Because no calibrated heat convection coefficient was available for the cooling system, the convection coefficient on the container surface was chosen as \(h_{\mathrm {conv}} = 1000\; \mathrm {W}/(\mathrm {m}^2\, ^\circ \mathrm {C})\), from the existing data for active liquid cooling on steel [67]. This value has been shown to yield reasonable heat dissipation and temperature profiles.

5 Numerical results and discussion

In this section, the numerical results obtained from the SPH model are presented and validated by the experimental data. Specifically, the material flow, strain, temperature, and extrusion force during the FSE process are investigated.

SPH is much more computationally intensive than conventional mesh-based methods because additional particle–neighborhood searching is needed to identify the kernel function and the bucket-sorting algorithm is required for the node-to-surface contact. To obtain a reasonable computational cost while maintaining accuracy of the results, we used a mass scaling technique to increase the time increment by artificially increasing the material density in the model. To this end, the target minimum step size chosen for the mechanical analysis was \(2\times 10^{-5}\)s, which requires a mass scaling factor of about \(3.91\times 10^5\). A fully implicit time integration was used for the transient thermal analysis, and an unconditionally stable thermal time step was set at \(2\times 10^{-3}\) s. Along with mass scaling, a time scaling technique was also used to further speed up the simulation by increasing the billet advancing and die rotational speeds by a factor of 15. With mass scaling and time scaling, the total simulation time can be reduced to about 52 hours running with 12 CPUs on a 12 core Intel Core i9-9920X with 64 GB memory.

It should be noted that to avoid significant deviation from the actual solution, the acceptable mass scaling factor should maintain the kinetic energy of the deforming SPH particles below 5–10\(\%\) of the internal energy caused by material deformation [68, 69]. Figure 7 shows the histories of the kinetic energy, internal energy, and their ratio for the billet material over the simulation time in the case of 300 rpm die rotational speed and 8 mm/min die advance speed. The figure shows that the kinetic energy is an order of magnitude smaller than the internal energy during the entire simulation. The maximum kinetic/internal energy ratio is about 4.2\(\%\) at the end of the simulation. This result indicates that the mass scaling factor used in this study is reasonable and has no significant effect on the predicted result.

Fig. 9
figure 9

Snapshots of marker material flow with 300 rpm die rotational speed and 8 mm/min die advance speed (inner and outer chambers and bottom die are transparent for a better visualization)

Full size image
Fig. 10
figure 10

Vertical point position tracking in the case of 300 rpm die rotational speed and 8 mm/min die advance speed (bottom view) at a \(t= 0\) s; b \(t= 7.5\) s; c \(t= 29.25\) s; d \(t= 36.75\) s; e cross section of extruded wire at 50 mm from wire start (black phase is 1/3r marker material); f cross section of extruded wire at 100 mm from wire start with 1/3r maker material (black phase is 1/3r marker material)

Full size image
Fig. 11
figure 11

Vertical point position tracking in the case of 300 rpm die rotational speed and 8 mm/min die advance speed, isometric view at a t= 0 s; b t= 7.5 s; c t= 29.25 s; d t= 36.75 s

Full size image
Fig. 12
figure 12

Horizontal point position tracking in the case of 300 rpm die rotational speed and 8 mm/min die advance speed, bottom view at a t= 0 s; b t = 7.5 s; c t= 10.5 s; d t= 18 s

Full size image
Fig. 13
figure 13

Horizontal point position tracking in the case of 300 rpm die rotational speed and 8 mm/min die advance speed, isometric view at a t= 0 s; b t= 7.5 s; c t= 10.5 s; d t= 18 s

Full size image
Fig. 14
figure 14

Illustration of strain measuring locations in the extruded wire: a strain measuring locations; b 1/3r location

Full size image

5.1 Material flow pattern

The billet material flow pattern during FSE is first investigated. Figure 8 presents snapshots of the billet deformation and effective plastic strain at different times predicted by the simulation with die rotational speed of 300 rpm and die advance speed of 8 mm/min. As shown in the figure, the rotating die gradually pushes against the billet with advance speed. The compressive extrusion forces and high shear forces from the die face help heat and deform the billet top in the high shear region. The plastic strain was first developed along the circumference of the top billet face where shear stress is large, as shown in Fig. 8a, b. With time increased, the billet temperature increased with reduction in material flow stress, and more billet materials in the high shear zone were plasticized, as shown in Fig. 8c, d. The softened and plasticized material was then extruded through the die orifice to form a wire. The total length of the simulated extruded wire is 795.5 mm for a 51-second processing time. Even though the billet material has undergone a large plastic deformation, no numerical problems were encountered because SPH is meshless. Moreover, since a penalty-based contact between particles and solid elements with restricted penetration check was employed, particle penetrations into the die were prevented. No particle splashing or obvious large gaps between particles were observed, indicating appropriate selection of time steps, material properties, and thermal boundary conditions for both mechanical and thermal analyses.

Figure 9 provides snapshots of the center and 1/3r marker material flows at different times in the case of 300 rpm die rotational speed and 8 mm/min die advance speed, from which distinct material flow pattern can be observed for different marker materials. The marker material at the billet center (yellow) was directly extruded from the die orifice without much lateral deformation. This is because they do not experience shear force applied by the rotating die. In contrast, the marker material at 1/3r from the billet center (green) was pushed against the die with pressure. With the die rotation, high shear force created by frictional contact between billet and die was applied to the 1/3r marker material, forcing it to rotate with the die. As the billet container advanced, the 1/3r marker materials close to the die orifice were extruded first, and the remaining marker materials continued to rotate until they were extruded.

Although the material flow pattern can be explicitly tracked by displaying the particle locations as shown in Fig. 9, only scattered data points are available, the number of which is restricted by the spatial density of the SPH particles. To obtain a smoother trajectory of the material flow, a cubic spline data interpolation is used to interpolate an additional 1000 query data points among the existing data points in the 3-D domain. The red dots in Fig. 10a to d represent the locations of the 1/3r vertical materials from the bottom view at \(t=\) 0, 7.5, 29.25, and 36.75 s, respectively. In the right-hand-side plot of each subfigure, both the original tracked points (red hollow circles) and the cubic spline fitted curve (black solid line) are given in the 2-D coordinate system. The coordinates of each point in these plots indicate the actual particle location in the billet and extruded wire. The same information in Fig. 10 is depicted in Fig. 11 but in a 3-D isometric view. From these figures, a clearly helical material flow pattern can be observed for the vertical 1/3r marker material. The prediction by the numerical model also indicates the helical material flow inferred by experimental observation of different cross sections of the extruded wire, as indicated by the black phases in Fig. 10e, f.

In addition to tracing the material vertically in the billet, the materials were also tracked horizontally and interpolated at different times. To this end, two lines of material, in the y direction and in the second SPH layer from the billet top surface, were selected, as shown in Figs. 12a and 13a. The material flow patterns at \(t = 7.5, 10.5\), and 18 s are plotted in Figs. 12b–d and 13b–d, respectively. Since there are two lines of horizontal materials in this case, the fitted spline curve is based on the average location of the two lines. The helical material flow pattern is observed again but is more obvious in these figures. One can also see that, although the particles close to the billet rim are highly deformed, they are unlikely to be extruded because they are far from the die orifice. The flow trajectories of these particles can help to explain the experimentally observed “dead zone” at the billet container corner where the material is difficult to be extruded.

Fig. 15
figure 15

Strain comparison between simulation results and experimental data, 300 rpm die rotational speed and 8 mm/min die advance speed

Full size image

5.2 Strain validation

A more quantitative comparison of the material deformation predicted by the simulation and experimentally observed was carried out by looking at the strain measurements. In the experiment, the strains at 1/3r of the wire cross section at 60, 100, 200, 300, 400, and 500 mm from the wire start were estimated based on the spiral shape of the 1/3r marker material at these locations [44]. The cylindrical coordinate system was used for the estimated experimental strains. To compare these to results from simulation obtain the same strain data must be converted into cylindrical coordinates. The strain tensors in the Cartesian coordinate system were extracted from the same locations in the wire, as indicated by the red dashed boxes in Fig. 14a. At each location, a 5-mm-long wire segment is selected and the strain tensors of all the particles at 1/3r of the wire segment are extracted, as indicated by the small white window in Fig. 14b. After the strain tensor and the displacement history are obtained for each particle within the small white window, a coordinate system transformation is performed to transfer the strain data from the Cartesian to cylindrical systems. After calculating all the strain tensors, the averaged values at each location are used to represent the simulated 1/3r strain in the cylindrical system.

Fig. 16
figure 16

Temperature history comparison between simulation and experiment: a 300 rpm die rotational speed with 8 mm/min and 12 mm/min die advance speeds; b 8 mm/min die advance speeds with 300 rpm and 150 rpm die rotational speeds

Full size image

Comparisons of longitudinal strain (\(\varepsilon _{zz}\)), tangential strain (\(\varepsilon _{\theta \theta }\)), radial strain (\(\varepsilon _{rr}\)), and the sum of strains (\(\varepsilon _{rr}+\varepsilon _{\theta \theta }+\varepsilon _{rr}\)) in the case of 300 rpm die rotational speed and 8 mm/min die advance speed obtained from simulation and experiment are shown in Fig. 15, where z, \(\theta \), and r denote the longitudinal, tangential, and radial directions, respectively. It is shown that the longitudinal strains at different locations match well, with magnitudes around 4. There are some mismatches in radial and tangential strains, especially at the locations closer to the wire start. This may be because at the start of the process, a steady state has not yet been reached and material has not been extruded from the entire plane of the billet. Nevertheless, the strain predictions are overall in good agreement with the experimental data, demonstrating the validity of the SPH model for predicting large strains during FSE.

5.3 Temperature validation and investigation

The proposed model was further validated by comparing the temperature data obtained from the simulation and experiment at three different processing parameters. The temperature variation in the experiment was recorded by a K-type thermocouple inserted 0.5 mm behind the die face and 5 mm from its central axis, as indicated in Fig. 1. Figure 16a shows the comparison results with 300 rpm die rotational speed, 8 mm/min (black) and 12 mm/min (red) die advance speeds, while Fig. 16b shows the ones with 8 mm/min die advance speed, 300 rpm (black) and 150 rpm (red) die rotational speeds. In both figures, the experimentally measured temperatures as a function of processing time are plotted as the dashed curves, and the simulated temperature variations at the same location relative to the die face are depicted as the solid curves.

As indicated in Fig. 16, the temperature variations on the die face predicted by the SPH model are generally in good agreement with the experimental measurement for all three sets of process parameters. The SPH model can capture the rapid temperature increase in the early extrusion stage (between 0 and 15 s) attributable to high material plastic deformation and frictional sliding energy. The temperatures increase slower in the middle and late stages (after 15 s), which is a result of material softening and more heat loss by active cooling at high temperature. When the die rotational speed is fixed at 300 rpm, a faster die advance speed leads to a faster temperature elevation on the die face as shown in Fig. 16a. This is because the billet material is plasticized more and extruded faster with a higher die advance speed, which leads to more heat generation during the same amount of time compared to the case with lower die advance speed. When the die advance speed is fixed at 8 mm/min, a faster die rotational speed results in a faster temperature development as shown in Fig. 16b. This is due the fact that the frictional shear stress applied on the billet is severer with a higher die rotational speed. This higher shear stress plasticizes the billet and generate heat faster compared to the case with slower die rotational speed.

Fig. 17
figure 17

Predicted temperature profiles on the central cross section at a 4.5 s, b 9 s, c 13.5s, d 22.5s, e 31.5s, and f 40.5 s, 300 rpm die rotational speed and 8 mm/min die advance speed

Full size image

Figure 16 shows that the simulation underestimates temperature in the late stage. This discrepancy may arise for two reasons: (1) the current model uses an approximated heat flux boundary condition on the container surface; (2) in the actual experiment, sometimes “flash” was extruded through the small gap between the outer surface of the die and the inner wall of the billet “ring.” The friction between this flash and the rotating die generated additional heat. A thermal boundary calibrated to the active cooling system could improve the accuracy of the predicted temperature variation, but that was not pursued here. This temperature validation indicates that the heat generation, conduction, and dissipation are reasonably captured and accounted for in the thermomechanically coupled SPH model. A high-fidelity temperature field again yielded more accurate material mechanical properties for use in simulating the FSE process.

Fig. 18
figure 18

Simulated temperature profiles on billet container surface and bottom die at a 3 s, b 9 s, c 18 s, and d 26 s, 300 rpm die rotational speed and 8 mm/min die advance speed

Full size image

5.4 Extrusion force validation

Figure 17 presents temperature distributions on the central cross section of the simulation domain at different times in the case of 300 rpm die rotational speed and 8 mm/min die advance speed, which are difficult to measure by experiment. The temperatures in the inner billet container and the top die are presented only on virtual SPH particles for better visualization. In the early stage (\(t = 4.5 \hbox {s}\)), most of the simulation domain remains at room temperature, while heat generation and conduction from the high-shear billet–die contact region begin. With increased processing time, heat is continuously generated by the frictional force and plastic deformation and is conducted to the rest of the domain, as shown in Fig. 17b–f. Our results show that a higher temperature always occurs at the billet–die contact. In the early stage, the highest temperature occurs close to the highly confined corner, as shown in Fig. 17b. As the temperature increases, the outer billet container dissipates more thermal energy, cooling down the billet particles that are close to the chamber. Temperatures gradually increase closer to the extrusion hole, as shown in Fig. 17c–f. Moreover, temperature transitions between the SPH billet and virtual SPH particles, and between virtual SPH particles and solid parts, are smooth, indicating that the virtual SPH particles effectively bridge the heat transfer between SPH parts and solid parts. The effectiveness of heat transfer by virtual SPH particles is observed in Fig. 18, which shows the simulated temperature profiles on the outer billet container and the bottom die at different times.

The temperature validated SPH model was then used to predict thermal histories at the point around the die hole and the central point on the billet bottom, whose locations are indicated by the white boxes in Fig. 17a. Figure 19 gives the predicted results with all three process parameters. It is shown that the cases with faster die rotational and advance speeds have faster temperature elevation and higher final temperature for points at both die hole and billet bottom. Temperatures at the die hole increase faster to higher values compared to the ones at the billet bottom. The final temperature at the die hole for the case of 300 rpm die rotational speed and 12 mm/min die advance speed reaches to 626.6 \(^\circ \)C, which is just below the solidus temperature of AA1100-O (643 \(^\circ \)C). These numbers for the cases of 300 rpm die rotational speed and 8 mm/min die advance speed, and 150 rpm die rotational speed and 8 mm/min die advance speed are 560.4 \(^\circ \)C and 498.6 \(^\circ \)C, respectively. The final temperatures at billet bottom for all three cases range from 479 to 525 \(^\circ \)C.

Fig. 19
figure 19

Temperature predictions at the die hole and billet bottom

Full size image
Fig. 20
figure 20

Extrusion force comparison between simulation and experiment: a 300 rpm die rotational speed with 8 mm/min and 12 mm/min die advance speeds; b 8 mm/min die advance speeds with 300 rpm and 150 rpm die rotational speeds

Full size image

5.5 Extrusion force validation

Because constant die advance speeds were considered in FSE, the extrusion force varies. Therefore, extrusion force is also validated since it directly relates to the thermomechanical behavior of billet material. Figure 20a shows the comparison results with 300 rpm die rotational speed, 8 mm/min (black) and 12 mm/min (red) die advance speeds, while Fig. 20b shows the ones with 8 mm/min die advance speed, 300 rpm (black) and 150 rpm (red) die rotational speeds. Dashed curves indicate the experimental data and the solid curves denote the simulation results. As shown in Fig. 20, typical elastic stages are captured by the model, which agrees well with the experimental data. After 5 s, the extrusion force starts to decrease as the temperature increase to soften the billet material. This decrease in extrusion force obtained from simulation matches the experiment well, indicating the thermal field and material softening behavior were predicted accurately. A mismatch is observed in the peak load between 18 and 25 s for the case of 300 rpm die rotational speed and 8 mm/min die advance speed, where a large decrease in extrusion force occurs in the experimental data. A possible reason for the mismatch is that the contact between the sample and the backing block was imperfect in the experiment, while it is assumed to be perfect in the SPH model. Despite that, an overall good agreement between model predictions and experimental measurements was achieved.

6 Conclusions

In this work, a meshfree SPH model was developed for predicting the extreme thermomechanical conditions during the FSE process. In this model, the billet container and die tool were considered as rigid bodies, while the billet material was modeled as a thermo-elastoplastic solid described by Johnson–Cook material model. A virtual SPH approach was used to enable heat transfer between SPH parts and rigid solid parts. Contact between the billet and the solid parts was modeled by the Coulomb–Tresca friction law, and a heat dissipation thermal boundary condition was applied on the surface of the outer billet container. The developed SPH model was validated by comparing the numerically predicted material flow, strain, temperature, and extrusion force results with experimental observations at three different processing parameters. The comparisons indicated a good agreement between the numerical results and experimental data, demonstrating that the thermomechanical phenomena during FSE have been accurately simulated by the SPH model. Our numerical results reveal that the temperature elevates faster with faster die rotational speed and die advance speed due to more heat generation. The extrusion force is proportional to the die advance speed per rotational evolution.

The validated SPH model can be used to investigate the influences of various process parameters (e.g., billet material, die tool geometry, die rotational and advance speeds) on the field variables of the extrudate (e.g., material flow, temperature, strain, stress). Based on the field variables’ evolution histories, the material microstructure (e.g., grain size, microhardness, texture) evolution can be better understood through a weakly coupled post-processing metallurgical simulation. This model can also be used to optimize process parameters. A surrogate model could be constructed and calibrated to the input-response data generated by the SPH model. The trained surrogate model could be used to determine the optimum FSE process parameters to achieve certain desirable properties in the final extrudate.

In the current model, the Coulomb–Tresca friction law with a constant friction coefficient was used to describe the contact between the billet and rigid parts, and the phenomenological Johnson–Cook model was used for the thermomechanical response of the billet material. Future development of the SPH model should incorporate physics-based contact models by investigating the material attributes from lower length scales. Metallurgical aspects such as the effects of initial grain size, grain growth and recovery, and recrystallization can be included in the material constitutive model to fully describe the phenomena occurring during FSE. Moreover, detailed comparisons between SPH and grid-based FEM for FSE in terms of numerical accuracy and computational cost can be conducted in our future work.