Abstract
A meshfree modeling technique of material flow in the three-dimensional multiphysics thermomechanical friction stir welding process is presented. In this numerical model, the discretization in space is derived by the Element-Free Galerkin method using a Lagrangian meshfree convex approximation. The discrete thermal and mechanical equations are weakly coupled as the time advances using a forward difference scheme. A mortar contact algorithm is employed to model the stirring effect and heat generation due to frictional contact. Heat conductance between contacting bodies is considered as a function of contact pressure. A two-way adaptive procedure is introduced to the coupled thermomechanical system to surpass potential numerical problems associated with the extensive material deformation and spatial discretization. In each adaptive phase, a consistent projection operation utilizing the first-order meshfree convex approximation is performed to remap the solution variables. Finally, a three-dimensional multiphysics thermomechanical coupled friction stir welding problem is analyzed to demonstrate the effectiveness of the present meshfree numerical procedure.
1. Introduction
Friction stir welding (FSW) is an innovative welding method [1] that is viable for jointing aluminum alloys, cooper, magnesium, and other low-melting-point metallic materials. It combines frictional heating and stirring motion to soften and mix the interface between two work pieces yielding a solid and fully consolidated weld. Since FSW is a solid state jointing process, a high quality weld can be achieved with the absence of solidification cracking, porosity, oxidation, and other defects typical to traditional fusion welding. FSW offers high levels of repeatability, limited energy consumption, and ease of automation and thus gains its popularity in aerospace, automotive, railway, and nuclear industries. The FSW process consists of four basic phases, namely, the plunging, stirring, welding, and retraction. The plunging phase starts with plunging the rigid cylindrical spinning tool into the joint line until the shoulder contacts the top surface of work piece. In stirring phase, the heat is generated by means of work from frictional contact and material plastic deformation. This heat is dissipated into the welding zone and results in an increase of temperature and material softening. In particular, the heat generated from frictional contact plays an important role in determining the material rheological behavior and the success of the deposition welding process. The welding phase is initiated by moving either the tool or the work piece and enables the materials of two work pieces to mix together. Finally, the welding process stops and the tool retracts from the work piece. These four FSW phases constitute complicated multiphysics thermomechanical conditions which are very hard to determine experimentally [2]. Although many successful experimental investigations have already been conducted to the adjustment of tool profile and input weld parameters such as tool speed, feed rate, and tool depth, several aspects of the FSW are still poorly understood and require further study [3]. In particular, understanding the temperature distributions in the work piece during FSW process is an essential subject [4] due to its effects on the material flow, grain size, residual stresses, and, subsequently, the strength.
Finite element modeling provides an effective approach to better understand and visualize the influence of tool profile and input weld parameters in the welding process. It is recognized that a correct numerical model of the FSW process should avoid any unnecessary assumptions. A list of requirements [3] for multiphysics FSW analysis code should include () rotational boundary condition; () frictional contact algorithm; () support for very high levels of deformation; () elastic-plastic or elastic-viscoplastic material models; and () support for complex geometry. The main difficulty in finite element modeling of FSW process consists in dealing with high levels of deformations involving in the complex material flow due to frictional heating and stirring motion. The Eulerian description in Computational Fluid Dynamics (CFD) models [5–7] can be easily adapted to FSW process to circumvent the mesh distortion problem encountered in the Lagrangian formulation. In this scenario, the material is analyzed as a viscous incompressible non-Newtonian fluid flowing across the Eulerian mesh and interacting with the rotating tool. Since Eulerian representation of material flow does not capture free surfaces, additional algorithms [8] are required to introduce free surfaces and thus the sticking or sliding conditions. The numerical complexity in Eulerian formulation is also concerning the needs to model the welding heat source [9] in FSW process. The heat dissipation due to frictional contact necessitates the use of the analytical [10, 11] or particular thermal models [12] for the heat source calculation that considers the rotational speed of the tool and shear stress at the interface of tool and work piece. Recently, a solid mechanics-based Eulerian model [13] was developed as an alternative to describe the heat source and correspondingly the mechanical power in the Eulerian approach. Using this modeling technique, the relationship between the rotating speed and mechanical power can be well characterized. Another numerical challenge in Eulerian modeling of FSW process is the inherent complexity in advecting material histories. When the Eulerian description of FSW process is modeled using the flow formulation, difficulties arise in the followed elastic unloading and residual stress analyses [2, 3]. An improvement of mesh distortion problem in finite element modeling of FSW process can be obtained by an -adaptive remeshing scheme such as the Arbitrary Lagrangian Eulerian (ALE) method [14–16]. The ALE algorithm advances the mesh independently with material flow and makes it possible to take into account the movements of free surfaces while reducing mesh distortions. An 8-node trilinear element with reduced integration and hourglass control is usually utilized to avoid the volumetric locking problem in von Mises materials. Since advancing the ALE mesh does not involve additional degrees of freedom, it is found to be efficient in the steady-state analysis [17] of FSW process. Nevertheless, the ALE approach may fail in the transient analysis when the tool travels a long distance and the high-quality mesh cannot be maintained in the stirring zone. Unlike the -adaptive remeshing technique in ALE formulation, the global or combined rh-adaptive remeshing schemes [18] can be employed to effectively minimize the mesh distortion in the stirring zone. Buffa et al. [19] presented a thermomechanically coupled finite element model based on a combined rh-adaptive scheme for the FSW analysis. The main difficulty of three-dimensional mesh generation comes from the existence of silvers (highly distorted tetrahedral) with extreme dihedral angles. Although the Delaunay refinement scheme can remove silvers by adding some amount of Steiner point [20], it is still an open problem to determine a nontrivial lower bound on the dihedral angles of the output tetrahedral. In addition, the rigid-viscoplastic material assumption demands a postprocessing to extract the nodal temperature history for the calculation of residual stress field induced by the thermal cycle and the mechanical action applied by the tool.
An alternative approach to model large deformation problem is the use of meshfree methods. Meshfree methods based on the Galerkin approach have been shown to effectively deal with the von Mises and pressure dependent materials in large deformation analyses [21–24]. The Moving Least-Squares (MLS) approximation in Element-Free Galerkin method [25] and the Reproducing Kernel (RK) approximation in Reproducing Kernel Particle [26] are two commonly used approximations in meshfree methods. The MLS approximation was achieved by performing a position dependent weighted least squares fit and the RK approximation was achieved by introducing a correction function to the Smooth Particle Hydrodynamics (SPH) approximation [27] to impose the reproduction of basis functions. However, both approximations do not bear the Kronecker-delta property, which is a convenient property for imposing the Dirichlet boundary conditions. A different way of constructing an effective meshfree approximation for the solid mechanic analysis is to incorporate the convex approximants [28, 29]. The meshfree convex approximation guarantees the unique solution inside a convex hull with a minimum distributed data set and admits a weak Kronecker-delta property at the boundaries and therefore avoids the special treatments [30] on the Dirichlet boundaries. Wu and Koishi [28] provided a unified approach that can generate specific convex approximations and reproduce several existing meshfree approximations which are referred to as the Generalized Meshfree (GMF) approximations. From the dispersion analysis, it has been shown [31] that meshfree convex approximation exhibits smaller lagging phase and amplitude errors than conventional MLS and RK approximations in the full discretization of the wave equation. The eigenvalue analysis also reveals that a larger critical time step can be used in the explicit dynamic analysis. The meshfree convex approximation has been introduced to the low-order finite element to yield an inf-sup stable Meshfree-Enriched Finite Element (ME-FEM) [32] for the particulate rubber composites in the near-incompressible micromechanical analysis [33]. Recently, an immersed meshfree Galerkin formulation was proposed by Wu et al. [34] using the nonconforming meshfree convex approximation to solve the interface constraint problems in composites and was shown to satisfy an optimal error estimates in the energy norm. Most recently, a general framework of the meshfree-enriched multiscale finite element formulation was developed by Wu and Hu [35] for the acoustic wave propagation analysis. In their multiscale finite element model, the fine-scale meshfree convex approximation satisfies the homogenous Dirichlet boundary condition and requires no a priori knowledge of fundamental solution for approximating the oscillatory type of Helmholtz equations. While those meshfree methods using meshfree convex approximation provide high fidelity predictions of material responses in the solid mechanics analyses, their robustness has not yet been fully investigated in the severe deformation analysis of thermomechanical problems such as FSW process. Although an alternative solution using SPH method [36] has demonstrated the meshfree capability in modeling severe material flow in FSW process, the fundamental difficulties associated with the tensile instability [37] and the spurious modes [38] remain challenging numerical issues in solid applications. In particular, the explicit dynamic nature of SPH formulation and the nontrivial essential boundary condition treatments prohibit the adoption of an implicit version for the postwelding residual stress analysis.
The objective of this study is to present a robust numerical procedure based on the Element-Free Galerkin method [25] using the Lagrangian meshfree convex approximation and a two-way adaptive procedure for the multiphysics simulation of FSW process. The reminder of the paper is organized as follows: the governing equations in FSW problem and their variational formulations are given in Section 2. Details on the spatial discretization of the variational formulations using Lagrangian meshfree convex approximations are given in Section 3. The numerical treatment of contact condition is also discussed in the same section. Section 4 describes the time discretization for the coupled thermomechanical equations. The computational strategies for the two-way adaptive procedure and the remap algorithm are given in Section 5. Finally, a numerical example is presented in Section 6 and conclusions are made in Section 7.
2. Governing Equations
2.1. Thermal Problem
We consider the transient heat transfer response of a FSW work piece in a three-dimensional case. We assume that the domain of the work piece is a bounded polygon with disjointed boundary . The notation describes a Dirichlet boundary imposed by a temperature . is the Neumann boundary prescribed by a heat flux. We also assume that and do not vary with time. The boundary denotes the contact surface with a thermal exchange due the conduction between the tool and work piece. The boundary is the surface with the thermal exchange due to convection and radiation. We further assume that the heat generation in the work piece is only due to the plastic deformation and frictional contact between tool and work piece. Giving an internal heat generation rate per unit deformed volume from the plastic deformation, the strong form of the thermal energy conservation equation readssubject to boundary conditionsand initial condition at time where is the mass density and is the heat capacity. Equation (2) is known as the Fourier law [39] with denoting the isotropic thermal conductivity. is the gradient operator with respect to current position and “·” denotes the divergence operator. The Taylor-Quinney [40] coefficient in (3) takes into account the fraction of heat generated by the plastic deformation energy dissipation. and are the deviatoric parts of Cauchy stress and the rate of plastic straining, respectively. in (4) is the temperature imposed on the Dirichlet boundary. in (5) is the normal heat flux imposed on the Neumann boundary. , , and are heat transfer coefficients for conduction, convection, and radiation, respectively. The notation defines the outward unit normal vector on . The second term on the right-hand side of (6) represents the rate of frictional energy dissipation where is the fraction of heat generated by the frictional contact. is the Cauchy contact traction and is the contact slip rate which represents the jump in velocity across the contact. is the ambient temperature and is the temperature of the tool. in (8) is the position vector in the reference configuration. It is also necessary to match the initial condition with the Dirichlet boundary condition at time :
The variational formulation of the thermal energy conservation equation can be written to find the temperature field such that for arbitrary variation the following equation is satisfied:Integration by parts the second term of left-hand side in (10) leads to
2.2. Mechanical Problem
The motion of the FSW work piece is governed by the equation of motion with the prescribed boundary and initial conditions. The mechanical problem is given as follows:subject to Dirichlet and Neumann boundary conditionstogether with unilateral contact conditionsand Coulomb friction conditionsand initial conditionswhere is the body force vector and is the Cauchy stress obtained from the constitutive law. The notation describes a Dirichlet boundary imposed by a displacement and is the Neumann boundary prescribed by a surface traction with . The notation in (14) is the gap function defined in Section 3.2. denotes the outward unit normal to contact surface . is the tangential component of contact traction . Equation (15) states that the vector of slip rate should be in the direction of tangential traction. is the temperature dependent friction coefficient. stands for the Euclidean norm. and are the initial displacement and velocity, respectively.
Subsequently, the variational equation for the mechanical problem in FSW process is formulated using the integration by part to find the displacement field , such that for arbitrary variation , the following equation is satisfied:
Now the FSW problem is stated by coupling the mechanical weak form in (17) with the thermal weak form in (11) and subjecting to the prescribed Dirichlet boundaries and initial conditions.
3. Spatial Discretization
3.1. Thermomechanical Equations
The standard Galerkin method is formulated on a finite dimensional space employing the thermal weak form of (11) to find such thatwithwhere and is an index set. are meshfree shape functions. In order to prevent the tensile instability caused by the Eulerian kernel, the Lagrangian kernel approach [21, 22] is considered in this development. In Lagrangian formulation, the motion of a material point originally located at a position vector in the reference configuration is described by a mapping , where is the spatial position of the material point at time . Therefore, the material displacement can be defined by
Consequently, with meshfree discretization, discrete points that carry the primary unknown variables are attached to the same set of material points throughout the course of deformation. Under this consideration, the node set is the set of nodes defined in the reference configuration. In practice, the set of meshfree nodes can be taken from the finite element nodes created by a finite element mesh generator. We also let be the interior of the supp . We assume that is star-shaped with respect to a ball and there exists a constant such thatLet denote the nodal support radius and assume the following overlapping condition [29]:Often, a Lagrangian meshfree shape function is associated with a meshfree node and nodes are distinct.
Similarly, we have the mechanical weak form of (17) to find such thatwithThe material displacements are also approximated using the Lagrangian meshfree shape functions [21] aswhere is called the “generalized” displacement of node . In other words, the material displacement is considered as an interpolant of in a generalized sense. In general, conventional meshfree approximations are not interpolants; that is, . Subsequently, the temperature field is approximated using the Lagrangian meshfree shape functions and is given byWe also have in the approximation of semidiscrete thermal equations. For this reason, various numerical treatments using conforming or nonconforming methods can be chosen to impose the Dirichlet boundary conditions of the coupled thermomechanical problem. One may refer to [30] for a bibliographical study of those resolution techniques. In this study an alternative meshfree approximation that restores a weak Kronecker-delta property at the boundary, the meshfree convex approximation [28, 29], is utilized to allow the direct treatment of Dirichlet boundary conditions for FSW problem. We employ the GMF [28] method to obtain the meshfree convex approximation. The first-order convex GMF approximation is constructed using the inverse tangent basis function and the cubic spline window function is chosen to be the weight function in GMF method. Giving a convex hull of a node set defined bythe GMF method is introduced to construct a convex approximation of a given (smooth) function in the form of (25) and (26) such that the shape function satisfies the following linear polynomial reproduction property:With the meshfree convex approximation, we can define a -conforming subspace for the approximation of displacement field to be . Accordingly, the same approximation space definition applies to the temperature field. More detailed information about the derivation of GMF method and the corresponding mathematical properties can be found in [29].
Substituting (25) and (26) into (18) and (23) using meshfree convex approximation, the semidiscrete equations of the coupled thermomechanical problem can be expressed by the following algebraic equations:wherewhere denotes the material gradient operator. For the implementation purpose, all terms associated with the volume integration are evaluated on the reference configuration and the surface integration terms are computed on the current configuration. Using Nanson’s formula, the Cauchy stress that appears in the internal force term is transformed to the first Piola-Kirchhoff stress tensor . Note that we have used the chain rule for the spatial gradient operation in the above derivations. and are vectors of generalized nodal values for the interior nodes to be solved for the temperature and displacement fields, respectively. With the meshfree convex approximation, the unknown nodal vectors of temperature and displacement fields for boundary nodes are denoted by and , respectively. For convenience, the finite element mesh is taken as the integration cells [21, 22] for the domain integration. Each integration cell occupies an initial volume needed in the domain integration and the one-point integration rule is used for each integration cell in the computation. The integration cells also provide a set of boundary nodes information for the contact traction calculation as described in the next subsection.
3.2. Thermal Contact Model
The contact traction in (33) and (36) is computed by imposing the contact constraints through the mortar contact algorithm, which is briefly summarized as below. The interested readers may refer to [41, 42] for the implementation details.
Let and denote the current configuration of deformable slave contact surface and rigid master contact surface, respectively. Choosing any contact point , the gap function in (14) is defined with respect to asUsing the above definition, the two-body mortar contact virtual work defined on the nonmortar slave contact surface can be written bywhere is the mortar multiplier representing the Cauchy contact traction in (6) and (14). and are the time dependent displacement fields of slave and rigid master contact surfaces, respectively. The discrete version of mortar contact virtual work is defined by introducing appropriate shape function expansions to the contact surface fields and mortar multiplier. The mortar contact interpolation is performed as follows:where and are numbers of nodes on the slave and master surfaces, respectively and and are finite element interpolation functions defined on the current configuration of piecewise contact segments of and master contact surface , respectively. The associated discrete nodal values are denoted by , , and . is considered as an appropriate projection from slave contact point to master surface. Substituting (47) into (46) leads to the following discrete form:where and are referred to as mortar integrals:In order to impose the normal and tangential contact constraints, the mortar multiplier is decomposed bywhere is unit averaged normal vector at slave node . The definition of on a discretized contact surface can be found in the literature [41]. Substituting (50) into (48) yieldsAccording to (45), the discrete mortar normal contact gap and tangential slip gap at slave node are given [41] bywhere is unit tensor. The normal part of mortar contact traction is subject to discrete form of Kuhn-Tucker conditions from (14) viaImposing (54) by penalty regularization leads to the following contact pressure:where is normal penalty parameter.
In mortar contact algorithm, the penalty regularization [43] of augmented Lagrangian scheme is often used to impose the frictional contact constraints, which requires the computation of the incremental tangential slip gap. From (53), an incremental tangential slip gap which is frame indifferent is given [41] byAccordingly, a trial state-return map strategy is employed to determine the Coulomb frictional tractions. In the tangential contact direction, we first assume no slip from time to , and trial frictional traction can be expressed bywhere is frictional penalty parameter. The frictional contact traction is then corrected based on the slip condition at to givewhich can be used to describe the discrete stick/slip conditions of (15). When the tool surface is in contact with the work piece, the standard Fourier law cannot be used to fully describe the heat transfer phenomena. This is because the contact surfaces do not physically match perfectly leading to the gas trap between the contact surfaces. In general, this heat resistance decreases as the contact pressure increases. For this reason, the heat conductance is computed as a function of contact pressure in this study.
4. Time Discretization
The highly coupled and nonlinear system in the FSW thermomechanical equations is difficult to solve by simultaneous time-stepping algorithm. The large and unsymmetric system in the fully coupled thermomechanical equations inevitably involves the convergent problem and is expensive to solve, particularly in the presence of large deformation, severe contact conditions, and contact-induced thermal shock. In the application of interest, explicit and staggered time-stepping schemes [40, 44] are considered in this study.
The thermal equation in (29) is marched through time using the forward difference algorithm [39] to giveIt also suffices to integrate the mechanical equation (30) by the central difference integration algorithm to yieldwhere the capacity matrix and mass matrix are advantageously replaced by the lumped matrices ad in thermal and mechanical equations, respectively, in the explicit analysis. Using the Lagrangian meshfree shape functions, the position vector for integration point is updated bywhere the integration point position is initially located at the centroid of the integration cell and moves with material flow to its current position . In contrast to the finite element method, the current position of integration point needs not reside in the integration cell. Following the notation in (26), is defined to be the “generalized” current position of node . Subsequently, the deformation gradient at is computed using (39) and (61) to giveFor the computational efficiency in explicit time integration method, the material derivatives of meshfree shape functions are always computed and stored at the first Lagrangian time step and we reuse them during the time stepping.
In the staggered time-stepping scheme, the thermomechanical system is partitioned into two phases, the isothermal mechanical phase and the rigid conduction phase, during each time increment. The isothermal mechanical phase assumes a constant temperature during the mechanical computation and the rigid conduction phase considers the constant heat generation in the thermal computation at fixed configuration. A numerical stability requirement limits the maximum time increment in the explicit method for each phase. The overall stable time increment is then defined as the smaller of the two phases. In practice, the mechanical equations always set the critical time step for stability due to much smaller time scale associated with the mechanical problem.
5. Two-Way Adaptive Procedure and Remap Algorithm
As mentioned earlier in Introduction, the major difficulty arises from the numerical modeling of FSW process when it comes to large deformation simulation of the complex material flow in frictional heating and stirring motion. Strictly Lagrangian approach based on a fix mesh in finite element method experiences difficulty in dealing with mesh entanglement due to unconstrained material flows. Although the meshfree Galerkin method using Lagrangian kernel helps improve the mesh entanglement problem in finite element method, it is prohibitively extending the range of meshfree applicability to model severe deformation that goes beyond the Lagrangian description; that is, the discretized deformation mapping in (63) ceases to be injective:
One way to sidestep this numerical difficulty is to consider the reordering of the neighboring nodal information [45, 46] such that a reconstruction of the Lagrangian meshfree shape functions is allowed and the possibility of nonpositive Jacobin determinant in (63) is suppressed. Another way to evade the nonpositive Jacobin determinant problem is to adopt an adaptive procedure similar to the combined rh-adaptive remeshing or global remeshing techniques in the finite element method. When compared with the approach of reordering of the neighboring nodal information [45, 46], the adaptive method is able to refine the nodal density and generate accurate free surfaces for a better simulation in manufacturing problems. In this study, the adaptive method based on the concept of global remeshing is adopted and an anisotropic unstructured tetrahedral mesh is pursued to model the complex geometries evolved in the formation of stirring zone. A set of meshfree nodes is extracted from the unstructured tetrahedral mesh created by a mesh generator and used for the reconstruction of the Lagrangian meshfree shape functions. The set of contact boundary nodes can also be obtained from the generated tetrahedral mesh. While most research in adaptive error estimation has focused on the development of an accurate error estimator on linear analysis, formidable difficulties still remain in nonlinear analysis. For example, the well-known Z-Z error [47] developed based on superconvergence properties for linear problems cannot quantitatively be estimated in nonlinear problem. Therefore Z-Z error estimator can only be considered as an error indicator in nonlinear problems. In engineering practice the use of an efficient error indicator turns out to be more desirable particularly in the explicit time integration method. In this study, a pointwise error indicator based on the shear deformation is used to trigger the adaptive procedure for the von Mises material. The global mesh generation in this study is comprised of an anisotropic surface Delaunay triangulation [48] and an Advancing Front tetrahedral mesh generation [49]. An adaptive solution is then completed by mapping the solution variables between the old and new spatial discretization. Apparently, the element or patch-based remap procedures [50] are not suitable for the meshfree method. In this study, we introduce a second-order accurate projection operation based on the meshfree convex approximation [29] to transfer the solution valuables. If the global remesher fails to decompose the computation domain into tetrahedrons, an alternative meshfree adaptive procedure is taken to reconstruct the neighboring nodal information using the old discretization. A sketch of the two-way adaptive procedure is illustrated in Figure 1.
The basic steps in an adaptive solution strategy are summarized as below.(a)Compute the initial solution.(b)Estimate the pointwise error indicator in solution and interactively trigger the mesh generator as error indicated.(c)Input the deformed mesh at the beginning of each adaptive step and use it to define the surface domain.(d)Approximate the surface domain by a generalized Farin algorithm [51].(e)Perform the anisotropic surface triangulation based on a metric map with Delaunay kernel [48] and yield the initial front of triangular faces. If there are any non-Delaunay facets in the constraining surfaces, insertion of Steiner points [52] will be made and surfaces’ edges and faces will be recovered by a series of swaps or flips. In this study we exclusively consider the nonintersecting facets as triangle.(f)Conduct the Advancing Front method which consists of filling the empty space defined by the initial front, creating tetrahedrons one at a time, and updating the front with created faces. A quadtree data structure [53] is utilized for the proximity search and insertion or deletion of the front points.(g)Reconstruct the neighboring nodal information using the old discretization if any step in (d)~(f) fails. Go to step (i).(h)Perform the projection operation, transfer the solution valuables, and update all nodal and internal valuables.(i)Initialize the solution and construct the approximations based on the new nodal distribution using the Lagrangian meshfree shape functions.(j)Resume the numerical solution procedure.
A customized pointwise error indicator introduced in step (b) is based on a shear deformation measure and is computed from the deformation gradient in (39). This error indicator is defined bywhere denotes the collection of integration points at time . The value 0.5 in (64) is made in average sense since the deformation gradient is generally not symmetric. When the pointwise error indicator exceeds an acceptable level, the adaptive procedure is triggered to create a new discretization in the deformed configuration. The implementation details of the anisotropic surface Delaunay triangulation and the Advancing Front tetrahedral mesh generation can be found elsewhere [48, 49] and therefore are omitted in this paper. After a new mesh is generated, the remap procedure is performed to transfer the solution valuables from the old spatial discretization to the new spatial discretization. We denote the variables before and after the each remap to be superscripted with “−” and “+,” respectively. Subsequently, we denote the unstructured tetrahedral meshes before and after remap at time = by and . For nodal value its remap is defined by the following projection operation:From the definition in (25), we haveor equivalentlywhere . Substituting (67) into (65) yieldswhere in (68) is the remap function which satisfies the following interpolation property:as well as the following linear polynomial reproduction property:
Since the remap algorithm preserves linearly varying nodal values, it is second-order accurate in space. In each remap procedure, , are the newly constructed meshfree shape functions evaluated on the current configuration based on the old spatial discretization . The tilde symbol in (65) stands for a “generalized” nodal quantity as defined in (25). If any boundary node does not reside in the old spatial discretization , an appropriate projection is performed to find the closet point on the boundary of for the subsequent remapping. We proceed to show that the above remap procedure is consistent in the sense that if the new discretization is identical to the old discretization, then all nodal quantities will remain unchanged after the remap:
Similarly, another consistent remap procedure is performed for the state valuable such as effective plastic strain, stress components, and other internal valuables sampled at the integration point per integration cell:where The matrices are computed on the set of integration points defined on the current configuration using the GMF method. The summation “” denotes the total number of integration points in the old discretization. Since the stress field is smoothed in the remap procedure, the pressure oscillation in the conventional displacement-based Element-Free Galerkin formulation can be effectively suppressed. On the other hand, it is known that the Advancing Front technique usually encounters difficulty when fronts merge. As a result, highly distorted tetrahedrons can occur in the generated mesh and greatly affects the accuracy of finite element solution. Since the construction of Lagrangian meshfree shape functions in this study does not rely on the finite element mesh, a distinct advantage of the proposed method in the adaptive procedure is its insensitivity to the existence of highly distorted tetrahedrons in the mesh.
A noteworthy addition to the proposed method is its flexibility in reconstructing the neighboring nodal information without remeshing. In conventional finite element analysis, when the global remesher is unable to generate the desired discretization at time under certain geometrical conditions, the adaptive procedure is aborted and causes termination in the simulation. With the current method, the second way for adaptive procedure will step in and replaces the global remeshing step. Under this scenario, a meshfree nodal reconstruction step is taken which maintains the material quantities for all nodal and integration points in the Lagrangian setting but requires a new neighbor search [45, 46]. Using the chain rule, the calculation for the deformation gradient becomeswhere is the decomposed deformation gradient due to the reconstruction of meshfree shape functions and is given byHere, we define to be a position vector on a new reference configuration at time . Since this meshfree nodal reconstruction step does not involve remeshing, the remap procedures are not needed. The updated deformation gradients together with the reconstructed shape functions and their derivatives are used to evaluate the discrete terms in Section 3.1. This flexibility in choosing different adaptive procedures provides an incentive for the use of adaptive Element-Free Galerkin method in large deformation analysis of manufacturing problems.
6. Numerical Simulation
In this section, we consider a FSW problem as shown in Figure 2. The light blue part is the tool, and the dark blue one is the work piece which is made of aluminum. The green part is rigid to provide the global constraints. The detailed dimensions are shown in Figure 3. The initial mesh of work piece is plotted in Figure 4, which is going to be refined with respect to the contact surface curvature of the tool. The tool has a conical shoulder surface to shape the material pushed out by two sequential stages of the FSW process:(1)plunge stage: the rotating tool plunges into the work piece along vertical direction;(2)traverse stage: the rotating tool travels along horizontal direction.
The material model of the work piece is assumed to be temperature dependant ideal plasticity. The material parameters are reported in Table 1. The temperature dependant yield stresses of the work piece are listed in Table 2. The melting point of work material is around 1080°C. The tool has a constant rotating speed 125 rad/s. The plunge displacement curve and traverse speed curve can be found in Figure 5. The thermal mortar contact is defined between the tool and work piece with the Coulomb frictional coefficient 0.7. A constant temperature 20°C is applied to all the parts as initial condition.
(a)
(b)
The normalized support size of the meshfree GMF approximation is 1.1. Since the density of nodal distribution varies dramatically throughout the domain due to local refinement, the actual nodal support size for every node is adjusted according to its surrounding nodal distribution to improve overall computational performance. By using local adaptivity, the integration cell size varies between 1 mm and 8 mm. The present numerical procedure was implemented into the LS-DYNA code [54] and the analysis was done by MPP double precision with 4 Xeon E5520 cores. There were around 400 adaptive steps that gradually increased the number of integration cells from ~4,000 to ~130,000.
Figure 6 plots the deformation in two stages, where the top view is to the left and the central cross-section view is to the right. The material flow and free surface are well captured and represented by the local adaptive remeshing. The red zone with large effective plastic strain clearly shows the stirring region during the welding process. Figure 7 shows the temperature results of both work piece and tool at three different processing steps. The model is plotted by cutting through the central cross section to better provide the temperature distribution inside the parts. At the end of plunge stage at s, the maximum temperature is around the front edge of the tool contact surface where the most heat sources from frictional contact are generated. In the traverse stage at s, there is a “V” shape contour due to high gradient temperature distribution in the stirring zone, and the maximum temperature is close to the melting point. Figure 8 gives the von Mises stress results, where material softening can be clearly observed as the temperature increases from the plunge stage to traverse stage. The energy curve of the work piece and the tool contact force curve are plotted in Figure 9.
(a) Plunge stage
(b) Traverse stage
(a) s
(b) s
(c) s
(a) s
(b) s
(c) s
(a)
(b)
7. Conclusions
The multiphysics thermomechanical complexity and severe material flow of friction stir welding process make analytical models incapable of capturing all the details needed for the satisfactory quantitative prediction of temperature and stress fields generated in the work piece. Via numerical modeling techniques, various finite element models of the friction stir welding process have been developed to help visualize and study the fundamental thermomechanical behavior of work piece and input parameters of the tool. However, there still exist great numerical challenges for the finite element simulation of friction stir welding process when the issues of handling severe mesh distortion, modeling contact and traction-free boundary conditions, performing accurate state variable remap, and maintaining the quality of adaptive mesh are simultaneously presented in the model. This paper attempts to provide an alternative approach using the adaptive Element-Free Galerkin method to overcome those numerical challenges. From the best of authors’ knowledge, this is the first paper that describes a three-dimensional adaptive Element-Free Galerkin method for the thermomechanical analysis.
In our approach, the coupled thermomechanical equations using a staggered explicit time integration scheme are modeled within a Lagrangian framework. This approach facilitates a direct incorporation of the Element-Free Galerkin formulation with the adaptive procedures to circumvent the severe mesh distortion problems. The concept of the two-way adaptive procedure is introduced to bypass the numerical difficulty caused by the abortion of adaptive mesh generation. The Lagrangian meshfree convex approximation plays a key role in simplifying the boundary condition enforcement, suppressing tensile instability, minimizing the adaptivity-induced discretization sensitivity, and offering an accurate and consistent projection operation in the remap procedure. In comparison to the existing numerical methods, the present method shares the advantages of Lagrangian approach in solid mechanics applications as well as resolves the mesh distortion problems in extremely severe deformation analyses. Subsequently, the residual stress analysis based on the present numerical procedure demands no special numerical treatments and could be easily solved by a spring-back solution of the implicit formulation which will be addressed in the near future.
The proposed method meets all the computational requirements suggested by D. M. Neto and P. Neto [3] for FSW analysis including (1) consideration of rotational boundary condition; (2) consideration frictional contact; (3) support for very high levels of deformation; (4) support for elastic-plastic or elastic-viscoplastic material models; and (5) support for complex geometry. A further comparison with the experimental data requires a remodeling of the problem with a more realistic FSW environment such as the tool profile, feeding angle, friction coefficients, and the temperature dependent material constants, and the results will be presented separately in another paper.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors wish to thank Dr. John O. Hallquist of LSTC for his support to this research. The authors would also like to thank Dr. Art Shapiro for the technical discussions.