EULERIAN FORMULATION WITH VOLUME OF SOLID (VOS) APPLICATION TO METAL FORMING by Khaled S. Al-Athel B.Sc., M.Sc. King Fahad University of Petroleum and Minerals (KFUPM), 2005 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in The Faculty of Graduate Studies (Mechanical Engineering) The University of British Columbia (Vancouver) September 2010 © Khaled S. Al-Athel 2010 ii ABSTRACT The volume of fluid (VOF) method has been used extensively in computational fluid dynamics (CFD) applications. The method is used to track the unconstrained flow of the fluid at the free surface in an Eulerian formulations framework. In this work we propose the use of the VOF method in nonlinear solid mechanics and metal forming applications, therefore the new method is called the Volume of Solid (VOS) method. The VOS method is to be used with Eulerian finite element (FE) formulations. In Eulerian formulation, the mesh is fixed in space and the material flows freely within the domain. This type of formulation overcomes many difficulties that are associated with other FE formulations such as mesh distortion and the need for re-meshing or a proper mesh motion scheme. The main drawback of the Eulerian formulation is the difficulty in tracking the unconstrained flow of the material at the free boundaries; therefore an additional scheme is needed to overcome this problem. The derivation of the VOS equations and its implementation in an Eulerian formulation along with applications in metal forming and metal cutting are the main goals of this research. In this work, we present the derivation of the Eulerian quasi-static and dynamic equilibrium equations. The complete development of the VOS method is presented for the commonly used uniform mesh using a piecewise linear interface construction method, and for a general case of a non-uniform mesh using a new filling pattern scheme. The formulation of the VOS method includes the calculation of the fractional volume value of the solid in each element and finding a directional vector to locate the interface of the free surface. The implementation of the VOS method is discussed by covering important topics that include the connectivity of the free surface interfaces in addition to introducing and modeling a new algorithm for the rigid-flexible contact problem. The formulation presented in this work is implemented in a VOS FE-code and used to simulate several metal forming and metal cutting applications. The VOS simulations are validated by comparing the results with other experimental and numerical work presented in the literature. iii TABLE OF CONTENTS ABSTRACT ........................................................................................................................ ii TABLE OF CONTENTS ................................................................................................... iii LIST OF TABLES ............................................................................................................. vi LIST OF FIGURES .......................................................................................................... vii NOMENCLATURE .......................................................................................................... xi ACKNOWLEDGEMENTS ............................................................................................. xvi DEDICATION ................................................................................................................ xvii CHAPTER 1 INTRODUCTION .................................................................................... 1 1.1 Motivation ............................................................................................................ 1 1.2 Metal Forming Processes ..................................................................................... 2 1.3 Numerical Methods in Metal Forming ................................................................. 3 1.3.1 Traditional Large Strain Finite Element Formulations ................................. 3 1.3.2 Arbitrary Lagrangian-Eulerian (ALE) Formulation ..................................... 6 1.3.3 Meshless Methods ......................................................................................... 7 1.4 Scope of Research ................................................................................................ 8 CHAPTER 2 SURFACE TRACKING METHODS .................................................... 11 2.1 Traditional Free Boundary Methods .................................................................. 11 2.2 Volume of Fluid (VOF) Method ........................................................................ 13 2.2.1 The VOF Method ........................................................................................ 13 2.2.2 VOF in the Literature .................................................................................. 14 2.3 Transition from VOF to VOS............................................................................. 21 CHAPTER 3 EULERIAN FORMULATION .............................................................. 23 3.1 Quasi-Static Analysis ......................................................................................... 24 3.1.1 Incremental Decomposition ........................................................................ 25 iv 3.1.2 Linearization ............................................................................................... 27 3.1.3 Treatment of Convective Terms ................................................................. 28 3.1.4 Eulerian Equilibrium Equation for Quasi-static Analysis .......................... 30 3.2 Dynamic Analysis .............................................................................................. 31 3.2.1 Virtual Work Done by Inertia Forces ......................................................... 32 3.2.2 Decomposition of Velocity and Acceleration ............................................. 32 3.2.3 Linearization of the Referential Inertia Term ............................................. 33 3.2.4 Linearization of the Convective Inertia Term ............................................. 34 3.2.5 Eulerian Equilibrium Equation for Dynamic Analysis ............................... 36 3.3 Discretization of the Eulerian FE Equations ...................................................... 37 3.3.1 Discretization of the Quasi-Static Eulerian Equation ................................. 37 3.3.2 Discretization of the Dynamic Eulerian Equation ...................................... 40 CHAPTER 4 VOLUME OF SOLID (VOS) METHOD............................................... 44 4.1 VOS for a Uniform Mesh ................................................................................... 46 4.1.1 Solution of the Transport Equation ............................................................. 46 4.1.2 Free Surface Directional Vector ................................................................. 47 4.1.3 Construction of the Free Surface using a Geometrical PLIC Scheme ........ 52 4.2 VOS for a Non-Uniform Mesh........................................................................... 56 4.2.1 Modification to the Solution of the Transport Equation ............................. 56 4.2.2 Generalization of Free Surface Directional Vector .................................... 57 4.2.3 Construction of the Free Surface using a Filling Pattern Scheme .............. 59 CHAPTER 5 SOLUTION CONTROLS AND IMPLEMENTATION ........................ 62 5.1 Time Step Size and Stability .............................................................................. 62 5.2 Updating the Velocity Field ............................................................................... 63 5.3 Tracking Material Point Properties .................................................................... 64 v 5.4 Free Surface Connectivity .................................................................................. 65 5.5 Contact Analysis ................................................................................................ 69 5.6 Heat Transfer Analysis ....................................................................................... 72 5.7 Solution Procedure ............................................................................................. 73 CHAPTER 6 METAL FORMING APPLICATIONS .................................................. 76 6.1 Punch Indentation Problem ................................................................................ 76 6.2 Compression between Wedge-Shaped Dies ....................................................... 81 6.3 Sheet Metal Extrusion ........................................................................................ 86 6.4 Backward Extrusion ........................................................................................... 91 6.5 Local Loading of a Simply-Supported Plate ...................................................... 98 6.6 V-Bending ........................................................................................................ 104 6.7 Three-Point Bending ........................................................................................ 111 CHAPTER 7 METAL CUTTING APPLICATION ................................................... 117 7.1 Experimental Work and Numerical Setup ....................................................... 121 7.2 VOS Simulation of Metal Cutting .................................................................... 122 7.3 Results and Discussion ..................................................................................... 124 CHAPTER 8 CONCLUTIONS AND FUTURE WORK ........................................... 132 8.1 Summary and Conclusions ............................................................................... 132 8.2 Future Work ..................................................................................................... 134 REFERENCES ............................................................................................................... 137 vi LIST OF TABLES Table 4.1 Summary of the free surface interface parameters calculations ................. 55 Table 6.1 Contact pressure and strain comparison between VOS and [80] ................ 97 Table 6.2 von-Mises stress and strain comparison between VOS and ANSYS for the V-bending problem ................................................................................... 110 Table 7.1 List of experimental and numerical metal cutting setups ......................... 121 Table 7.2 Material properties for the workpiece and tool used in metal cutting ...... 122 vii LIST OF FIGURES Figure 1.1 Initial and final configuration of a punch indentation problem using (a) Lagrangian formulation and (b) Eulerian Formulation ............................ 4 Figure 2.1 Distributionof„f‟atfreesurfaceandneighboringelements.................. 13 Figure 2.2 Free surface averaging using (a) SLIC and (b) PLIC ............................. 15 Figure 2.3 Schematic of regular and hanging nodes in adaptive grid refinement ... 18 Figure 4.1 Definition of free surface element parameters ....................................... 47 Figure 4.2 Definition of the free surface directional vector ..................................... 51 Figure 4.3 Construction of fictitious elements at the right end of the domain ........ 51 Figure 4.4 Four quadrant categorization of the free surface .................................... 52 Figure 4.5 Free surface shape possibilities in a uniform mesh ................................ 54 Figure 4.6 3x3 non-uniform stencil with outward facial normals ........................... 57 Figure 4.7 3x3 non-uniform stencil with unit vectors .............................................. 59 Figure 4.8 Free surface sectioning in the filling pattern scheme ............................. 60 Figure 5.1 Velocity field update for newly added nodes ......................................... 63 Figure 5.2 Tracking material points at the boundary using the free surface directional vectors .................................................................................. 65 Figure 5.3 Default free surface connectivity scheme ............................................... 66 Figure 5.4 Contact-modified free surface connectivity ........................................... 68 Figure 5.5 Effect of weighting factor on the free surface connectivity process ...... 68 Figure 5.6 Penetration detection and surface repositioning in contact analysis ...... 71 Figure 5.7 Flowchart of the VOS FE-code solution procedure ............................... 75 viii Figure 6.1 Geometrical setup for the punch indentation problem ........................... 78 Figure 6.2 Evolution of the workpiece deformed shape in the punch indentation problemusing VOS ................................................................................. 78 Figure 6.3 von-Mises plastic strain distribution for the punch indentation problem 79 Figure 6.4 von-Mises stress distribution for the punch indentation problem .......... 80 Figure 6.5 Setup for the compression between wedge-shaped dies problem .......... 82 Figure 6.6 Evolution of the workpiece deformed shape in the compression between wedge-shaped dies problem using VOS ................................................. 83 Figure 6.7 von-Mises strain distribution in the compression between wedge-shaped dies problem ........................................................................................... 84 Figure 6.8 Load-strain comparison curves for the compression between wedge- shaped dies problem ............................................................................... 85 Figure 6.9 Geometrical setup for the sheet metal extrusion problem ...................... 88 Figure 6.10 Evolution of the billet deformed shape in the sheet metal extrusion problem using VOS ................................................................................ 89 Figure 6.11 von-Mises strain distribution for the sheet metal extrusion problem ..... 89 Figure 6.12 Longitudinal stresses at different lateral positions for the sheet metal extrusion problem ................................................................................... 90 Figure 6.13 Longitudinal stresses along the thickness of the billet in the sheet metal extrusion problem ................................................................................... 90 Figure 6.14 Geometrical setup for the backward extrusion problem ......................... 93 Figure 6.15 Evolution of the workpiece deformed shape in the backward extrusion problem using VOS ................................................................................ 94 ix Figure 6.16 Load-Stroke curve for the backward extrusion problem ........................ 95 Figure 6.17 Contact pressure (MPa) between workpiece and die at different strokes in the backward extrusion problem ........................................................ 95 Figure 6.18 von-Mises strain distribution for the backward extrusion problem ........ 96 Figure 6.19 Geometrical setup for the plate indentation problem ........................... 100 Figure 6.20 Deformed shape of the plate indentation problem using VOS ............. 101 Figure 6.21 Load-Displacement curve for the flat plate indentation problem ......... 101 Figure 6.22 Load-Displacement curve for the hemispherical plate indentation problem ................................................................................................. 102 Figure 6.23 Comparison of the load-displacement curves of the flat and the hemispherical indenters ........................................................................ 102 Figure 6.24 von-Mises strain distribution using the flat-face indenter in the plate indentation problem .............................................................................. 103 Figure 6.25 von-Mises strain distribution using the hemispherical indenter in the plate indentation problem ..................................................................... 103 Figure 6.26 Geometrical setup for the V-bending problem ..................................... 107 Figure 6.27 Evolution of the deformed shape from Experiment and VOS in the V- bending problem ................................................................................... 108 Figure 6.28 Load-Stroke curve for the V-bending problem .................................... 108 Figure 6.29 von-Mises stress distribution for the V-bending problem .................... 109 Figure 6.30 von-Mises strain distribution for the V-bending problem .................... 110 Figure 6.31 Geometrical setup for the three-point bending problem ....................... 113 x Figure 6.32 Evolution of the three-point bending deformed shape from the Experiment and VOS ............................................................................ 114 Figure 6.33 Load-Stroke curve for the three-point bending problem ...................... 115 Figure 6.34 von-Mises stress distribution for the three-point bending problem ...... 115 Figure 6.35 von-Mises strain distribution for the three-point bending problem ...... 116 Figure 7.1 Schematic view of metal cutting process with chamfered tool ............ 120 Figure 7.2 Evolution of chip deformation in the VOS metal cutting simulation ... 124 Figure 7.3 Velocity fields for sharp tool, -10o chamfered tool, -25o chamfered tool and a blunt tool ..................................................................................... 127 Figure 7.4 Effective stress distribution for sharp tool, -10o chamfered tool, -25o chamfered tool and a blunt tool ............................................................ 128 Figure 7.5 Plastic strain rate distribution for sharp tool, -10 o chamfered tool, -25 o chamfered tool and a blunt tool ............................................................ 129 Figure 7.6 Comparison of maximum normal and tangential stresses between VOS and ALE simulations ............................................................................ 130 Figure 7.7 Comparison between the cutting forces obtained using VOS, ALE and experimental results .............................................................................. 130 Figure 7.8 Comparison between the thrust forces obtained using VOS, ALE and experimental results .............................................................................. 131 xi NOMENCLATURE t Quantity occurring at time t t Quantity referred to time t a Element incremental material acceleration vector t ia Components of referential material acceleration vector ia Components of incremental referential material acceleration vector 1BAt Element shape function derivative matrix related to first convective stiffness 2BAt Element shape function derivative matrix related to second convective stiffness 1BEt Element shape function derivative matrix related to Eulerian material stiffness 2BEt Element shape function derivative matrix related to Eulerian geometrical stiffness Ct Element elastic-plastic material constitutive matrix 1CAt Element first convective velocity-stiffness matrix 2CAt Element second convective velocity-stiffness matrix 3CAt Element third convective velocity-stiffness matrix 4CAt Element fourth convective velocity-stiffness matrix t ijklC Components of fourth order material constitutive tensor c Specific heat t klD Components of rate of deformation tensor xii D Effective strain rate pD Effective plastic strain rate tdS Elemental surface area t CdS Contact surface area ds Thickness of small area in the filling pattern scheme tdV Elemental volume E Young‟smodulus TE Tangent modulus t e Element strain vector t ije Components of strain tensor ft Element internal force vector extft Element internal force vector f Fractional volume value crf Critical fractional volume value t f Arbitrary function at time t t f Material time derivative of arbitrary function t f Grid derivative of arbitrary function t B if Components of body force per unit mass t S if Components of surface traction vector xiii t C if Components of contact traction vector H Element interpolation matrix kh Element shape function at node k h Convection heat transfer or film coefficient Kt A Element convective stiffness matrix 1Kt A Element first convective stiffness matrix 2Kt A Element second convective stiffness matrix Kt E Element Eulerian stiffness matrix 1Kt E Element Eulerian material stiffness matrix 2Kt E Element Eulerian geometrical stiffness matrix k Conductivity coefficient Mt A Element convective mass matrix Mt E Element Eulerian mass matrix N Number of nodal points IN Nodal shape function at free surface node I n Unit outward normal vector t in Components of unit normal vector q Rate of heat generation per unit volume ri Free surface directional vector 1St A Element stress matrix related to first convective matrix 2St E Element stress matrix related to second convective matrix S Total thickness of free surface area in the filling pattern scheme xiv yS Yield stress T Temperature t Time u Element incremental material displacement vector u g Element incremental grid displacement vector iu Components of incremental material displacement vector g iu Components of incremental grid displacement vector V Element total volume v Element incremental material velocity vector vg Element incremental grid velocity vector iv Components of incremental material velocity vector g iv Components of incremental grid velocity vector t iv Components of material point velocity vector t g iv Components of grid point velocity vector t iv Components of material acceleration vector t xkv Components of nodal incremental material velocity vector vi Average velocity pv Penetration vector Kv Corrective velocity constraint at free surface point K t extW External work w Width of small area in the filling pattern scheme xv m iX Components of material point initial coordinates g iX Components of grid point initial coordinates t ix Components of material point position vector t ix Components of nodal coordinate vector s Element edge length t Time increment Virtual or variation Equivalent strain Poisson‟sratio Free surface angle t Material density σt Element stress vector t ij Components of Cauchy stress tensor t T ij Components of Truesdell stress rate tensor Effective stress xvi ACKNOWLEDGEMENTS I would like to express my sincere gratitude and appreciation for my advisor Dr. Mohamed S. Gadala for his continuous support, encouragement and guidance throughout the stages of this research and my course work. I would be in debt for his invaluable support for all of my future academic life. I would also like to thank the members of my examining committee, Dr. Yusuf Altintas, Dr. Farrokh Sassani and Dr. Akram Alfantazi for their helpful discussion and their valuable time in reviewing this work A special thank for Dr. A. M. Arif for his support and his knowledge which sparked my interest in the area of Finite Element Method (FEM). I also express my gratitude for Dr. M. Sunar who served as my advisor during my Master‟sdegreewho influenced me greatly during the early stages of my academic career. I am also grateful for the support of Dr. Y. Al-Nassar and for his encouragement. xvii DEDICATION To my mother for her continuous love and support throughout my whole life. To my father who never stopped believing in me and for his guidance and help during all of my academic life. I will always be in debt to them and will try to make them proud for the rest of my life… 1 CHAPTER 1 INTRODUCTION 1.1 Motivation Manufacturing, in its broadest sense, is the process of converting raw materials into products. It encompasses the design and fabrication of goods by means of various production methods and techniques. Manufacturing processes includes a wide variety of metal forming and machining applications, each with a different problem definition and running conditions. Before numerical methods became an important tool in solving problems in the field of engineering mechanics, attaining an exact solution for metal forming applications was quite difficult to accomplish due to the complexity of these applications. Pairing theoretical solutions with empirical relations did help in the progress of this field but the problem was that these solutions were limited to basic cases with simple geometries and material models. Even if some of those solutions were tabulated for a range of different material properties and geometrical variables, they were not very practical, especially in prototype design stages. Moreover, empirical relations are usually derived based on experiments performed under certain conditions and it is not clear how changing some of the parameters would impact the accuracy of these relations. Therefore it was clear that additional tools are needed in order to obtain broad and reliable results for metal forming applications. The vast evolution in numerical techniques and Finite Element Methods (FEM) in theearly1980‟sresultedinaninnovativeadvancementin the metal forming industries. The ability to simulate different metal forming processes while covering different operating conditions, geometrical setups and material models reduced the time and cost of having to perform a different experiment for each set. Once the solution from a numerical simulation is verified by experimental or available theoretical/empirical results, it can be used in the prototype design by studying the effect of the process 2 conditions on the final product in addition to predicting potential forming defects. Numerical simulations can also be used to achieve a better and more comprehensive understanding of the process mechanical properties, effect of the forming loads on the final geometry of the workpiece and the development of residual stresses in addition to different damage and failure concerns depending on the process. Even with the ongoing significant improvement innumericaltechniquesandFEM‟s, metal forming simulations still face some problems in terms of accuracy and convergence problems for certain applications. Additionally, it is still very challenging to develop a general FEM approach that can be used to simulate different processes with the same level of reliability without having to simplify the problem on the expense of the accuracy. 1.2 Metal Forming Processes Three classes of manufacturing processes are of interest in this work; bulk deformation processes, sheet metal forming and metal cutting. In bulk deformation applications (extrusion, compression between dies, coining…etc), the shape of the workpiece is changed by plastic deformation under loads enforced by tools. Usually the workpiece is confined within a set of dies which controls the final shape of the workpiece. Sheet metal forming (plate bending, deep drawing…etc) applications are usually performed for workpieces with high ratio of area to thickness. The tool, in these types of applications, is what controls the final shape of the workpiece by pressing against it forcing it to take a certain shape. Metal cutting is a process that involves removing material from the workpiece. In some cases, the final manufactured products often require additional processing and machining to improve certain aspects such as the surface finish or dimensional accuracy. All metal forming processes are considered to be complex problems because they involve both geometrical and material nonlinearities. Large deformations are usually considered to be one of the main sources of difficulty when it comes to simulating metal forming. Large boundary motion is also expected in such processes and requires accurate treatment to maintain the shape integrity of the workpiece. Modeling the interaction 3 between the workpiece and the tools/dies is very important to produce realistic operating conditions. The nonlinear material behavior is also very important in simulating metal forming processes and sometimes affects the interaction between the workpiece and the tool depending on the values of the material properties. Therefore, finite element simulations require large strain formulations that can handle both material and geometrical nonlinearities as well as nonlinear boundary conditions. 1.3 Numerical Methods in Metal Forming 1.3.1 Traditional Large Strain Finite Element Formulations Two numerical formulations are commonly used in finite element (FE) simulation of large deformation solid mechanics problems: the Lagrangian and the Eulerian formulations. Each has its own advantages and disadvantages. One of the earliest large strain finite element formulations can be traced back to Hibbitt et al. [1] which was referred to as Total Lagrangian formulation. The method was later developed into a modified version named Updated Lagrangian formulation by McMeeking and Rice [2]. The basic difference between the two formulations is in the reference configuration. In Total Lagrangian, the initial configuration of the problem is used in the linearization of the incremental equation of motion. On the other hand, Updated Lagrangian uses the last calculated configuration as a reference for the calculation at the next incremental step. Generally, both formulations should provide the same results. The choice between the two formulations mostly depends on their relative effectiveness and numerical efficiency [3]. In Lagrangian formulation, the finite element mesh is attached to the material points and possesses the same velocity of the material, therefore undergoing the same deformation. Tracking the free surface is easily achieved since the material is attached to the mesh point but this leads to mesh distortion and convergence problems in applications that undergo large deformations of the material. Figure 1.1 shows a simple punch indentation problem where a rigid tool presses against the workpiece. The deformation of the material using Lagrangian formulation can be seen in Figure 1.1a. 4 To remedy the mesh distortion problem in Lagrangian formulation, researchers and commercial software programs have relied on „Re-Meshing‟ procedures [4]. In this technique, the whole domain (or part of it) of the problem is re-meshed after a specified number of solution steps and/or after some mesh distortion criteria has been reached. Then, the program transfers problem parameters such as stress and strain to the new mesh through interpolation schemes. This re-meshing approach has other major drawbacks. One problem is that it is time consuming because of the data transfer process and the calculation of all element and global equations for the new mesh. Another problem is the accuracy degradation due to the transfer of the material point properties and solution variables unless the interpolation scheme is chosen properly to ensure accurate and consistent data transfer from the old mesh to the new one. Another problem in using Lagrangian formulation is the difficulty in implementing nonlinear boundary conditions that change during deformation. Punch Initial configuration (a) Lagrangian mesh (b) Eulerian mesh Figure 1.1 Initial and final configuration of a punch indentation problem using (a) Lagrangian formulation and (b) Eulerian Formulation 5 The Eulerian formulation is suitable for problems with fluid-like material behavior, as in metal forming applications. The adaptation of this formulation, which is commonly used in the area of fluid mechanics, to large strain and metal forming problems has been used and presented in the literature [5-8]. Figure 1.1b shows the final configuration of the punch indentation problem using Eulerian formulation. The finite element mesh is spatially fixed and the material is allowed to flow freely within the fixed domain, hence no mesh distortion occurs. However, the Eulerian approach suffers from the difficulty in handling deformation dependent boundary conditions and in tracing history dependent properties as in the case of plasticity. Another drawback of this formulation is its incapability in tracking the free boundaries and modeling the unconstrained flow of the material within the fixed mesh causing elements to be partially occupied by the material. This drawback makes this approach limited to the simulation of steady-state processes, excluding the initiation of the process. Different attempts were made to overcome this problem of the free boundaries by coupling the Eulerian formulation with a scheme that can handle the unconstrained flow of the material. One of the most popular methods in the literature is the pseudo- concentrations technique [9-13]. The pseudo-concentrations technique is based on assigning a pseudo-concentration (ci for i = 1,…, N) for each node that defines the presence of the material in any node (ci > 0 indicates the presence of the material at node i, ci < 0 indicates the absence of the material at node i and ci = 0 means that node i is on the free surface). The pseudo-concentrations are then interpolated between the nodes within each element to estimate the free surface which requires a smoothing algorithm. One of the drawbacks of this method is that the formulation depends on the number of free surfaces, which should be specified before deriving the finite element equations of this method. This means that a different formulation is required for each application with a different number of free boundaries. The method also cannot handle the contact between the material and another body within the Eulerian domain. The contact edge must be modeled prior to the simulation with the second contact body being modeled outside of the Eulerian domain as a rigid body or a Lagrangian mesh. Another drawback is the fact that each node with no present material is actually given an artificial low viscosity rather than excluding it from the analysis. In addition to the nodal interpolations 6 within each element, this may lead to numerical inaccuracies and convergence problems depending on the size of the problem domain. Other methods include an operator split method [14] and particle level set [15]. Both methods are limited to cases with only uniformly constructed mesh. In addition to that, both suffer from the difficulty in handling the contact within the Eulerian domain. 1.3.2 Arbitrary Lagrangian-Eulerian (ALE) Formulation It is generally conceived that neither the Lagrangian nor the Eulerian approaches is well-suited for modeling deformation processes involving both large deformations and unconstrained flow of material with free boundaries. Each formulation excels at handling the weak points of the other method since the Lagrangian approach can handle the unconstrained flow of the material by definition, whereas the Eulerian approach does not suffer from the problem of mesh distortion. The Arbitrary Lagrangian-Eulerian (ALE) formulation emerged as a new formulation that combines the merits of both the Lagrangian and Eulerian ones [16, 17]. In ALE formulation, the finite element mesh is neither attached to the material (i.e. Lagrangian) nor fixed in space (i.e. Eulerian) but rather assigned an arbitrary motion independent of the material deformation [18, 19]. As the material deforms, the elements should move in a way that would optimize their shape to minimize mesh distortion while following the material flow on the boundaries to satisfy the boundary constraints that prevents the relative motion between material and mesh points normal to the boundary. This ensures that the mesh would have the same boundary as the material throughout the whole analysis. The implementation of the ALE formulation in metal forming processes has become a major tool in the simulation of many applications. Different rolling and extrusion cases were tested using ALE and compared with results from the literature to validate the efficiency of the method once a proper mesh motion scheme is implemented [20-22]. Also ALE was used in simulating more complex and advanced applications such as the metal cutting process [23-25]. Some attempts have been made to employ ALE 7 formulation in 3D applications [26, 27], although these attempts utilize an operator split approach where the mesh motion is uncoupled into an Updated Lagrangian step and an Eulerian step rather than“true”coupledALEapproach[19-22]. In contrast to the Lagrangian re-meshing technique, in ALE, the material point variables are transferred from the grid points to the material points using proper continuum mechanics equations and therefore the accuracy of the results is not compromised by this transformation. The ALE formulation entails, however, some implementation difficulties and drawbacks. One of these difficulties is to devise a proper mesh motion scheme for each application since a generalized procedure for a class of applications is difficult to attain. This makes the efficiency of the ALE approach dependable of the mesh motion scheme advised by the user and the type of application used in the simulation. Other problems include the difficulty in modeling contact in metal forming simulations, the inherent convergence problems that comes from the complexity of implementing consistent and coupled ALE formulations especially in dynamic applications, and the fact that the number of unknowns is increased, enclosing both the material and mesh displacements [28]. 1.3.3 Meshless Methods A new family of numerical computational methods, called meshless methods [29], emerged in the last two decades to mainly overcome mesh related problems in simulating metal forming processes. These methods are based on a constructed domain of nodes that represents the material without the need for any elements. The analysis is carried by solving partial differential equations on the domain of the irregular distribution of nodes. Although these methods eliminate any difficulties in dealing with the mesh that arise from using finite element methods, they do suffer from several drawbacks. The main issue with this approach is that the number and initial distribution of nodes have a great influence on the final shape of the workpiece. Following the unconstrained flow of the material without any elements is very difficult to accomplish, especially in problems with very large plastic deformations like the process of backward extrusion 8 [30]. Guan et al. [31] studied the effect of increasing the number of nodes on the final deformed shape of the workpiece. A major drawback observed from this method is the material volume loss during the simulation which means that a large number of nodes are required to minimize the volume loss. Since the method does not have any elements, a special interpolation process must be applied to present a smooth distribution of the data in addition to imposing essential boundary conditions to keep the outer boundary intact during the analysis [32]. Nodes must be laid out properly at the beginning of the simulation when dealing with curved tools and dies. Since there are no elements, there must be a node at each corner and center of a curve to have an accurate presentation of the shape [33, 34]. Different attempts were made to develop a coupled FE-meshless approach [35, 36]. The coupled approach starts with a regular finite element setup and when part of the mesh reaches a critical stage or gets distorted (based on a condition defined by the user) it gets converted into a meshless region and remains that way for the duration of the simulation. Although this approach is more efficient than complete meshless methods, it still shares many of the difficulties and drawbacks. In summary, meshless methods are usually limited to either steady-state metal forming processes or very confined bulk deformations (e.g. closed die forging). That is due to the fact that the contact boundary conditions must be defined and applied prior to the simulation otherwise the method will not be able to properly follow the unconstrained flow of the metal. Additionally, the number of nodes not only affects the accuracy of the results but also the final shape and the volume loss of the material. 1.4 Scope of Research Numerical simulation is a powerful and effective tool in simulating metal forming applications, particularly complex systems. The use of such tools can help in predicting design flaws and identifying problems associated with each process before reaching the prototype testing stage. A wide range of defects can be detected using numerical simulations. In bulk deformation processes, such defects include surface cracking, internal cracking, buckling and alligatoring among others. In sheet metal forming, factors 9 that can affect the process include the elongation of the yield point, development of residual stresses, springback effect and wrinkling. In metal cutting, the chip curl is a very important subject and requires many test/simulations to study. Excessive curling may cause the chip to fold into the uncut part of the workpiece or into the cutting tool which will cause problems due to the entanglement. Studying the temperature distribution is also very important for failure and fatigue analysis of the final products. All of the mentioned defects and factors can be studied and observed with the help of finite element analysis and numerical simulations It is clear that traditional finite element formulations proved to be incapable of simulating metal forming processes. That led to the derivation of the ALE formulation combining the strength of both approaches while eliminating most of their drawbacks. On the other hand, combining two different formulations did raise several difficulty issues when using ALE. Devising a proper and robust mesh motion scheme for each application is not an easy task to attain, especially for non-steady-state processes. One of the very well established free surface tracking schemes in the literature is the Volume of Fluid (VOF) method [37]. The method was presented to be used with Eulerian formulation to capture the free surface of the fluid by calculating the fractional volume value of the fluid at each cell, then constructing the free surface based on these values. In this work we are proposing the use of the Eulerian formulation in simulating metal forming processes with the implementation of a tracking scheme that is based on the VOF method. The resulting method, to be called “VolumeofSolid(VOS)method”, is to be used in simulating large deformation and metal forming applications. As mentioned above, the use of Eulerian approach has the advantage of not dealing with mesh distortion and eliminating the need of modeling a proper mesh motion scheme. The implementation of the VOS method will overcome the free boundary problem by properly tracking the location of the free surface of the solid. A critical review of the VOF method is presented in Chapter 2 highlighting the important development of the method, since it was first presented by Hirt and Nichols [29], and the different approaches in constructing the free surface. This gives a better 10 understanding of the method which lays the foundation behind transforming the method from fluid mechanics to solid mechanics. Based on the principle of virtual work and continuum mechanics, a complete derivation of the Eulerian formulation is presented in Chapter 3 for quasi-static and dynamic analyses. Treatment of the convective terms is also included in the formulation. A section on heat transfer analysis is also included in this chapter. The proposed Volume of Solid (VOS) method is discussed in details in Chapter 4. The discussion follows three major steps in the formulation of the VOS method; calculating the fractional volume value of the solid, defining the location and orientation of the free surface, and constructing the free surface. The method is first presented for uniformly constructed Eulerian grids. After that, a new scheme is presented to generalize the method for applications with non-uniform Eulerian grids. Other important features of the method are also discussed including tracking material point related properties and updating newly added nodes in the analysis. Chapter 5 covers topics related to the implementation of the VOS method. A simple in nature yet effective averaging process is incorporated in the method to avoid any smearing of the data and to eliminate the need for any smoothing function. A contact algorithm is presented for rigid-flexible contact cases. The algorithm allows the modeling of a moving rigid body within the Eulerian domain and tracking the change of the contact boundary conditions at the free surface of the solid during the course of the simulation. A summary of the method and its implementation in a 2-D computer code is presented at the end of the chapter. In Chapters 6 and 7, various metal forming and machining applications are simulated using the proposed VOS method. The results are compared to experimental measurements and/or numerical results from the literature. The cases simulated in these Chapters covers a wide range of processes to highlight different aspects of the VOS method. Finally, Chapter 8 concludes this work with a summary of the method and a discussion on its implementation in solid mechanics and metal forming problems. 11 CHAPTER 2 SURFACE TRACKING METHODS The use of Lagrangian or ALE formulations in computational fluid dynamics (CFD) problems was a common choice of formulation due to the simple fact that the fluid surface is confined within the elements at all times. This led to relatively straightforward computations. However, these formulations still have difficulties dealing with mesh distortion and complexity in the shape of the outer surface in some problems due to the folding and merging of the surface. With the difficulties in using Lagrangian and ALE, Eulerian formulation became the main tool in modeling CFD problems. Even though the fluid flows freely within a fixed mesh, it is customary to identify the fluid in an Eulerian mesh as elements on which computations are performed. The fluid moves between elements based on the velocity calculations performed at each step. Therefore, the fluid flow requires averaging of the flow properties for the fluid which may lead to some smearing and discontinuity in the free surface, and that is where the need for a special treatment of the free boundaries arises. Free surface tracking schemes are necessary to have an acceptable simulation of the fluid flow and representation of the free surface. The goals of these schemes are to be able to identify Eulerian cells that contain a free surface, describe the locations and shape of the free boundaries, enforce desired surface boundary conditions on the computational mesh and the newly constructed surfaces, and finally to compute the advancement of the free surface after each step. 2.1 Traditional Free Boundary Methods When using an Eulerian mesh, the problem of the free surface is still a concern, but this type of formulation allows different free boundary capturing/tracking schemes to be incorporated into the main formulation to capture the free surface shape. At the same 12 time, this formulation can model the complex shape of the fluid free surface. Using Eulerian formulation, free boundary tracking schemes may be categorized into two groups: surface and volume tracking methods. One of the earliest and simplest surface tracking methods is the Height Functions [38]. In this method, the free boundary is defined by its distance from a reference line as a position function along that reference line. This means that the method can only be used for a vertical or horizontal direction at a time. As a standalone scheme, Height Functions are very limited to simple shapes and small slope changes. A more generalized form of the Height Functions is the method of Line Segments [39]. This approach uses a series of short lines, or point, that are connected by line segments. Compared with Height Functions, Line Segments present a better resolution of the free surface but on the expense of the computational cost. Another surface tracking method is the Marker Particles [40], which, instead of defining the free surface directly, defines particles at the fluid-occupied regions in the surface with each particle specified to move with the fluid velocity at its location. A free surface is identified when a mesh element contains markers but neighbors an element with no markers. Additional computation must be performed for the location and intersection of the free surface with the element. The main problem with surface tracking methods occurs when there is an intersection between two or more surfaces. Additionally, these methods are not able to handle cases where the fluid exhibits severe surface folding and merging over on itself without undergoing rigorous mesh distortion. Another problem with surface tracking methods is that they require large computer storage to store all the point coordinates. Volume tracking schemes are based on defining a tracer or a marker particle at each fluid region of the domain and not only at the surface. One of the most successful methods in this category was the Marker-and-Cell (MAC) method [41]. It was the first method to use pressure and velocity as the primary variables. The method is based on employing a distribution of marker particles all over the fluid region, and set free surface pressures at the centers of the surface cells. The method was improved with time to 13 include the pressure boundary condition at the location of the boundary within the cell that defines the free surface. The MAC method still requires large computer storage, plus it experiences a problem with non-natural creation of high and low marker number densities due to the irregularity of the flow field in some cases. 2.2 Volume of Fluid (VOF) Method 2.2.1 The VOF Method The Volume of Fluid (VOF) was first introduced by Hirt and Nichols [37] based on the MAC method, but with major improvements. VOF is based on defining one new variable, fractional volume value f, for each element. The fractional volume is set to have a value of 1 when the mesh cell is fully occupied by the fluid, 0 when the cell is empty with no fluid and a value between 0 and 1 when the cell is partially occupied by the fluid which represents a free surface (Figure 2.1). Figure 2.1 Distribution of ‘f’ at free surface and neighboring elements 14 The VOF method is based on three major steps. First is to calculate the fractional volume f of each cell at each time step. These fractional values are stored for each element and are updated after each time step to check if the state of an element has changed (e.g. fully occupied to partially occupied cells). The second step is to predict the shape of the free surface at the boundary cells. It is very important to distinguish between the location and the direction for the free surface at this step because for each location and fractional value combination, more than one free surface shape possibility exists. Finally, to construct the free surface based on the values of the fractional volumes and the shape of the free surface. The natural computational domain for each element is a 3x3 stencil with the element of interest being in the center of that stencil. Different VOF schemes utilize the 3x3 stencil in different manners. Some of the less common schemes use a different stencil size which will be discussed in the next section. The beauty of this method is that it does not depend on specific variables such as the fluid pressure, which allows for different techniques to be used in the second and third steps when trying to predict and capture the shape and location of the free surface. 2.2.2 VOF in the Literature With the use of the VOF method, the free surface problems became easier to deal with in fluid mechanics applications, but there are still some important points to consider with this approach. The main point is the smearing of f due to the averaging process and discontinuity which leads to loosing the interface between the elements [37]. This was particularly noticeable in the early development stages of the VOF method when a simple line interface calculation (SLIC) was utilized for free surface construction. In SLIC, the free surface interface is reconstructed using a sequence of segments that are aligned with the grid (either vertically or horizontally). The SLIC-VOF reconstruction method is relatively crude and would cause severe smearing of the fractional volume distribution. Therefore, the piecewise linear interface construction (PLIC) method was introduced to reduce the smearing and to smooth the representation of the free boundary. One of the critical aspects of the PLIC method is that 15 it does not perform the surface reconstruction as a series of joined segments for the whole free boundary, but as a discontinuous chain of segments with relatively small discontinuities. Figure 2.2 shows the averaging of the free surface using: (a) SLIC and (b) PLIC [42]. An improvement of the PLIC-VOF was suggested by Martinez et al. [43] to improve the curve of the constructed free surface segments. By incorporating the calculation of the free surface tension forces into the PLIC-VOF equations, a better curve estimation can be achieved. The use of this improved curvature technique eliminated the use of smoothed color functions to connect the chain of segments at the free surface. 1.01.0 1.0 0.0 0.1 0.4 0.1 0.875 0.4 1.01.0 1.0 0.0 0.1 0.4 0.1 0.875 0.4 (a) SLIC (b) PLIC Figure 2.2 Free surface averaging using (a) SLIC and (b) PLIC Because of the smearing problem and the difficulty in defining exactly the free surface, different free surface capturing techniques have been developed and implemented in the VOF method. One of the common methods in literature is the donor- acceptor method which is used to overcome smearing and to represent the free surface in a smooth manner [37, 44-47]. The concept behind the donor-acceptor method is that for each interaction between neighboring elements, one acts as a donor element as the fluid flows out of that element and the other one as an acceptor as the fluid flows into it. This 16 donor-acceptor scheme is achieved by calculating the change in the fluid volume fluxes between each cell and its neighbors, then average those values so that a sharp interface at the free boundaries can be maintained. Based on the distribution of f, and with the help of an orientation vector that points towards the free surface, the shape and location of the free surface may be determined. The donor-acceptor method can be also implemented in conjunction with a SLIC or PLIC or any other surface construction scheme. Nobari et al. [48] worked on a modified definition of the donor-acceptor method called the defined donating and accepting region (DDAR) with a SLIC-VOF scheme. In their work they developed an algorithm that uses a fully multi-dimensional boundary flux integration technique. The algorithm predicts the VOF flux to not only vertical and horizontal neighboring elements, but also to diagonal neighboring elements. The results obtained from their work indicate that it is necessary to predict the VOF advection for all of the eight neighboring elements to achieve maximum accuracy. The modified SLIC-VOF 8-node advection calculation lead to accurate results compared to some of the conventional PLIC-VOF methods. A major observation of the donor-acceptor method is that it is very difficult to implement in non-uniform grids, which limits it mostly to cases with uniformly constructed grids. Another interesting scheme when using the VOF method is the implementation of height functions in the construction of the free surface. Height functions have been used in 2D [49] and 3D [50] cases as well. In most VOF methods, a 3x3 stencil is used as the computational domain where the element of concern is located at the center and is neighbored by 8 elements. When height functions are employed, a 7x3 or 3x7 stencil is used for the calculation of the fractional volume and the free surface orientation. The idea is to calculate the normals of the free surface based on the horizontal and vertical fluxes, then, based on these values, if the vertical normal is bigger than the horizontal one, a 7x3 stencil is considered and vice versa. Fractional volumes are then summed along the direction that is more normal to the interface. Some methods start with the common 3x3 stencil to examine whether the initial estimate of the normal to the free surface is more horizontal or vertical before implementing a 3x7 or a 7x3 stencil [51]. This method has the advantage of having accurate curvatures and surface tension forces at contact lines in 17 addition to enforcing a contact angle boundary condition at a contact line because they utilize seven columns/rows instead of three. The use of the height functions is, however, more complex than using regular 3x3 stencils, and mesh refinement at the contact interface is usually needed to insure accurate results. In addition to that, height functions may lead to poor results if the stencil configuration includes a free surface other than the interface located at the element of interest. To overcome this problem, any partially occupied element (other than the element of interest) located in the 3x7/7x3 stencil is assumed to have a fractional volume value of one at that instance. After the computation is done for that stencil, the fractional volume value of that element is changed back to its original value. This scheme is also limited to uniformly constructed grids since elements of a 3x7/7x3 stencil will not be aligned on the same level in non-uniform mesh. The use of mesh refinement at the free surface is utilized in many VOF methods to minimize the need for smoothing color functions or modified free surface calculations. The approach is an adaptive grid refinement. The idea is to refine the mesh at the free surface by dividing each free surface cell (parent cell) into four child cells (also called baby or leaf cells), then keep dividing each free surface child cell into another four child cells until a user specified mesh refinement criterion is met [51]. Using adaptive grid refinement leads to better shape and more accurate curves at the free surface interface, depending on the level of cell division used in the computation. Mesh refinement can also be used with any other scheme (e.g. VOF-PLIC, donor-acceptor…etc) but it suffers from the same problems of the pertinent method and it is difficult to implement it into a non- uniform mesh. The reason behind this is because the element division process is usually set prior to the computation to a certain size which cannot be set for irregular element since they have different sizes. Additionally it is very challenging to transform the free surface interface from the parent cell to the child cells for a random grid. The algorithm is usually developed in a systematic fashion to ensure the same level of accuracy along the interface [52]. One of the main drawbacks of using adaptive grid refinement is that with every division step weendupwithwhatiscalled“hangingnodes” at the faces between neighboring elements of different sizes, in which a corner node from the child cell lies at the face-center of a larger cell (Figure 2.3). This is a problem because the connectivity of the FE is no longer satisfied. To overcome this issue, the larger cell is considered to have 18 two faces neighboring the child cells. Then, linear interpolation is used to estimate the values of the fluxes at the face-center in order to divide the fluid flow from the bigger cell into the smaller cells accurately. Although mesh refinement produces better results regarding the shape of the free surface, it is on the expense of the process variables computational accuracy since these variables must be interpolated back and forth between the original grid elements and the newly constructed ones[53]. Regular nodes Hanging nodes Figure 2.3 Schematic of regular and hanging nodes in adaptive grid refinement Another approach presented in the literature to capture the evolution of the free surface interface with the VOF method is to use an Eulerian-Lagrangian advection scheme [54]. In this method, the advection of the free surface interface is carried out using a Lagrangian step, and then the fractional volume values are calculated for the Eulerian grid. Based on the values of the fractional volumes, the fluid is remapped from the Lagrangian grid to the Eulerian grid. Mashayek and Ashgriz [55] discretized all the 19 cell wetting possible cases between the element of interest and its neighboring elements. Based on these cases, the nodes at the outer boundary get altered in a Lagrangian fashion to capture the fluid surface. Although this approach, by definition, can be easily implemented into non-uniform grids, it is unfeasible to summarize all the possible cases for an irregular grid. Later, Ashgriz et al. [56] developed a more general approach by moving the free surface elements, and not just the surface nodes, in a Lagrangian sense. The mesh does not get altered in this approach; therefore the interface after the advection step has to be mounted on the Eulerian grid. The fractional volume values must be determined at each neighboring element after the advection step. The free surface interface is constructed geometrically by knowing the coordinates of the Lagrangian elements after the advection and the original Eulerian elements. This modification allows this approach to be used in cases with non-uniform grid. Gueyffier et al. [57] developed a free surface interface tracking scheme for three dimensional flows using PLIC-VOF with a split Lagrangian advection scheme. They used a Lagrangian step for the advection of the free surface interface. Then, a PLIC technique is used to construct the free surface of the Eulerian grid based on the intersection between the Lagrangian and the Eulerian grids. Surface tension was implemented using the continuous surface stress or force method. This combination of the PLIC-VOF with the modified split Lagrangian scheme is very complex compared to other methods but can be upgraded to three dimensional problems without any significant increase in the difficulty. An Unsplit Lagrangian Advection (ULA) scheme for the VOF method was presented by Yang et al. [58]. Instead of using the operator split method, where more than one interface reconstruction might be required at each step, a complete and true Lagrangian advection is utilized in this approach. When the interface from one cell moves into the neighboring cells, then the flux between those cells is evaluated geometrically. In the ULA scheme, instead of calculating the outflow and inflow fluxes for the elements separately, the outflow flux for any cell is considered as an inflow for another cell. This way the conservation of mass is not compromised during the simulation. 20 Janssen and Krafczyk [59] presented a hybrid VOF algorithm that combines the advantages of different schemes. They used a PLIC-VOF with mesh refinement in conjunction with an averaging process called cell fill level approach. The goal of the averaging process is to determine, without using geometrical surface reconstruction, the filled face level (i.e. the fraction of a cell face that is wetted by the fluid) of a cell by averaging the fill level of the two cells between which fluid mass exchange takes place. This approach ensures mass conservation between two free surface elements but not between free surface and fully occupied elements. To overcome this issue, a case distinction is used to set the face fill level value to one if either of the two cells is a fully occupied, otherwise and average value is used between them. One of the common methods in computational physics is the level set (LS) method [60]. In recent years the coupling between the VOF and LS methods became very popular in CFD problems [61-63]. In the LS method, a smooth function (LS function) is used for the free surface interface representation. Generally, the LS function is defined as a signed distance function based on the shortest distance from a reference point to the interface where the sign indicates if the surface is continuous or dispersed. The function is reinitialized after each time step to update the location of the new interface. In the hybrid VOF and LS method, the LS function is used to calculate the interface normal and curvature, where VOF is used for the fluxes calculations and the free surface interface construction. The VOF method corrects the volume loss produced after reinitializing the distance after each step, where the LS function improves the accuracy of the interface curvature [64, 65]. Another approach used in the VOF method is the filling pattern [44, 45]. In this method the filling pattern follows the estimated direction of the free surface and starts gradually filling the cell with the fluid until it reach the desired fractional volume value. The filling scheme is very popular in the simulation of mould filling type of problems [66] since the filling pattern is, to some extent, steady and follows the same path during the whole process. Jeong and Yang [67] generalized the filling pattern to be used in more complex problems. They classified the free surface interfaces into eight different cases that cover all general possibilities. Each case is assigned a unique pattern depending of 21 the shape of the free surface. The filling process then fills the elements according to its pattern until the element is filled with the desired amount of fluid. To improve the results, an adaptive mesh refinement was employed in addition to a free surface smoothing function. Since the filling pattern is only used for the construction of the interface, it does not prevent the use of other VOF schemes which makes this approach very easy to modify and improve. 2.3 Transition from VOF to VOS One of the main reasons that the VOF method has been used extensively in the literature is because it allows for a wide range of schemes to be incorporated within its formulation. The main goal of this research is to apply a method that can handle the free surface of the unconstrained flow of the material in an Eulerian FE mesh for large deformation applications. The VOF method is used for the same purpose in CFD analyses; hence we propose the adaptation of the method into solid mechanics and therefore, would be named Volume of Solid (VOS) method. Rather than having a literature review of the VOF method and its development in a chronological sense, a general study of the various schemes and modifications was presented in the previous section. Understanding the advantages and disadvantages of the VOF methods presented in the literature is the basic ground for this work. Some of the methods might be superior to others in the field of fluid dynamics but that is not necessary true when adapting VOF into solid mechanics. Due to the difference in the mechanics between metal forming applications and fluid dynamic problems, certain methods can be utilized in a more advantageous way. Understanding the dynamics of metal forming would also lead to new modifications in the VOF method that may improve the efficiency of the method. In Chapter 4, a detailed discussion of the proposed Volume of Solid (VOS) method is presented. The calculation of the fractional volume values is presented in a similar way to the available VOF method in the literature and no major modification is needed when transitioning to VOS. A new method for identifying the direction and orientation of the 22 free surface interface is proposed in this work by fully utilizing a 3x3 stencil in a way similar to the height functions. The construction of the free surface is developed for a uniformly constructed Cartesian grid and for a general case of non-uniform grid. For the uniform grid case, a classification of all the possible free surface shapes is used with a PLIC scheme whereas a modified filling scheme is developed for the non-uniform case. The implementation of the VOS method is presented in Chapter 5 with a discussion on all related issues such as the connectivity of the free surface, handling the contact and tracking the material point related properties. 23 CHAPTER 3 EULERIAN FORMULATION In this section we present the governing Eulerian virtual work equilibrium equations used in finite element analysis of quasi-static and dynamic systems. The derivation will follow the standard notation used by many researchers, e.g. [3]. The derivation presented in this work is based on the ALE formulations developed by Bayoumi and Gadala [68, 69]. In any finite element formulation there are two sets of configurations, one refers to the material particles and the other one refers to the grid points. The two sets share the same velocity in Lagrangian formulation and therefore one set of equations is required. In ALE formulation, since the mesh moves independently from the material, the two configurations have different velocities and, therefore, two sets of equations are required. Since we are using Eulerian formulation in this work where the grid is fixed in space, and hence has a zero velocity, we only consider the material velocity. When deriving the Eulerian equations, time derivative of several quantities will be required. There are two different time derivatives in finite element analysis, one is the material time derivative and the other is the grid time derivative. The material time derivative is required to track the change in the function with time for the material point. For any arbitrary function f, the material time derivative is defined as the rate of change of the function while the material particle is fixed: m i t t X f f t (3.1) Another definition for the function rate of change is required as the grid point is held fixed. This is essential to track the value of the function within the grid as material particles move with time. The grid time derivative is defined as the rate of change of the function while the grid point is fixed: 24 g i t t X f f t (3.2) The relation between the material and grid time derivatives can be written by simply defining the progression of function f according to the following transport equation: t t t t i t i f f f v x (3.3) where t vi is the material point velocity at time t. This relation is similar to the relation used to derive ALE equations [70] except for the grid velocity which is set to zero in Eq. (3.3). The formulation of the Eulerian equation in this work is divided into two sections. In the first section, the Eulerian equilibrium equation is derived for quasi-static analysis which concerns cases where inertia effects may be neglected as in low speed metal forming processes. In the second section, the derivation is extended to dynamic analysis where the effect of the inertia forces is included into the equilibrium equation. 3.1 Quasi-Static Analysis In quasi-static analysis, the natural way to derive and write the equilibrium equation is to use the principle of virtual displacements since the inertia effects are not considered. The equilibrium of the body at time t t using the principle of virtual displacements is written as: t t t t t t t t ext ij t t ij V e dV W (3.4) where represents the virtual or the variation of a variable, t t ij are the Cauchy stress tensor components at time t t and t t ije is the conjugate strain tensor which is defined as: 25 1 2 ji t t ij t t t t j i uu e x x (3.5) The external virtual work in Eq. (3.4) is expanded by writing the components of the body force and the surface traction: t t t t t t ext t t t t B t t t t S t t i i i i V S W f u dV f u dS (3.6) where t t B if and t t S if are the components of the body force per unit mass and surface traction at time t t , respectively. Once the variables at time t t are computed then Eq. (3.4) can be expanded to get the complete Eulerian equilibrium equation. 3.1.1 Incremental Decomposition Since we are referring variables to the grid configuration, we assume that the value of any variable at time t t is defined to be its value at time t plus an increment that is calculated by multiplying the grid time derivative by the time increment t . Following the stated assumption, the material density is written as: t t t t t (3.7) The conservation of mass at time t is given by: t t t i t i v x (3.8) Using Eq. (3.8) in Eq. (3.3), the continuity equation is written as: t t t t ti it t i i v v x x (3.9) 26 Now we substitute Eq. (3.9) in Eq. (3.7) to get the full decomposition of the material density: t t t t t k kt t k k u u x x (3.10) Similar to Eq. (3.7), the stress components at time t t are expressed as: t t t t ij ij ij t (3.11) Using Eq. (3.3), we can expand Eq. (3.11) to the following: t ijt t t t ij ij ij k t k t u x (3.12) The material rate of Cauchy stresses t ij is calculated from the material constitutive relation: ˆ t t t ij ijkl klC D (3.13) where t ijklC is the forth order tensor of the material properties, and t klD is the rate of deformation tensor given by: 1 2 tt jt i kl t t j i vv D x x (3.14) ˆt ij is an objective stress rate tensor (Truesdell, Green-Naghdi, Jaumann…etc).In our case the Truesdell stress rate will be used based on the comparison presented by Gadala and Wang [71]: tt t jt T t t t tk i ij ij ij ik jkt t t k k k vv v x x x (3.15) 27 The variation in the strain components t t ije can be derived in a similar way: t t ij t ij t ije e e t (3.16) Since the grid in Eulerian analysis is fixed in space with zero velocity, the grid time derivative of the variation in the strain components t ije is equal to zero, so Eq. (3.16) is reduced to: t t ij t ije e (3.17) The incremental decomposition of the elemental volume and elemental surface area at time t t is derived exactly as the variation in the strain components as there is no change in volume or surface area in an Eulerian element: t t tdV dV (3.18) t t tdS dS (3.19) 3.1.2 Linearization After expanding the variables using the incremental decomposition, substituting these decompositions into the equilibrium equation, Eq. (3.4), and neglecting all incremental quantities with higher orders; we can end up with the linearized equation. First we will expand the LHS of Eq. (3.4), by using Eq‟s. (3.12), (3.17) and (3.18): t t t t t t t t t t t t t ij t t ij ij t ij ij t ij V V V t ij t k t ijt kV e dV e dV t e dV u e dV x (3.20) The external virtual work in Eq. (3.6) consists of two parts, body forces and surface tractions. The surface traction term of Eq. (3.6) can be expressed by Eq. (3.19): 28 t t t t t S t t t t S t i i i i S S f u dS f u dS (3.21) The body force can be expanded and linearized using Eq‟s. (3.10) and (3.18): t t t t t t t t t B t t t t t B t i i i i V V t t t t B t t t B tk i i i k it t k kV V f u dV f u dV u f u dV f u u dV x x (3.22) 3.1.3 Treatment of Convective Terms Convective terms refer to the spatial derivative of variables t t if x as the last term in Eq. (3.3). The presence of convective terms in any equation complicates the numerical solution. One of the common methods used to deal with convective terms in Eulerian FEA is the operator split method. In operator split method, Eq. (3.3) would be divided into two equations, non-advective and advective equations. In this work we will modify the method presented by Bayoumi and Gadala [68, 69] for dealing with convective terms in ALE formulation to utilize it in Eulerian formulation. The advantage of this method is that it does not require two sets of equations to be solved separately as in the operator split method, in addition to avoiding the need to compute the spatial derivative of the variables. The first convective term that we have to deal with is the last term in Eq. (3.20) which requires the computation of the spatial derivative of the stresses. Because the stresses are usually evaluated at the integration points rather than the nodal points, the stress field is generally discontinuous along the element edges. This leads to the stress gradients not being reliably computed on the element level when evaluating the element matrices. That is where methods such as the operator split method are used to overcome this problem. In this formulation, we will use Gauss divergence theorem to deal with the convective term. First we re-write the last term in Eq. (3.20) as follows: 29 t t t t tt k ij t ijij t ijt t t t k t ij k ijt t t k k kV V V t tk ij t ijt kV u e e u e dV dV u dV x x x u e dV x (3.23) Since we use the VOS method to construct the interface at the free surface elements, then at that point the material in the new free surface elements will be given the velocity values through the surface nodes, thus making both, the material and the grid, share the same velocity at that instant. This is because in the VOS method we define the free surface within the element cell that contains the free surface, making the material and the grids share the same surface/boundary all times. So to ensure that the material and the grid have the same boundary at all times, we write the following boundary constraint by enforcing a fictitious grid boundary velocity: on the boundary 0gt t ti iiv v n (3.24) Now we include the grid displacement term in the first term on the RHS of Eq. (3.23) and apply the divergence theorem with the use of Eq. (3.24): 0 t t g t k k ij t ij t g t t t k k ij t ij kt kV S u u e dV u u e n dS x (3.25) Substituting Eq‟s. (3.23) and (3.25) into (3.20), the internal virtual work is written as: t t t t t t t t t t t t t t ij t t ij ij t ij ij t ij V V V t ijt t t tk k ij ij t ijt t k kV V e dV e dV t e dV e u u dV e dV x x (3.26) Now the other convective term that we have to deal with is the convective body force term which is the last term in Eq. (3.12). Re-writing the body force convective term: 30 t t t t t t t B t t i k it t B t t i k it t k kV V t t B t t t t B t ti k k i i it t k kV V t t B t ti i k t kV f u u f u u dV dV x x f u u u dV f u dV x x u f u dV x (3.27) Applying the divergence theorem to the first term on the RHS of Eq. (3.27) and using the boundary constraint in Eq. (3.24) we get: 0 t t t t B g t i k ik t t t B g t t t i k i kkt kV S f u u u dV f u u u n dS x (3.28) Substituting Eq‟s. (3.27) and (3.28) into (3.22), the body force is written as: t t t t t t t t t B t t t t t B t i i i i V V t t B t t t t t B ti i k i i kt t k kV V f u dV f u dV f u u u dV f u dV x x (3.29) 3.1.4 Eulerian Equilibrium Equation for Quasi-static Analysis To write the final form of the Eulerian equilibrium equation, we substitute the internal virtual work (3.26), the body force (3.29) and the surface traction (3.21) into Eq. (3.4). This gives us the final form of the principle of virtual displacements at time t t referred to time t: t t t t t ijt t t t t tk ij t ij ij t ij k ijt t k kV V V t t ext t t ij t ij V eu t e dV e dV u dV x x W e dV (3.30) where the external virtual work is: 31 t t t t t B t t ext t t t B ti i k it kV t t t B t t t S ti i k i it kV S f W f u u dV x u f u dV f u dS x (3.31) Finally we use the constitutive relations introduced in Eq‟s. (3.13), (3.14) and (3.15) to the first integral in Eq. (3.30) gives us the final form of the equilibrium equation: t t t t t t t t t ijkl t kl t ij ij t ij V V t ijt t t tk i k ij ijt t t k j kV V t t ext t t ij t ij V C e e dV dV e u u u dV dV x x x W e dV (3.32) where 1 2 k k t ij t t i j u u x x (3.33) Eulerian equation (3.32) is linear in the incremental displacement iu and is referred to an Eulerian grid since the velocity and displacement of the grid does not show in the equation. 3.2 Dynamic Analysis In this section we expand the derivation of the quasi-static Eulerian equilibrium equation to dynamic analysis by including the inertia effects. By including the inertia forces we will have to deal with the material time derivative of the material velocity t t iv , and the grid time derivative of the material velocity t t iv . As the Eulerian mesh has zero velocity, we end up with the time derivative of just the material velocity, and for 32 clarity we will refer to the referential material acceleration (grid time derivative of the material velocity) as t t ia . 3.2.1 Virtual Work Done by Inertia Forces Using the relation between the material and grid time derivatives in Eq. (3.3), we can write the virtual work done by the inertia forces as: t t t t t t t t t t t t t t t t t t i i i i V V t t t t t t t ti j it t jV v u dV a u dV v v u dV x (3.34) The first term on the RHS of Eq. (3.34) is the referential inertial term, where the second term is the convective inertia term. Once we linearize both terms we then can incorporate them in the Eulerian virtual work Eq. (3.32). 3.2.2 Decomposition of Velocity and Acceleration The material velocity and acceleration at time t t can be related to the previous value at time t by simply using the relations: t t ti i ia a a (3.35) t t ti i iv v v (3.36) where the incremental acceleration and velocity ia and iv depend on the time integration scheme that is used in the analysis. Note that higher orders of these quantities will be neglected to uphold a linear increment in their values. 33 3.2.3 Linearization of the Referential Inertia Term Following the same procedure used in quasi-static analysis, using the incremental decompositions and linearization gives us the following: t t t t t t t t t t t t t t t i i i i V V t t t t t t t t t tk i i k i it t k kV V a u dV a u dV v a u t dV v a u t dV x x (3.37) The last term on the RHS of Eq. (3.37) can be written as: t t t t t t t t t k i it t t t t k i it t k kV V t t t i it t t t t t tk i i kt t k kV V v a u v a u t dV t dV x x a uv a u t dV v t dV x x (3.38) Applying Gauss divergence theorem to the first integral in RHS of Eq. (3.38) and using Eq. (3.24) we get: 0 t t t gt t t k i ik t t kV t gt t t t k i i kk S v v a u t dV x v v a u n t dS (3.39) Substituting Eq‟s (3.38) and (3.39) in Eq. (3.37), we get: t t t t t t t t t t t t t t i i i i V V t t i it t t k t kV a u dV a u dV a u v t dV x (3.40) Using Eq. (3.35), we can write the referential inertia term as follows: 34 t t t t t t t t t t t t t t t t i i i i i i V V V t i it t t k t kV a u dV a u dV a u dV a u v t dV x (3.41) 3.2.4 Linearization of the Convective Inertia Term The convective term (the second term on the RHS) in Eq. (3.34) can be expanded into the following: t t t t t t t t t t t t t t t t t t ti i j i j it t t j jV V t t t t t t tk i j it t k jV t tt t t t ti k j it t k jV v v v u dV v u dV x x v v v u t dV x x v v v u t dV x x (3.42) We-writing the third term on the RHS of Eq. (3.42): t t t t t tt t t t ti k j it t k jV t t t t t t i k j it j t t kV t t t t t t tk i j it t k jV t t t t jt t ti k it t k jV t t t t t t i k j it t k j v v v u t dV x x v v v u x t dV x v v v u t dV x x v v v u t dV x x v v v u x x t t V t dV (3.43) 35 Applying Gauss divergence theorem and the boundary constrain to the first integral in RHS of Eq. (3.43) we get: 0 t t t tt gt t t i k j ik t j t t kV t tt gt t t ti k j i kk t jS v v v v u x t dV x v v v v u n t dS x (3.44) Substituting Eq‟s. (3.43) and (3.44) into Eq. (3.42) we get: t t t t t t t t t t t t t t t t t t ti i j i j it t t j jV V t t t t jt t ti k it t k jV t t t t t t ti k j it t k jV v v v u dV v u dV x x v v v u t dV x x v v v u t dV x x (3.45) Using Eq. (3.36), we write the convective inertia term as follows: t t t t t t t t t t t t t t t t t t ti i j i j it t t j jV V t t t t t ti i j i j it t j jV V t t jt t ti k it t k jV t t t t ti k j it t k jV t t t i k j t t k v v v u dV v u dV x x v v v u dV v u dV x x v v v u t dV x x v v v u t dV x x v v v x t t i jV u t dV x (3.46) 36 3.2.5 Eulerian Equilibrium Equation for Dynamic Analysis Incorporating Eq‟s. (3.41) and (3.46) into Eq. (3.32) we get the Eulerian equation of motion that includes the inertia effects: t t t t t t t t t t t t t t i i ijkl t kl t ij ij t ij V V V t t t t t ti i j i j it t j jV V t t t ti k j it t k jV t ijt t t tk i k ij ijt t t k j kV V t t ext t ij a u dV C e e dV dV v v v u dV v u dV x x v v v u t dV x x e u u u dV dV x x x W t t t t t t t t t t t ij i i V V t ti it t t t t ti k j it t k jV V t t jt t ti k it t k jV t t t t ti k j it t k jV e dV a u dV a u v v t dV v u dV x x v v v u t dV x x v v v u t dV x x (3.47) Equation (3.47) is referred to an Eulerian grid that is fixed in space and has zero velocity and acceleration. The equation is linear in the incremental displacements ui, incremental velocities vi and incremental accelerations ai. 37 3.3 Discretization of the Eulerian FE Equations 3.3.1 Discretization of the Quasi-Static Eulerian Equation The matrix form of the Eulerian quasi-static equilibrium equation (3.32) is written as follows: 1extK K u f fi it E t A t t t t (3.48) where 1 2K K Kt E t E t E (3.49) 1 2K K Kt A t A t A (3.50) 1Kt E and 2Kt E are the Eulerian material and geometrical stiffness matrices given by: 1 1 1K B C B d t T t E E t E t t t V V (3.51) 2 2 2 2K B S B d t T t E E t E E t t t V V (3.52) 1Kt A and 2Kt A are the convective stiffness matrices given by: 1 1 1K B S H d t T t A A t A t t V V (3.53) 2 2 2 2K B S B d t T t A A t E E t t t V V (3.54) 1 f it t is the internal force vector given by: 1 11f B σ d t Ti it t E t t t t t t V V (3.55) 38 and extft t is the vector of the externally applied loads. Matrix Ct is the elastic-plastic material constitutive matrix. Defining kh to be the element shape function corresponding to nodal point k and N is the number of element nodal points, matrices H, 1BE , 2BE , 1BA , 2BA , 2St E and 1St A are given by: 2 2 0 H 0 k k N h h (3.56) 1 1 4 2 0 0 B 0 k t k t E k kt t t k N t j j j N h x h y h h y x h h x (3.57) 2 1 5 2 0 0 0 B 0 0 k t k t k tE t k t k N t j j j N h x h y h x h y h h x (3.58) 39 2 2 2 2 2 2 2 1 2 2 2 2 8 2 0 0 0 0 B 0 0 1 0 1 0 x t x t x t t x t A t x t x t t x k t t t x t t N h x h y h x y h x h y h x y h h x x x h x y (3.59) 2 1 5 2 0 0 0 B 0 0 k t k t k tA t k t k N t j j j N h x h x h y h y h h x (3.60) 40 2 5 5 0 0 0 0 0 0 S 0 0 0 0 0 0 0 0 0 0 t t xx xy t t xy yy t E t t xx xy t t xy yy t zz (3.61) 1 2 8 0 0 0 S 0 0 0 T t t t t t xx xy xy yy zzt A t t t t t xy xx yy xy zz (3.62) 1σ it t is the stress vector defined by: 1 1 1 1 1 σ σ σ σ σ i i i i it t t t t t t t t t xx yy xy zz (3.63) 3.3.2 Discretization of the Dynamic Eulerian Equation The matrix form of the Eulerian dynamic equilibrium equation (3.47) is written as follows: 11 2 ext 1 11 3 4 M a C C v K K u f f M M a C C C v i i i it E t A t A t E t A t t t t i it E t A t t t A t A t A t t (3.64) Matrix Mt E is the Eulerian mass matrix given by: M H H d t t E t T t V V (3.65) 41 Matrix Mt A is given by: 2 1,2 1 2 1,2 2 ,2 1 2 ,2 2 2 M t A t A i j i jt A t A t A i j i j N N M M M M (3.66) in which i and j indicate node numbers from 1 to N, and 1 2 1,2 1 2 ,2 1 d t N j ti i j k xkt t kt A t A t t i j i j N j tV i i j k ykt t k h h h h h v x x M M t V h h h h h v y y (3.67) 2 1,2 2 ,2 1 0 t A t A i j i jM M (3.68) Matrix 1Ct A is given by: 1 1 2 1,2 1 2 1,21 1 1 2 ,2 1 2 ,2 2 2 C t A t A i j i jt A t A t A i j i j N N C C C C (3.69) in which 1 12 1,2 1 2 ,2 1 1 d t N N j jt A t A t t t t i j i j i k xk k ykt t k kV h h C C h h v h v V x y (3.70) 1 1 2 1,2 2 ,2 1 0 t A t A i j i jC C (3.71) 42 Matrix 2Ct A is given by: 2 2 2 1,2 1 2 1,22 2 2 2 ,2 1 2 ,2 2 2 C t A t A i j i jt A t A t A i j i j N N C C C C (3.72) in which 2 2 1,2 1 1 2 2 1,2 1 2 2 1,2 1 1 2 2 1,2 1 d d d d t t t t N t A t t tk i j i j xkt kV N t A t t tk i j i j xkt kV N t A t t tk i j i j ykt kV N t A t t tk i j i j ykt kV h C h h v V x h C h h v V y h C h h v V x h C h h v V y (3.73) Matrix 3Ct A is given by: 3 3 2 1,2 1 2 1,23 3 3 2 ,2 1 2 ,2 2 2 C t A t A i j i jt A t A t A i j i j N N C C C C (3.74) in which 3 3 2 1,2 1 2 ,2 1 1 1 1 1 1 d t t A t A i j i j N N N j jt t t tk k i xk yk k xkt t t t k k kV N N N j jt t t tk k xk yk k ykt t t t k k k C C h hh h h v v h v x x y x h hh h v v h v t V y y x y (3.75) 43 3 3 2 1,2 2 ,2 1 0 t A t A i j i jC C (3.76) Matrix 4Ct A is given by: 4 4 2 1,2 1 2 1,24 4 4 2 ,2 1 2 ,2 2 2 C t A t A i j i jt A t A t A i j i j N N C C C C (3.77) in which 4 4 2 1,2 1 2 ,2 22 2 1 2 1 1 2 2 1 2 t t A t A i j i j N j jt ti k i xkt t t t kV N N j j j t ti i i k xk k ykt t t t t t k k N j j ti i k ykt t t k C C h hh h h v x x x x h h hh h h h v h v x y y x x y h hh h h v y y y 2 dtt V (3.78) 4 4 2 1,2 2 ,2 1 0 t A t A i j i jC C (3.79) 44 CHAPTER 4 VOLUME OF SOLID (VOS) METHOD Finite element Eulerian formulation in metal forming has mostly been used to simulate steady state problems [5-8]. Only limited trials in simulating unsteady state problems may be found [9-14]. In these trials, only the middle stage of the forming process is modeled and the expected shape of that stage is meshed with Eulerian grid. This pseudo analysis approach has been considered to deal with the problem of free surface detection and identification. The drawback is that it is limited to certain number of processes, and it does not give a full picture of the whole process. Using Eulerian formulation to simulate all stages of unsteady state metal forming process suffers from the detection of the free surface and the treatment of the convective terms in the analysis. In this section we discuss the development of the VOS theory and its application to problems in solid mechanics and metal forming. Also, we discuss the integration of the method in a FE context for a uniform grid. First, we define a function f which represents the fractional volume of the solid in each cell that has a value between zero (at an empty cell) and one (at a cell fully occupied by solid). A value of f between zero and one signals a cell with a free surface (Figure 2.1). Hirt and Nichols [37] presented a transport equation for the evolution of the function f in a VOF method as: 0 t t t i t i f f v x (4.1) The above equation states that f moves with the material and, in a Lagrangian formulation, it indicates that f remains constant in each cell. Also, it indicates that the material time derivative of f is zero as given by Eq. (3.3). Solving Eq. (4.1) gives the distribution of f and locates the elements in which free surface occurs. The next step would be to find the direction in which the free surface is located in each element. Different techniques were proposed in literature and they mostly present some 45 modification to the basic method of VOF-PLIC. Some work attempted solving equation (4.1) directly [31, 32] by considering f as an unknown variable that is discretized on the computational grid. Other methods are based on tracking the change of f through the fluid fluxes [37, 43]. One of the common methods used is the cubic-interpolated propagation method (CIP). This method is based on the fact that a variable and its spatial derivative propagates with the velocity and the spatial profile within each grid and is interpolated with a cubic polynomial [54]. Some methods use a semi-Lagrangian conservative scheme to solve the transport equation and to find the shape of the free surface. In addition to these techniques, mesh refinement and re-meshing techniques were also introduced to VOF solutions at the free surface elements [49-56]. The numerical solution of equation (4.1) tends to give some smearing problems; the reason being that the gradient of f should be physically singular at the boundary but is finite in the numerical solution. So we need a method that can specify the shape and direction of the free surface in a consistent and continuous way along the element boundaries. Once the free surface is defined along the boundary, we need a scheme to update the current free surface to the new one in the Eulerian grid. The theory of the VOS method is presented through three major steps. First, equation (4.1) is solved in order to obtain the distribution of the fractional volume value f on the domain. Next, a direction vector that is normal to the free surface interface and points toward the free boundary is calculated based on the distribution of f that was obtained in the previous step. In the final step, the free surface interface is constricted by combining the free surface directional vector and the values of f at each element. The presentation in this chapter is divided into two parts. The first part presents the theory of the VOS method for a uniformly constructed grid, whereas the second part extends the method for a generalized non-uniform grid. The first step of the VOS method is almost the same between the two parts except for some minor differences. The calculation of the free surface directional vector in the second part follows the same approach of the uniform mesh but with a major modification in its equation in order to extend it for a generalized case of a non-uniform mesh. The construction of the surface is where the two parts differ. For the uniform mesh, we present a method based on the PLIC 46 scheme and a classification of all the possible interfaces. For the non-uniform mesh, we present a new scheme based on the filling pattern that can be used regardless of the shape of the element. 4.1 VOS for a Uniform Mesh 4.1.1 Solution of the Transport Equation In order to find the location of the free surface we first need to calculate the fractional volume of each element. In what follows, this is calculated based on the volume change between each element and its neighboring elements. In solving Eq. (4.1), f is defined for the cell or the element, not the node. This means that in contrast to the discretized virtual work equation, another scheme is required to integrate the transport equation of f. In the present analysis, the equation is first integrated over each element (cell) as follows: 0 t t t i t iV f f v dV x (4.2) Separating the integrals and applying the divergence theorem, Eq. (4.2) becomes: t t t i i A f V f v n dA (4.3) Explicit time integration for Eq. (4.3) has been implemented in our analysis. Applying the central difference method, Eq. (4.3) becomes [44]: v n j new old i i j i j t f f f A V (4.4) j solid j j A f A (4.5) 47 where i is the element number and j is the element face number, new if and old if are the fractional volume values for the i th element at the new and old time steps, Vi is the i th element volume, Aj is the j th face of the i th element, v and n are the velocity and unit outward normal vectors at the j th face of the i th element. Figure 4.1 shows the definition of the parameters of Eq‟s.(4.4) and (4.5). Since the elements here are considered as control volumes, the fractional volume calculation is based on the facial net change in each element, which is presented by the last term in Eq. (4.4). Also, j f , defined as the facial fractional value, indicates what part of the j th face of the i th element is in contact with the material and is defined by the fraction between the facial area Aj solid in contact with the material and the total facial area Aj as shown in Figure 4.1. Flow Velocity i th element n j th face Figure 4.1 Definition of free surface element parameters 4.1.2 Free Surface Directional Vector One of the main problems in most methods used to find the location of the free surface is that they do not ensure a consistency of the free surface among all the free surface elements. This means that the free surface can be calculated accurately for each 48 element separately, but these free surfaces may not be connected. Solving Eq. (4.4) gives value for f in each element, but there are many cutting planes within an element that satisfies a given f value. Eliminating this problem is based on the idea that a normal vector together with the fractional volume value of each element determines a unique planar interface cutting the element (cell). Different methods have been implemented to calculate the free surface normal vector [44, 45, 51, 58, 66]. These methods differ with respect to the influence of the fractional volume distribution in the sub-domain around the element of interest. The most common sub-domains used in the literature are 3x3, 5x3, 3x5, 7x3 and 3x7 stencils. All of these, except of the 3x3 stencil, are mostly used with height function or line segment methods to calculate the vector normal to the free surface. In this study we use the 3x3 stencil but with a technique that captures the advantages of the height function and the regular 3x3 method of calculation. To ensure optimum connectivity of the free surface, we calculate the location for all the free-surface elements simultaneously with the same level of accuracy. To do this, we first calculate the direction of a vector that is perpendicular to the free surface, and then construct the free surface for all the elements simultaneously. In most of the methods that use a 3x3 stencil, the directional vector is calculated by computing the horizontal and vertical fluxes of f between the element of interest and the neighboring elements in the horizontal and vertical directions, respectively. This is considered an acceptable approximation but it takes into account only two neighboring elements in each direction which may give rise to inaccurate indication of the directional vector. In methods that use the height function, a 3x7 or a 7x3 stencil is used, depending on the angle of the free surface, and the number of utilized neighboring elements increases. This bigger stencil gives better results for the interface curvature than the traditional 3x3 stencils but has the drawback of assuming that one component of the directional vector is a unit normal to the side with the 7 elements. Therefore, we need to develop a procedure that calculates the directional vector more accurately without such assumptions. In this analysis we calculate the direction of the vector by taking into account the fractional volume values for each neighboring element (horizontally, vertically and diagonally). We present a new approach that computes the fluxes in a similar way to the 49 height function method but with some modifications. First we use a 3x3 stencil with the element of interest being at the center of the stencil. Second, we utilize the horizontal and vertical directions simultaneously without the assumption of the unit normal that is used in height functions. Therefore, we are measuring the fluxes in f between neighboring rows/columns rather than neighboring cells. This gives a better picture for the distribution of the fractional volume values and their effect on calculating the directional vector while using a smaller stencil size than the height function in the same time. For our 3x3 stencil, we consider the element of interest is located in row k and column l. The fractional volume difference between the three vertical elements on the l+1 column and the three on the l-1 column is used for the horizontal component of the direction vector, and vice versa. The free surface direction vector is calculated as follows: 1 1 , 1 , 1 1, 1, 1 1 ,i i i R R n n R r k l m l m l x k n k n y m k n l i f f f f (4.6) where r is calculated for the i th element, f is the fractional volume of the element in the k th row and l th column depending on the summation, and nx , ny correspond to the facial normal vectors in the x- and y-direction, respectively. Note that Eq. (4.6) is carried out locally rather than globally, so that for each element we focus on the 3x3 stencil and perform the calculation within that stencil only. The process is then repeated for all the elements along the free boundary. Because interior elements have a fractional volume that is larger than that of elements at the surface, making the free surface directional vector point to the inward direction, the negative sign is used to make the directional vector point to the outward direction of the interface. Equation (4.6) indicates that the directional vector is normalized and is based on the distribution of the fractional volume along the domain. Figure 4.2 shows the definition of the directional vector. This procedure is applied for a uniform grid with equal size square elements. The advantage is the ease of mesh generation and the ease of calculating r mathematically, since the normal vectors are set before we start the solution and defined as i, j, -i or –j. 50 In some situations, the free surface element might be located at one end of the grid domain. This means that the element is located next to a fictitious wall and might be in contact with a fixed rigid object, and in this case the 3x3 stencil will be missing three elements (either horizontal or vertical). To overcome this problem and to ensure having a 3x3 stencil all the time, we assume to have a fictitious element on the side adjacent to face j of the free surface element, where the grid domain ends, with two fictitious elements neighboring it in the other opposite directions. We assume that this fictitious element has a value of f that follows a linear variation with the adjacent elements. If the linear progression of f leads to an overshoot of the fractional volume value then an f value of one is used for that element and the amount of the overshoot is assigned to the fictitious element neighboring it in the interface progression direction. The other fictitious element is assigned an f value of one. If the linear progression of f leads to undershoot of f then a value of zero is used for both fictitious elements and the undershooting is applied to the remaining fictitious element that neighbors it in the direction if the interface progression. Considering the case where the free surface is adjacent to the right end of the grid domain, the construction of the fictitious elements is performed as follows: 1, 1 , , , 1 , , , 1 , 1 1, 1 1, 1 , , , 1 , 1 , , , 1 1, 1 , 1 if 1 then 1 1 0 if 0 1 then 1 if fic k l k l k l k l fic k l k l k l k l fic k l fic k l fic k l k l k l k l k l k l k l fic k l k l k f f f f f f f f f f f f f f f f f f f f 1, 1 , , 1 , 1 1, 1 , , , 1 0 0 then 0 1 fic k l fic l k l k l fic k l k l k l k l f f f f f f f (4.7) Figure 4.3 shows the three different cases for the construction of the fictitious elements based on the linear progression of f. The three cases correspond to the case 51 where a fictitious column is missing in right end of the domain. Other cases where fictitious rows/columns are needed at other sides of the domain follow the same rule used in Eq. (4.7) and shown in Figure 4.3. i th element k-1 k k+1 l-1 l l+1 nx r ny fk-1,l+1 Figure 4.2 Definition of the free surface directional vector (b) Within range(a) Overshoot (c) Undershoot Figure 4.3 Construction of fictitious elements at the right end of the domain 52 4.1.3 Construction of the Free Surface using a Geometrical PLIC Scheme Based on the value of f and the free surface directional vector, we advance to locate the free surface. In the following analysis we divide the possible free surface configurations into four main categories, depending on the quadrant the free surface directional vector pointing to. Figure 4.4 shows the classification of the four categories based on the free surface directional vector r. In the following discussion we only cover the case with directional vector components, rx and ry, having positive signs. The rest of the quadrants follow the same procedure. In any quadrant there are three interface shape cases. Two of those cases are unique and the third one is a general case that consists of several possibilities. These cases are identified by the angle 𝜃 = tan−1 𝑟𝑦 𝑟𝑥 as illustrated in Figure 4.4. The general case has 0 o < θ < 90o and the unique cases have θ equal to 0 o or 90 o . θ θ r+ry +rx θ θ r +ry -rx θ θ r-ry +rx θ θ r -ry -rx Figure 4.4 Four quadrant categorization of the free surface 53 Based on the value of the angle θ and the value of the fractional volume f we can define the location and shape of the free surface. Basically, we calculate the area of the element that is filled with the material and use the angle to specify which part of the element is filled. Consider the general case; there are six possible situations to fill the element as shown in Figure 4.5 (cases 1-3, 5-7). The lengths a1 and a2 are two unknown variables that locate where the free surface intercepts the element edges and are calculated for each case depending on θ and f. Cases 1, 2 and 3 are considered when the free surface cuts through the element at an angle θ ≤45o. Cases 5, 6 and 7 are considered when the free surface cuts through the element at an angle 45 o < θ < 90o. The difference between cases 1 & 3 and cases 5 & 7 is that in cases 1 & 3 a1> a2, while cases 5 & 7 have the opposite. The difference between cases 2 and 6 is that in case 2 the free surface cuts through two element sides that are horizontally opposite to each other, while in case 6 the free surface cuts through the vertically opposite edges. Since we know which set of cases we have (1, 2, 3 or 5, 6, 7) from the identified angle, we need to specify which case among the set is applicable for each element. To do that, we define two critical fractional volumes 𝑓𝑐𝑟 1 and 𝑓𝑐𝑟 2 for the first set and two critical fractional volumes 𝑓𝑐𝑟 3 and 𝑓𝑐𝑟 4 for the second set, based on the area and geometry of the material within the cell: 1 2 3 4 1 1 0 45 tan , 1 tan 2 2 1 1 45 90 cot , 1 cot 2 2 cr cr cr cr f f f f (4.8) By comparing the actual fractional volume value f with the critical values defined in Eq. (4.8) we choose the specific case from Figure 4.5 as follows: 1 1 2 2 Case 1 for 0 45 Case 2 Case 3 cr cr cr cr f f f f f f f (4.9) 54 3 3 4 4 Case 5 for 45 90 Case 6 Case 7 cr cr cr cr f f f f f f f (4.10) The last two special cases are straight forward, for θ = 00 and for θ = 90o (cases 4 & 8 in Figure 4.5). Now that we know the specific case for each element, we need to calculate the variables a1 and a2 to construct the free surface. This is done by calculating the area that is filled by the material and by relating the two variables through the angle θ. If Δs is the element edge length, then by using f, θ and Δs we can generate the free surface using a piecewise linear interface construction (PLIC). Considering cases 1 and 5 as an example, the a1 and a2 are calculated as follows: 2 11 2 2 1 1 2 tan 2 cot tan1 2 a a s fa a a a a s f (4.11) Similarly, the rest of the cases follow the same procedure. Table 4.1 summarizes the equations used in the calculation of a1 and a2 for all cases. θ θ a2 θ Case 1 Case 3Case 2 θθ θ Case 5 Case 7Case 6 Case 4 Case 8 a1 a1 a1 a1 a2a1 a1 a1 a2 a2 a2 a2 a2 Figure 4.5 Free surface shape possibilities in a uniform mesh 55 Table 4.1 Summary of the free surface interface parameters calculations Relation between parameters Equations Case 1 𝑎2 𝑎1 = tan 𝜃 1 2 𝑎1𝑎2 = ∆𝑠 2𝑓 𝑎1 = ∆𝑠 2𝑓 cot𝜃 𝑎2 = 𝑎1 tan𝜃 Case 2 𝑎2 − 𝑎1 ∆𝑠 = tan𝜃 𝑎1 + 𝑎2 ∆𝑠 2 = ∆𝑠 2𝑓 𝑎1 = ∆𝑠 2 2𝑓 − tan𝜃 𝑎2 = 𝑎1 + ∆𝑠 tan𝜃 Case 3 ∆𝑠 − 𝑎1 ∆𝑠 − 𝑎2 = tan 𝜃 1 2 ∆𝑠 − 𝑎1 ∆𝑠 − 𝑎2 = ∆𝑠 2 1 − 𝑓 𝑎1 = ∆𝑠 1 − 2 1 − 𝑓 tan𝜃 𝑎2 = ∆𝑠 1 − 2 1 − 𝑓 cot𝜃 Case 4 𝑎1 = 𝑎2 𝑎1∆𝑠 = ∆𝑠 2𝑓 𝑎1 = ∆𝑠𝑓 𝑎2 = 𝑎1 Case 5 𝑎2 𝑎1 = tan 𝜃 1 2 𝑎1𝑎2 = ∆𝑠 2𝑓 𝑎1 = ∆𝑠 2𝑓 cot𝜃 𝑎2 = 𝑎1 tan𝜃 Case 6 ∆𝑠 𝑎2 − 𝑎1 = tan𝜃 𝑎1 + 𝑎2 ∆𝑠 2 = ∆𝑠 2𝑓 𝑎1 = ∆𝑠 2 2𝑓 + cot𝜃 𝑎2 = 𝑎1 − ∆𝑠 cot𝜃 Case 7 ∆𝑠 − 𝑎1 ∆𝑠 − 𝑎2 = tan 𝜃 1 2 ∆𝑠 − 𝑎1 ∆𝑠 − 𝑎2 = ∆𝑠 2 1 − 𝑓 𝑎1 = ∆𝑠 1 − 2 1 − 𝑓 tan𝜃 𝑎2 = ∆𝑠 1 − 2 1 − 𝑓 cot𝜃 Case 8 𝑎2 = 𝑎1 𝑎2∆𝑠 = ∆𝑠 2𝑓 𝑎2 = ∆𝑠𝑓 𝑎1 = 𝑎2 56 4.2 VOS for a Non-Uniform Mesh Most of the work done on the VOF method is based on a uniform Eulerian grid with rectangular elements. This is convenient due to the fact that the VOF method depends on calculating the amount of fluid moving through cell boundaries (fluxes) based on the fractional volume and the velocity field in each cell. Also, in most problems there is no need for a non-uniform mesh since the grid is fixed and the material is allowed to move freely within the domain. Nevertheless, some applications and problems need to use a non-uniform grid. Some researchers handled this problem through modified VOF techniques. These techniques are based on a mixed-two-step solution: a Lagrangian step followed by an Eulerian step. Although this method is efficient in most cases, it poses some challenges in problems where the mesh is too small in some regions. It also introduces more intense calculations since the solution consists of two steps rather than one. In addition to the difficulty in mapping parts of two or more Lagrangian elements into one Eulerian grid when needed [54-56]. In the following sections, we present a new technique to deal with non-uniform mesh and to calculate the fractional volume of each cell without the need for a Lagrangian step. The solution of the transport equation and the calculation of the free surface directional vector are similar to the ones presented for the uniform mesh whereas the construction of the free surface is where the two formulations differ. 4.2.1 Modification to the Solution of the Transport Equation Similar to the development for the uniform mesh, we calculate the fractional volume of each cell using Eq. (4.4) by generalizing it for quadrilateral element. We modify the calculations of the element volume and the outward unit normal vector to each edge by defining a reference node (O) in each element and choose this node to always be at the bottom left corner. Referring to Figure 4.6, a typical cell would be OBCD. The cell volume is simply calculated by summing up volumes of two triangles (OBC and CDO). Denoting a typical position vector between two grid points, e.g., OB, by R and its corresponding outward normal by n, we may write the following relation: 57 , for 0 and 0 , for 0 and 0 where , for 0 and 0 , for 0 and 0 x y x yy x x y x y R R R RR R R R R R n R R i j (4.12) The facial area is calculated for each face of the local cell using the length of the position vector between the two grid points and the cell thickness. O C D B ith element nOB nBC nDO nCD jth face Figure 4.6 3x3 non-uniform stencil with outward facial normals 4.2.2 Generalization of Free Surface Directional Vector We use the same general definition for the free surface direction vector and modify Eq. (4.6) to handle quadrilateral elements. In this case, each cell face will have both horizontal and vertical flux components and, therefore, we use the unit normal vectors for all the elements of the 3x3 stencil. The normal vectors will be calculated as inward vectors (Figure 4.7) for the element where the fractional volume f is considered. This will eliminate the problem of having to identify the sign for each component of the vectors 58 when we calculate the difference in the fractional volume values, and we simply add the values of the fractional volumes multiplied by its inward normal vector. To get a rational picture for the distribution of the fractional volumes, we multiply each fractional volume by its element size. The modified version of Eq. (4.6) is written as: 3 ,3 ,3 ,3 ,1 ,1 ,1 1 3 3, 3, 3, 1, 1, 1, 1 ,i i i R R R n n n n r k k k k k k k l l l l l l l i f V f V f V f V (4.13) where r is calculated for the i th element and f, V and n are the fractional volume value, the volume and the unit normal vectors of the element in the k th row and l th column in a 3x3 stencil, respectively. The main difference between Eq. (4.6) and Eq. (4.13) is that the numerator in Eq. (4.13) is no longer separated into horizontal and vertical components prior to the calculation, but is separated into rows and columns with cells having both horizontal and vertical components. Note that corner elements have two unit normal vectors, so for the calculation of the fluxes between neighboring columns we use the unit vectors normal to the left face of the elements in the right column and vice versa. We do the same thing for the flux calculation between neighboring rows. Figure 4.7 shows the element numbering, in circles, for a 3x3 stencil. In cases where one or more of the free surface elements are located at the end of the grid domain, we create a fictitious element beyond the domain with a reasonable f value that follows a linear progression similar to the scheme that was presented for the uniform mesh in Section 4.1.2. The constructed fictitious elements are considered to be uniform and Eq. (4.7) is used for the calculation of the fractional volume value of the fictitious elements. 59 nr 21 23 31 32 33 12 13 11 Figure 4.7 3x3 non-uniform stencil with unit vectors 4.2.3 Construction of the Free Surface using a Filling Pattern Scheme It is not a straightforward task to find a general set of equations that may locate the free surface based on the value of the fractional volume and the free surface directional vector for a non-uniform mesh. Therefore, classifying the free surface interface shape possibilities, as in the case for the uniform mesh, is not practical and a different approach is needed for the free surface construction. In this work, we present a filling technique that would incrementally fill the element using small strips through a specified pattern until it reaches the desired value of the fractional volume f that was calculated using Eq. (4.4). First, we choose a reference node for the start of the filling process. The reference node is chosen based on the signs of the free surface directional vector components defined in Eq. (4.13). Taking for example a typical element as in Figure 4.6 and Figure 4.8, we choose node O as the reference node if the both horizontal and vertical components of the directional vector (rx and ry) are 60 positive (+ve). Similarly, we choose node B for -ve rx and +ve ry; node C for -ve rx and - ve ry, and node D for +ve rx and -ve ry. Each free surface element is then divided into three sections that are defined by drawing perpendicular lines to the free surface directional vector from the two nodes neighboring the reference node (Figure 4.8). We start filling each section with a small area of thickness ds and width w. The thickness is defined by the user prior to the solution while the width is calculated at each time step based on the parameters defined in Figure 4.8. The total thickness of all the filled areas is defined as S. r Section 1 Section 2 Section 3 O C D B ds w r r α β O B’ D’ O B” D’ B γ β O D” r B” D B γδ S Figure 4.8 Free surface sectioning in the filling pattern scheme 61 All the angles shown in Figure 4.8 are calculated using the dot product between the free surface directional vector and the vector connecting the reference node with the point where the strip is added: 1 1 1 1 cos cos cos cos r R r R r R r R r R r R r R r R OB OD OB OD OB OD OB OD (4.14) Once the angles are known, the width of the small area may be calculated as follows: tan tan Section 1 tan tan Section 2 tan tan Section 3 S w S S (4.15) For each section, we multiply the width w by the thickness ds and divide by the cell volume to compute the fractional volume value of the added strip. Then, we compare the fractional volume of the filled area with the fractional volume value f. If the filled area is less than f then we keep filling the element until we reach the desired value of f. The total thickness of the filled area S is updated after each filling step by adding ds to the previous value. A bisectioning routine is employed if the filled area exceeds the desired value beyond the tolerance limit. In the routine, the thickness ds of the last filling step is reduced into a smaller thickness and the last filling step is re-initialized using the new ds value. 62 CHAPTER 5 SOLUTION CONTROL AND IMPLEMENTATION 5.1 Time Step Size and Stability The choice of solution time step size is crucial for stability and convergence when using the VOS method. Based on the analysis presented in this work, the time increment must not allow the material to flow through more than one element in any direction. This condition is imposed for both cases of positive or negative net change for each element to ensure that an element will not be filled/emptied more that its volume. This is a very important stability concern since the VOS equations presented in this work are based on the assumption that fluxes occur only between adjacent elements. We calculate the critical time step size for all the elements at each time step, and compare it with the minimum critical value given by the following expression: 1 0 min , j j i i cr j j j j f V f V t f A f A v n v n (5.1) Equation (5.1) has two different time steps, one of them computes the critical time step for a cell with a positive net change in the material volume and the other one is for a cell with a negative net change in the material volume. Note that the denominator in Eq. (5.1) represents the material net change in an element. This ensures that no element is filled with material more than its volume, and no filled (or partially filled) element is emptied more than its volume of material in one time increment. 63 5.2 Updating the Velocity Field As the material flows within the fixed Eulerian grid, it is bound to occupy previously empty elements. These elements were not part of the computational domain in the previous time step but now must be included due to the recent presence of the material. The newly filled elements have no velocity values at their nodes and they need to be assigned a reasonable value that matches the velocity field distribution around them in order to be able to continue the simulation and advance the free surface. On the other hand, elements that are being emptied will be assigned a zero velocity at their nodes. In this work, we use a volume weighted averaging process [74] for the velocity calculation at the new nodes by averaging the velocity field of the neighboring elements: new node v v i i i V V (5.2) where 𝐯 𝑖 is the average velocity of the neighboring elements calculated by averaging the velocity of all nodes in that element and Vi is the volume of the neighboring elements. Figure 5.1 shows the newly added nodes at a free surface interface between times t and t+Δt. Current computational nodes Newly added nodes Time t+ΔtTime t Figure 5.1 Velocity field update for newly added nodes 64 5.3 Tracking Material Point Properties One of the main concerns in Eulerian and ALE formulations is how to deal with history dependant properties such as the stress and the strain. Because these physical parameters represent the history of the material during the course of the deformation, finding a proper formulation for them is essential. To do this, we apply Eq. (3.3) relating the grid and material time derivatives as follows; written for the equivalent strain as a typical quantity: t t t t i t i v x (5.3) The solution of Eq. (5.3) can be obtained globally by applying the Petrov-Galerkin approach and using the element shape functions of the Eulerian grid: 0 t t t t i i it iV V v h dV h dV x (5.4) where vi is the material velocity, hi are the shape functions and V is the volume of the elements. The equivalent strain rate 𝜀 can be obtained from the velocity calculations. The solution of Eq. (5.3) could be achieved in a simpler manner at the free surface of the material if the trajectories are known. In Sections 4.1.2 and 4.2.2, the free surface location was calculated by finding the free surface directional vector which is a representation of the trajectory of the material points within the free surface element. Therefore, if the value of the equivalent strain is known at a certain material point M 1 , then it can be calculated progressively along the directional vector following the same material point on the new updated free surface as shown in Figure 5.2. This is achieved using the following equation [75]: gradt t t ti i i i ix t x x x x (5.5) 65 S 1 S 2 S 3 M 1 M 2 M 3 S 1 S 2 S 3 Figure 5.2 Tracking material points at the boundary using the free surface directional vectors 5.4 Free Surface Connectivity Even though the VOS method might ensure the connectivity of the free surface interfaces since they are calculated based on the fractional volume values of neighboring elements, it is not always guaranteed to have the surface connected between all elements at the free boundary. The objective here is to ensure connectivity of the free surface since such condition is not guaranteed. Therefore, we perform a connectivity step after each free surface capturing step is completed. To guarantee the free surface to be connected at all times, we implement a default averaging scheme that connects the surface of two neighboring elements that share the same face. This is not the case for elements sharing a face with fictitious elements because their free surface remains the same from that side neighboring the fictitious element and no connectivity is required at that point. Once a connectivity step is performed, all subsequent surface interfaces must account for the change in neighboring free surface elements. Based on the definition of facial fractional value in Eq. (4.5) and Figure 5.3, the default averaging process is performed as follows: 66 2 if 0 1 and 0 1 1 if 1 or 1 if is fictiticious if is fictiticious j j j j j j j j j j j f f f f f f f f f f f (5.6) where the + and - signs correspond to the opposite sides of the facial fractional volumes that share the same element side. If one facial value is equal to 1.0, we assume both values to be 1.0 so that a full element will not turn into a partially occupied one, and if one of the sides is a fictitious element we keep the value the same as it is for the true element. The process starts from one element and is repeated for all other surface elements to average all the values and their after effect. The reason is because after connecting the free surface of two elements, the free surface of elements neighboring those two must account for the effect of the connectivity of those two elements to insure that the total volume of the solid is not changed during the simulation, i.e. to ensure that we do not lose any mass during the approximation (Figure 5.3). Free surface nodes before connectivity Free surface nodes after connectivity Figure 5.3 Default free surface connectivity scheme 67 Originally, in fluid mechanics, the VOF method has been used in applications where the shape of the fluid free boundary is not exactly known in many cases and the VOF method was then used to predict the shape of the free boundaries. On the contrary to CFD applications, the expected final deformed shape of the workpiece in metal forming applications is known or at least partially known. Taking advantage of the pre-known deformed shape of the material in some regions can improve the quality and efficiency of the results obtained using the VOS method. The first important modification we do to the default connectivity process is employed in the contact between the workpiece and rigid punch/dies. The connectivity process is modified so that it would not compromise free surface elements that are in contact with the punch/dies. The deformed shape of elements that are in contact with rigid objects will not be changed and will dominate the connectivity between them and neighboring elements that are outside of the contact region. Figure 5.4 shows the effect of having the contact-modified connectivity for a simple punch indentation problem compared to the default connectivity process. We know that the punch land area must remain flat due to the contact with the punch and not inclined, which is what happens with the default connectivity process as in Figure 5.4c. It is also clear that the default connectivity process will introduce some penetration between the punch and the material whereas the contact-modified process eliminates any penetration after effect (Figure 5.4d). The second important modification to the default connectivity process is by introducing a weighting factor, with a value between zero and one, to the default averaging process in Eq. (5.6) to change the connectivity between two adjacent surfaces from a 50% averaging process to any other desired value. This can be used to modify the shape of curves and decide if a steep change in the curve is required towards the center or a gradual change is more suitable. Figure 5.5 shows the difference between two curves one with default connectivity and the other one with higher weighting factor value at both end of the curve. The weighting factor can also be used in other cases such as dealing with sharp corners of a workpiece during the deformation so that the corners will not lose their identity through the connectivity process. 68 Flat punch (a) before contact (b) after contact (c) default connectivity (d) contact-modified connectivity Figure 5.4 Contact-modified free surface connectivity (a) Default connectivity (b) Weighting factor controlled connectivity Figure 5.5 Effect of weighting factor on the free surface connectivity process 69 5.5 Contact Analysis One of the main features when simulating metal forming processes is modeling the contact between the workpiece and the tool. This is a very complex task due to the need for tracking the movement of both bodies, defining the state in which contact occurs between them, and applying constraints and boundary conditions that simulate the contact. A common assumption in metal forming processes is to model the contact set as a rigid-flexible pair, where the tool is modeled as a rigid body and the workpiece as a flexible body. In this analysis, we use the rigid-flexible contact assumption since flexible- flexible contact is beyond the scope of this research. Usually modeling the punch or the dies as a flexible entity is essential where the elastic deformation of the tool is of main interest. In this work, the focus is on the large deformations of the workpiece. The two most common techniques used in the literature for the rigid-deformable contact are the Lagrange multiplier and the penalty methods. In the Lagrange multiplier method, a constrained minimization problem is solved where the constraint does not allow any penetration between the two contact bodies [28]. One of the problems with this method is that it increases the computational cost of the analysis because of the inclusion of the Lagrange multipliers into the system. Also, because special elements (gap elements) are used in this method, knowledge of where the contact will occur and a specification of the whole contact surface prior to the simulation are required [25]. The penalty is easier to employ but it allows some penetration between the contact bodies based on the penalty constant which also affects the accuracy and the stability of the numerical solution [28, 76]. In this work, a direct approach is used based on contact module presented for ALE analysis in [25, 28]. The contact module was used to track the contact state between a moving grid and a rigid tool, therefore it needs to be modified, so that it can track the movement of the material since the grid is fixed in space, and incorporated into the Eulerian formulation. This direct approach is originally based on a node-to-surface contact algorithm in which the tool is designed as geometrical segments, such as lines and arcs, which are connected together and dealt with as a whole surface. Since elements are fixed in an Eulerian grid, modifications are made to this contact module by adjusting 70 the way rigid segments are defined and changing the contact detection tracking scheme. The tool segments are now defined on the elemental level, which means that each tool segment is defined within an element by itself instead of defining all tool segments as one body. This allows the module to track the movement of the toll segments on the Eulerian grid. The contact is now detected at the free surface elements between rigid segments and free surface interfaces, and is now defined as a node-to-line which means the contact will be dealt with in each element separately rather than in the structure as a whole. In each incremental time step, the algorithm first tracks the movement of the tool segments and records the coordinates of its final position. Then it tracks the movement of the boundary nodes, which are located at the intersection between the free surface interface and the element edges. The module then compares the coordinates of the tool and the free boundaries to detect if any penetration has occurred between them. Once penetration is detected, the algorithm calculates the depth and the direction of the penetration between the free surface interface and the rigid segments based on the free surface location that was calculated using the VOS method. If penetration occurs at more than one surface element, the algorithm determines the element at which penetration occurred first and starts with it first and then updates other contact elements accordingly before checking for penetration again. Once the penetration depth is calculated geometrically, a corrective velocity boundary condition, based on the value of the penetration depth, is enforced on the tool to reset the tool surface on the free surface interface. The corrective velocity should compensate for the free surface movement as the workpiece deforms during the process. Taking that into account, the corrective velocity constraint is calculated as follows (Figure 5.6): K I I J J pv N v N v v (5.7) where N I and N J are the shape functions corresponding to nodes I and J, respectively, and v p is the penetration vector that is calculated geometrically based on the depth of the contact penetration. Once the tool segments are repositioned, appropriate contact boundary conditions are applied to the boundary nodes that are in contact. The friction force of the contact is calculated based on Coulomb or shear friction laws, depending on 71 the problem, and applied to the contact nodes. The contributions of the contact forces are added to the Eulerian equilibrium equation: t t t t C t t t t t t ext C t t C t t C ij t t ij i i V S e dV W u f dS (5.8) where 𝛿𝑢𝑖 𝐶 are the components of the virtual displacements on the contact surface defined at the boundary nodes and 𝑓𝑖 𝐶 are the components of the contact traction acting on the contact surface S C . For each boundary node, the contact algorithm will be performed for each surface element that includes that node separately. At the end, all the resulting forces and contact boundary conditions that were calculated for each element sharing that node are combined together to the shared node to simulate the effect of the contact. The algorithm also checks if the nodal force, normal to the contact surface, is compressive or tensile. If the force is compressive then the node is still considered to be in contact. However, if it is in tensile mode then the node is considered to be separated from the contact surface and these forces are set to zero. Figure 5.6 shows the penetration detection and the corrective repositioning steps. I J K v p Figure 5.6 Penetration detection and surface repositioning in contact analysis 72 5.6 Heat Transfer Analysis The internal heat generation in metal cutting process plays an important role and cannot be neglected. Therefore we present in this section a heat transfer module that is used in the VOS simulation of metal cutting in Chapter 8. Most of the work dissipated due to the plastic deformation and friction is converted into heat that will raise the temperature of the workpiece and tool significantly. In this work we use an incremental staggered approach in which we solve a deformation problem followed by heat transfer problem in each increment. The deformation problem is performed as an isothermal stress analysis one with known temperature distribution at the beginning of each incremental step. The work dissipated in the stress analysis step is used as the source for the internal heat generation in the heat transfer analysis part at the end of that same incremental step. This approach is considered acceptable for small time steps and when the dependency of the material properties on the temperature is weak [25]. The general form of energy equation for fixed set of coordinates is written as: , ,i i icT cv T kT q (5.9) in which v is the velocity, T is the temperature, ρ is the density, c is the specific heat, k is the conductivity coefficient and 𝑞 is the rate of heat generation per unit volume. The boundary conditions for this equation are: 1 2 3 1 2 1 i i s m mS S mS q n q h T T T T (5.10) which represents specified flux, convection and specified temperature on respective parts of boundary, S1 to S3. In Eq. (5.10), h represents the convection heat transfer or film coefficient, T is the temperature and n is the surface normal vector. The internal heat generation is defined by calculating the work done in the deformation step as follow: iij ij ij j v q D D x (5.11) 73 where 𝜎 and 𝐷 represents the effective stress and effective strain rate, respectively. Practically, not all of the dissipated work is converted into heat and only the dissipated plastic work is counted for [77, 78]. Therefore we rewrite the internal heat generation as: pq D (5.12) where 𝐷 𝑝 is the effective plastic strain rate and η is a factor that represents the amount of the plastic work converted to heat which is found to range between 85-95%. Similarly, using the same definition for η, the heat generation due to friction between the workpiece and the tool can be written as: S f fq v (5.13) where 𝜎𝑓 and vf are the frictional stress and the tangential velocity at the contact surface. 5.7 Solution Procedure The VOS formulation presented in this work is implemented into a finite element code. The code is divided into three major steps. First we need to define the parameters of the process which includes all parts of the FE pre-processor stage. Next we set up the solution control parameters and apply connectivity modification where needed. After the problem setup is completed, the simulation is started by solving the Eulerian equilibrium equation and all related VOS equations. The solution procedure is summarized as follow (Figure 5.7): (1) Problem Definition: Specify material properties Generate mesh domain and specify initial distribution of f Define tool geometry Apply boundary conditions for the workpiece and the tool 74 (2) Solution Setup: Create contact pairing between the workpiece and the tool(s) Define contact parameters Apply connectivity weighting factor if needed Specify solution time step (3) Solution Routine: Solve Eulerian equilibrium equation Obtain displacement and velocity fields Solve the transport equation and obtain the new distribution of the fractional volume values f Calculate the free surface directional vector Construct the free surface interface using the geometrical PLIC scheme (uniform mesh) or the filling pattern scheme (non-uniform mesh) Apply default connectivity process with specified weighting factors Apply contact-modified connectivity Update the velocity field for the newly added nodes at the free surface Obtain material point related properties (4) Repeat the solution routine (3) until final solution time is reached. 75 - Material properties - Mesh - Tool geometry - Boundary and initial conditions - Contact definition - Contact parameters - Connectivity weighting factor - Initial time step - Solve equilibrium equation - Obtain displ. and vel. fields - Solve transport equation - Obtain distribution of f - Calculate free surface directional vector - Construct free surface interface - Apply connectivity process - Update vel. Field for newly added nodes - Obtain material point related properties End of simulation? No Stop Yes Problem Definition Solution setup Solution Figure 5.7 Flowchart of the VOS FE-code solution procedure 76 CHAPTER 6 METAL FORMING APPLICATIONS 6.1 Punch Indentation Problem The punch indentation problem is simulated to show the advantages of using the VOS method compared to the traditional Lagrangian formulation used in commercial FE software, e.g., ANSYS [76]. The main issues involved in simulating such processes are the mesh distortion due to large deformations, the penetration between the punch and the workpiece during the process, and the convergence problems, especially in achieving large penetration ratios. The workpiece is compressed between two frictionless rigid dies moving with constant velocity until 50% height reduction is achieved. Due to symmetry, only one quarter of the problem will be modeled using plane strain conditions. Figure 6.1 shows the setup of the problem. The material model of the workpiece is assumed to be a bilinear elastic-plasticwithYoung‟smodulusofelasticityE =200GPa,Poisson‟sratiov = 0.3, initial yield stress Sy = 250 MPa and tangent modulus ET = 1000 MPa. First, the problem is solved with ANSYS using Lagrangian formulation with 140 equally sized elements. For the VOS simulation, the whole domain is meshed into 600 elements with the original workpiece filling initially 300 elements. A simple rigid- flexible contact algorithm is used for this problem which detects the contact at the free surface or free boundary of the material. Contact-modified connectivity is used to keep elements under the punch land flat during the process. Figure 6.2 shows the evolution of the workpiece deformed shape from the VOS simulation. We notice that there is no penetration between the punch and the workpiece especially at the punch corner since the contact area between the punch and the workpiece is kept flat during the entire process with the help of the contact-modified connectivity process. 77 Figure 6.3 and Figure 6.4 show the von-Mises strain and stress distributions in ANSYS and the VOS model, respectively. The distribution and values of the strain and stress in the VOS simulation are similar to those obtained with ANSYS, however, the advantage of the VOS method is that the distribution of the values is smoother and does not possess high jumps. Also, by comparing the final deformed shape between the two models we can see how the punch land is completely flat in the VOS simulation where on the other hand slight penetration occurs between the punch corner and the workpiece in the ANSYS model. The main purpose behind the simulation of this application was to validate the VOS FE-code by comparing with well-established commercial FE softwares such as ANSYS using a metal forming process that can be simulated using Lagrangian formulations without the need for any special treatment or re-meshing techniques. The results showed strong agreement between the VOS code and ANSYS. The VOS method is used to deal with the free surface tracking only for the elements at the boundary so that we do not have to deal with the mesh motion of the workpiece. Also, all of the empty cells are discarded from the solution until they are filled by the material, unlike traditional VOF methods where empty cells are considered to be filled by air and are included in the computational domain. Finally, by comparing the setup of the problem and the steps required at the preprocessor stage we see that in the VOS simulation no special contact options are required to be specified by the user. Also, in VOS there are no contact-tuning parameters to handle, such as the allowable penetration or the normal penalty stiffness parameters, which may need some tweaking when using ANSYS to get the optimum results. 78 20 mm 24 mm 60 mm Punch Workpiece Figure 6.1 Geometrical setup for the punch indentation problem (a) Original undeformed shape (b) 20% height reduction (c) 40% height reduction (d) Final deformed shape Figure 6.2 Evolution of the workpiece deformed shape in the punch indentation problem using VOS 79 (a) ANSYS (b) VOS Punch Figure 6.3 von-Mises plastic strain distribution for the punch indentation problem 80 (a) ANSYS (b) VOS Punch Figure 6.4 von-Mises stress distribution for the punch indentation problem 81 6.2 Compression between Wedge-Shaped Dies Plane strain compression of a pre-shaped workpiece between two wedge-shaped dies is simulated using the VOS method. The same process was investigated experimentally by Chitkara and Johnson [79] and numerically by Gadala and Wang using an ALE formulation [22]. The purpose of this example is to demonstrate the capability of the VOS method in dealing with metal forming applications and the importance of using a non-uniform grid to setup the FE model of the problem. Since the workpiece is pre- shaped at the beginning of the process, then using a uniform mesh is not possible and a hybrid grid of uniform and non-uniform elements is essential to model the problem. The process consists of two rough wedge-shaped dies pressing on two opposite sides of a pre-shaped workpiece. The geometrical setup of the problem is shown in Figure 6.5 with H defining the mean distance between the centers of both contact lines with a value of 24.07mm. Only half of the problem is modeled due to the symmetry. The workpiece is made of half-hard aluminum and assumed to have bilinear elastic-plastic material model withYoung‟smodulusofelasticityE =68.9GPa,Poisson‟srationv = 0.33, yield stress Sy = 276 MPa and a tangent modulus ET = 25.07 MPa. The process is simulated up to 45% reduction in the initial height using a prescribed displacement. The experiment in [79] was conducted on a 300-ton Denison hydraulic testing machine with a load less than 3.45 MPa. A 121.6 x 15 mm window is used for the Eulerian grid in the VOS simulation with 778 elements, 436 of these elements are initially filled by the workpiece in its original configuration. Since the dies are rough, sticking contact condition is applied between the die and the workpiece. The deformed shape of the workpiece is shown in Figure 6.6 for 15%, 30% and 45% reduction in height. The deformed shapes of the workpiece are very similar to the ones presented using ALE by Gadala and Wang [22]. We can see that by using the VOS method we avoided the mesh distortion problem that occurs in this type of applications when using traditional Lagrangian formulations. The plastic strain distribution during the deformation stages is shown in Figure 6.7. Comparing the values presented in here with the values obtained using ALE [22] we find that the strain values are very close to the 82 ones predicted using ALE. From the distribution of the strain at different stages, we can see that the plastic zone is formed around the die corners where the strain values are the highest. The central part of the interface has lower strain values, which is contributed to the sticking contact condition. The plastic zone starts forming under the corners and develops towards the center of the workpiece. No plastic zone occurs outside the dies and the strain values are very small. We conclude that the strain distribution is in agreement with the predicted behavior of such processes and is also similar to the behavior predicted by the ALE formulation. By defining h to be the mean height of the workpiece during the deformation, we identify the compressive strain as (ln H/h). Figure 6.8 shows the load variation on the die with respect to the compressive strain obtained from the VOS, ALE and experiment results. The load-strain curve from the VOS simulation matches very well with the curve obtained using the ALE formulation. Both, VOS and ALE, curves are very close to the experimental curve for strain values less than 0.35-0.4. A noticeable difference can be seen between the experimental and the VOS/ALE curves for higher strain values and this is due to the fact that a reduction in friction occurs in the experiment as the process continues, were in the numerical analysis the contact friction stays the same during the whole process. 25.4 mm 22.74 mm 38.1 mm 38.1 mm25.4 mm H Figure 6.5 Setup for the compression between wedge-shaped dies problem 83 (a) Undeformed shape (b) 15% Height reduction (c) 30% Height reduction (d) 45% Height reduction Figure 6.6 Evolution of the workpiece deformed shape in the compression between wedge-shaped dies problem using VOS 84 (a) 15% Height reduction (b) 30% Height reduction (c) 45% Height reduction Figure 6.7 von-Mises strain distribution in the compression between wedge-shaped dies problem 85 Figure 6.8 Load-strain comparison curves for the compression between wedge- shaped dies problem 0 4 8 12 16 0 0.1 0.2 0.3 0.4 0.5 0.6 P ( M P a) ln H/h VOS EXP ALE 86 6.3 Sheet Metal Extrusion Sheet metal extrusion (forward extrusion) is one of the metal forming processes where large strains are expected due to the change of the sheet thickness. When using traditional Lagrangian formulations, remeshing is required to overcome the problem of the mesh entanglement at the region of the thickness reduction. Updating the boundary conditions is also important because the outer boundary of the workpiece is in contact with the die during the whole process. In this example, we use the VOS method to demonstrate the ease of simulating such applications without remeshing and updating the boundary conditions since we are using an Eulerian grid in which all boundary conditions are pre-defined for each region of the problem for the whole process. In this application, a plane strain sheet metal extrusions is simulated for a 25% reduction in the sheet thickness. The die is shaped in the form of 5 th order polynomial curve with zero curvature and slope at both ends. An aluminum billet with a thickness of 2a is forced between two frictionless rigid dies by a rigid piston pressing against the rear face of the billet. The piston is given a prescribed displacement of 4a. The length of the reduction die is 1.2a. Figure 6.9 shows the geometry of the problem at the undeformed stage. The material model of the aluminum billet is assumed to be bilinear elastic-plastic withYoung‟smodulusofelasticityE = 6.8947x104 MPa,Poisson‟srationv = 0.3, yield stress Sy = 393 MPa and a tangent modulus ET = 1137 MPa. Due to the symmetry, only one half of the process is considered. The same problem was solved by Gadala and Wang [20], and Bayoumi and Gadala [68] using ALE formulations with the same material model, and by Maniatty et al. [7] using Eulerian formulation and a different material model. The difference between the listed references and this work is that they simulated the problem as steady state with part of the billet already being extruded. In our case we start from the undeformed state of the billet to show the capability of the VOS method to handle to process from the beginning. An 11.5a x 1.5a window is meshed into 1202 elements for the Eulerian domain with 352 elements being initially filled by the material of the undeformed shape of the billet. In this application we take advantage of using a combination between uniform and non- uniform elements to mesh the domain. By doing that we can improve the shape of the 87 deformed billet under the die since it is shaped as a 5 th order polynomial. Two contact sets were defined in this problem, the first one is between the fixed rigid die and the billet and the other one is between the rigid piston and the billet. Both contact sets are defined with no additional contact treatment. The deformed shape of the billet is controlled by the fixed rigid die during the whole simulation with the help of the non-uniform mesh under the die that fits the shape of the billet once it reaches the steady state stage. Figure 6.10 shows the FE model of the problem in addition to the deformed shape when the billet enters the region of the extrusion die and the final deformed shape. It is clear from Figure 6.10b that the irregularity in the mesh under the reduction die helps the billet material to flow according to the shape of the 5 th order polynomial accurately without any compromisation. The distribution of the von-Mises strains after a displacement of 4a is shown in Figure 6.11. The distribution and the values of von-Mises strain are in a very good agreement with the values obtained using ALE in [20, 68]. The values of the longitudinal stress component at different lateral positions from the mid- plane of the initial configuration are shown in Figure 6.12 starting from 1.5a before the reduction die up to 2a after the exit of the die. Generally, the values and the fluctuation behavior between compression and tension through the three regions of the problem are in agreement with the same figure that was presented in [20, 68]. The differences between the results are due to fact that in this work we simulated the process from the initial undeformed configuration which contributed to the slight difference in values at the region before the reduction die. Maniatty et al. [7] used an Eulerian formulation method for simulating steady state forming processes. In this work we used the VOS method with Eulerian formulation to simulate the process from the beginning rather than simulating only the steady state part. Maniatty et al. [7] showed in their work that the ratio between the length of the reduction die and the mean thickness determines the pattern of the residual stresses. For low ratios, they showed that residual stresses should be compressive towards the center of the billet, whereas on the surface they tend to be in tension mode. In this problem we have a ratio of approximately 1.37 which is considered low based on their work, so we expect the same kind of behavior for the longitudinal stresses in our case. For b = 0.75a, 88 Figure 6.13 shows the longitudinal stress along the thickness (y) of the billet at a distance 4a from the die exit. We can see that the pattern is as expected. Stresses are in compression towards the center of the billet and in tension towards the surface. The above example shows the advantage of using a non-uniform mesh in obtaining an accurate shape of a deformed metal sheet. Compared to other methods, there is no need to jump to the steady state stage and ignore the initiation of the process in order to simplify the problem. The process is simulated from the beginning which allowed us to get better results that gave a clearer picture of the process and its behavior. Figure 6.9 Geometrical setup for the sheet metal extrusion problem 89 (a) Initial configuration (b) After 0.8a stroke (c) Final Configuration Figure 6.10 Evolution of the billet deformed shape in the sheet metal extrusion problem using VOS Figure 6.11 von-Mises strain distribution for the sheet metal extrusion problem 90 -1500 -1000 -500 0 500 1000 1500 Lo n gi tu d in al S tr es s (M P a) 0.125a 0.375a 0.625a 0.875a 1.5a 1.2a 2a Figure 6.12 Longitudinal stresses at different lateral positions for the sheet metal extrusion problem Figure 6.13 Longitudinal stresses along the thickness of the billet in the sheet metal extrusion problem -1000 -750 -500 -250 0 250 500 750 1000 0 0.2 0.4 0.6 0.8 1 R e si d u al S tr e ss ( M P a) y/b 91 6.4 Backward Extrusion In backward extrusions, mesh distortion is a critical problem and it is not possible to use traditional Lagrangian formulation due to excessive distortion and folding of the elements unless an extensive re-meshing algorithm is used. Using the coupled Eulerian and VOS formulation eliminates the problem of the mesh distortion and easily overcomes the problems of contact and free boundary tracking. The following example shows that such processes can be smoothly simulated using an Eulerian formulation with the VOS method and a uniform mesh. This example is based on the same setup presented by Hur et al. [80] and the results from the VOS simulation will be validated by comparing them to the results presented by [80]. The geometrical setup is shown in Figure 6.14. The only difference in this model is that a rigid die is used instead of an elastic die, and a bilinear elastic-plastic material model is used. The workpiece is made from AISI-1010 steel with a flow stress relation: 𝜎 = 715.95𝜀 0.22. The bilinear material modelhasYoung‟smodulusE = 206754 MPa, Poisson‟sratiov = 0.3, yield stress Sy = 305 MPa and tangent modulus ET = 780 MPa. In the simulation, the friction coefficient between the punch and the die with the workpiece is 0.1. The problem is simulated as an axisymmetric problem and the simulation is performed for a total punch stroke of 10 mm. A total of 2,556 elements were used to mesh the entire domain, with 1,353 elements used for the initial configuration of the workpiece. The mesh used in the domain has a higher intensity in the bottom half of the domain to increase the accuracy of the curved workpiece shape around the contact area. Elements that are in contact with the punch land are defined to take the shape of the punch during the whole simulation with the help of the contact-modified connectivity process. The bottom edge of the workpiece is fixed in the vertical direction, but left free in the horizontal direction. The process is solved as an axisymmetric FE model by using only the right half of the model and applying axisymmetric boundary conditions on the left edge of the model. Figure 6.15 shows the original undeformed shape of the workpiece and the deformed shape at 2, 4, 6, 8 and 10 mm punch strokes. Note that at the region where the punch 92 contacts the workpiece, the material flows smoothly around that corner and goes into the groove between the punch and the die without any complications. Also, the curved part of the punch land is accurately represented using the contact algorithm in our VOS code in addition to the weighting factor that was used in the connectivity procedure. The load-stroke diagram for the backward extrusion is shown in Figure 6.16. The values are comparable to the ones presented in [80]. The differences between the two curves could be attributed mainly to the difference in the FE formulation of the two methods in addition to the approximation of the material model used in the VOS simulation. Usually a flow stress relation, like the one used in [80], has a curved shape stress-strain relation where in our model we used an elastic-plastic model that has one slope for the elastic region and one for the plastic region. This difference in the material model causes the results to differ at some regions, and come closer at others. The contact pressure at the initial contact region between the workpiece and the rigid die is shown in Figure 6.17 for different punch strokes. Compared to [80], the contact pressure in both models starts increasing rapidly as we progress, until we reach the last stages of the process where it starts decreasing gradually. We also notice that the region with nonzero contact pressure starts decreasing as we progress through the stages of the process. The reason behind this is that the material starts moving more smoothly as the material of the workpiece enters the plastic region. With the difference in the material model and the fact that a rigid die is used rather than an elastic one as in [80], the maximum values are within an acceptable range. The development of the von-Mises strain distribution is shown in Figure 6.18 for various punch strokes. The development of the strains is very logical throughout the simulation since the strains start to build up around the punch land corner where the material experiences high rotations and large deformation. The bottom right side of the workpiece, where it meets the die, has the lowest strain values during the whole process and this is because the material gets pressed against the die in that region and stays at the same state during the entire process. This is also one of the reasons why the contact pressure values are high at that region. Compared to the strain distribution at the final stage (punch stroke of 10mm) presented in [80], the values in our work are in very good 93 agreement, and the distribution is quite similar in the two cases. Table 6.1 compares the values of the maximum pressure at different punch strokes and the strain values at the final stage. Overall, values presented using VOS are very close to those in [80]. The VOS method was used in this application to demonstrate its ability in handling nonlinear problems with large deformation and rotations and. Backward extrusion is considered to be among the metal forming application that are difficult to simulate by means of regular finite element formulations. The process was simulated efficiently with no need for re-meshing techniques or any mixed formulation. The combination of the Eulerian formulation and the VOS method presented the expected deformed shape of the workpiece in such applications. The results obtained from the VOS method were very good compared to results from the literature, and the distribution of these values follows the common knowledge that is known from this application. Workpiece Punch Die 11 mm 1 5 m m 11.4 mm 7.75 mm 4 mm R1.5 7 o Figure 6.14 Geometrical setup for the backward extrusion problem 94 (a) Undeformed shape (b) 2mm punch stroke (c) 4mm punch stroke (d) 6mm punch stroke (e) 8mm punch stroke (f) 10mm punch stroke Figure 6.15 Evolution of the workpiece deformed shape in the backward extrusion problem using VOS 95 Figure 6.16 Load-Stroke curve for the backward extrusion problem 0 648.86 0.9 mm 0 891.33 3.4 mm 0 1051.9 7.3 mm 0 1066.5 8.6 mm 0 1054.8 10 mm Figure 6.17 Contact pressure (MPa) between workpiece and die at different strokes in the backward extrusion problem 0 50 100 150 200 250 300 350 400 0 2 4 6 8 10 12 P u n ch L o ad ( kN ) Punch Stroke (mm) [80] VOS 96 (a) 2mm punch stroke (b) 4mm punch stroke (c) 6mm punch stroke (d) 8mm punch stroke (e) 10mm punch stroke Figure 6.18 von-Mises strain distribution for the backward extrusion problem 97 Table 6.1 Contact pressure and strain comparison between VOS and [80] Value Method Max. contact pressure (0.9mm) Max. contact pressure (3.4mm) Max. contact pressure (7.3mm) Max. contact pressure (8.6mm) Max. contact pressure (10mm) Max. effective strain value VOS 648.86 891.33 1051.9 1066.5 1054.8 2.5 Hue et al. [80] 858.28 967.255 1037.95 1027.11 998.363 2.4 98 6.5 Local Loading of a Simply-Supported Plate The plate indentation in this example is simulated using two different indenters; a flat-faced and hemispherically-tipped indenters. The process consists of a circular steel plate that is loaded locally while the plate is resting on a ring support. The VOS method is used in this example to show the capability of dealing with sharp corner contacts (flat indenter) and curved shape contacts (hemispherical indenter) without the need to change the contact options or parameters between the two cases. Avoiding mesh distortion while accurately representing the contact face without any penetration is accomplished easily using the VOS with the Eulerian formulations. The two cases are compared with experimental results [81] and numerical results using ALE formulation [28]. The plate is 6 mm thick, 300 mm in diameter and is simply supported at a diameter of 248 mm. Both indenters are 12.7 mm in diameter and are loading the plate at a rate of 3 mm/min. The plate is made of mild steel and is assumed to be bilinear elastic-plastic withYoung‟smodulusofelasticity E =235GPa,Poisson‟srationv = 0.3, yield stress Sy = 305 MPa and a tangent modulus ET = 1000 MPa. The process was simulated for a total displacement of 20 mm for the flat-faced indenter and 30 mm for the hemispherically- tipped indenter. Figure 6.19 shows the setup with the flat and hemispherical indenters. An axisymmetric finite element model was used for this process. A 150 x 30 mm window, meshed into 1191 elements, is used for the flat indenter case whereas a 150 x 40 mm window is used for the hemispherical indenter with 1566 elements. Both Eulerian grids had 262 elements initially filled with the initial configuration of the steel plate. The contact is defined between the rigid indenters and the plate with the contact-modified connectivity process that ensures an accurate shape of the indenters land shape. The contact is assumed to have no friction and no separation between the indenter and the plate. The mesh is refined at the top right corner of the grid domain to have a better prediction for the shape for the sharp corner of the steel plate since the corner is connecting two perpendicular lines and the VOS method uses a single free surface line at each boundary element. The smaller elements allow the sharp corner to progress through the simulation while keeping the geometry intact. 99 Figure 6.20 shows the final deformed shape of the finite element model for the circular steel plate when loaded using the flat-faced and the hemispherically-tipped indenters, respectively. Compared with the deformed shapes presented in [28], the deformed shapes using the VOS method are in a very good agreement with those obtained using ALE. The load-displacement curves obtained from the VOS method, ALE [28] and the experiment [81] are shown in Figure 6.21 and Figure 6.22 for the two cases. The experimental results in [81] showed that the plate fails by perforation at approximately 20 mm under the flat-faced indenter and 30 mm under the hemispherically-tipped indenter which is why the simulation was performed for these two values. From the two figures it is clear that the results obtained from the VOS method are in a good agreement with the experimental and ALE results. Even compared with the ALE results, the results from the VOS method are closer to the experimental ones. Generally we can see that at higher deformation levels, the numerical curves starts to deviate from the experimental one, and that is due to the discrepancies between the actual material model and the bilinear elastic- plastic model used in the numerical simulations. Figure 6.23 compares the load-displacement curves between the flat and the hemispherical indenters. The flat indenter produces a slightly steeper curve than the hemispherical indenter, but with a lower failure value. From the figure we can see clearly that the indenter profile affects the indentation process to a great deal. The hemispherical indenter causes a higher degree of local indentation and more bending around the contact region than the flat indenter which the sheering effect to decrease allowing the plate to withstand higher load values before perforation occurs. Figure 6.24 and Figure 6.25 show the deformed shape and the von-Mises strain distribution around the contact region. The deformed shapes and the strain distribution values are very similar to those obtained using ALE [28] for the hemispherically-tipped indenter, whereas the flat-faced indenter has some differences in the strain distribution around the indenter corner. The reason is because we used a rigid indenter and modeled the contact between the indenter and the plate, where in ALE results the indenter was substituted with an equivalent boundary condition to simulate its effect. In the rest of the 100 region the values and the distribution were very similar. We can observe the indenter land is flat in the first case which shows that the VOS method handled the contact with no penetration. In the second case, the effect of the curved shape indenter was presented accurately with the VOS method. The figures show that the hemispherically-tipped indenter caused a greater degree of local indentation and bending around the indenter than the flat-faced indenter. The VOS method eliminated the need to deal with the mesh distortion under the hemispherically-tipped indenter. 248 mm 300 mm 6 mm 12.7 mm Flat Indenter Hemispherical indenter Figure 6.19 Geometrical setup for the plate indentation problem 101 (a) Flat Indenter (b) Hemispherical Indenter Figure 6.20 Deformed shape of the plate indentation problem using VOS Figure 6.21 Load-Displacement curve for the flat plate indentation problem 0 10 20 30 40 50 60 70 80 0 5 10 15 20 Lo ad ( kN ) Displacement (mm) VOS ALE Exp. 102 Figure 6.22 Load-Displacement curve for the hemispherical plate indentation problem Figure 6.23 Comparison of the load-displacement curves of the flat and the hemispherical indenters 0 10 20 30 40 50 60 70 80 90 100 0 5 10 15 20 25 30 Lo ad ( kN ) Displacement (mm) VOS ALE Exp. 0 20 40 60 80 100 0 5 10 15 20 25 30 Lo ad ( kN ) Displacement (mm) Flat Indenter Hemispherical Indenter 103 Figure 6.24 von-Mises strain distribution using the flat-face indenter in the plate indentation problem Figure 6.25 von-Mises strain distribution using the hemispherical indenter in the plate indentation problem 104 6.6 V-Bending In this application, the V-bending process is simulated to demonstrate the capability of the VOS method in dealing with sharp corner contacts and in handling complex changes in nonlinear boundary conditions during the stages of the simulation. The simulation results from the VOS method are compared to experimental results [28] and numerical results from ANSYS [76]. The material used for the sheet metal in this application is Aluminum 5052-H32. A tensile test following the ASTM standards is performed using a Tinius Olsen Testing Press to obtain the stress-strain curve of the material. The material behavior is approximated to a bilinear elastic-plasticmodelwithaYoung‟smodulusE = 70 GPa and Poisson‟sratiov = 0.33, initial yield stress Sy = 120 MPa and tangent modulus ET = 580 MPa. The punch and the die are made from a high strength Mild Steel 1018 that is machined to the required dimensions. The process consists of a 90 o V-shaped punch pressing against an aluminum sheet into a 90 o V-shaped die with a total punch stroke of 24.5 mm in the downward direction. Figure 6.26 shows the setup of the V-bending application with the dimensions of each part. The experiment is performed using the Tinius Olsen Testing Press. All the parts shown in Figure 6.26 are machined for this test. The punch load and displacements are recorded during the experiment to be later compared with the numerical results. The same process is simulated using ANSYS with the workpiece meshed into 144 uniform 4-noded elements and treated as a plain strain problem. Due to symmetry, the VOS simulation is performed for half of the model using a total number of 8,057 elements for the whole domain with 1,495 of these elements being initially filled by the metal sheet in its initial configuration. The Eulerian grid domain is meshed using a uniform grid at some regions and a non-uniform grid at other regions. This is done to optimize the number of elements at regions that are expected to be empty during the course of the deformation. Considering the right half of the V-bending model, the mesh intensity is set to be low at the top left and bottom right corners of the domain since those two regions are expected to be empty during the entire process. The die is modeled as a rigid body using line segments fixed in all directions. The punch is modeled 105 as a rigid body using three line segments, all fixed in the horizontal direction. Each line of the punch and the die is defined as smaller segments at the elemental level. Two contact pairs are defined in the simulation. The first one is between the punch and the workpiece and the second one is between the die and the workpiece with a coefficient of friction of 0.1 for both pairs. The contact algorithm in the VOS simulation defines the contact at the elemental level, and the contact is detected once the rigid line elemental segments contact a free surface interface. After the initiation of the contact at each free surface element, the coordinates of the free surface element line and the rigid line segment are monitored to detect any update in the contact condition. A connectivity weighting factor of 1 is given to the bottom boundary of the metal sheet at the region under the punch. This ensures that the sheet will maintain its thickness since the part in contact with the punch dominates the connectivity process in the upper boundary of the workpiece. All subsequent regions share the same weighting factor at both sides of the workpiece. Figure 6.27 shows the evolution of the deformed shape for the metal sheet during the experiment and the VOS simulation. Comparing the deformed shapes of the metal sheet, it is possible to see that the deformed shape is accurately presented in VOS. The contact with the corner point remained during the entire process with hardly any penetration which is difficult to maintain when using traditional Lagrangian or ALE formulations. The contact algorithm is part of the VOS method which incorporates the contact boundary conditions as a part of the free surface detection scheme. This arrangement allows the method to eliminate any penetration during the simulation process. Figure 6.28 shows the load-displacement curve obtained from the VOS simulation, ANSYS and the experimental results. In general, the load-stroke curve has a similar behavior in the three solutions. However, at the last stages of the process, the punch load starts to increase at a very high rate due to the fact that most of the workpiece boundary comes in contact with the punch and the die. The VOS and ANSYS curves have load values that are much higher than the experimental curve. This difference is attributed to the fact that the die in the experiment is designed with a stress relaxation gap to prevent overloading the load cell of the testing press at maximum punch displacement. Another 106 observation is that before the last stages of the process, the load values of the VOS method are slightly higher than in ANSYS and the experiment. This is due to the approximation of the material model that was used in the VOS simulation. The development of the von-Mises stresses and strains at different punch strokes using the VOS method are shown in Figure 6.29 and Figure 6.30, respectively. It can be seen from the distribution that the maximum values for the stress and strain occur at the punch land towards the center line (line of symmetry). This distribution is very logical since that region is where the sheet undergoes the maximum compression (top boundary) and tension (bottom boundary). Table 6.2 compares the maximum values for the von- Mises stress and strain from the VOS method and ANSYS. The results from VOS are very close to the values obtained using ANSYS. The slight differences in values are due to the approximation of the material modeling and the difference in the formulation between the Lagrangian formulation used in ANSYS, and the Eulerian formulation used in VOS. The VOS method was validated in the above example by comparison with experimental work [28] and numerical analysis using commercial FE software (ANSYS). The development of the sheet metal deformation was presented accurately using VOS, and the shape of the free boundary was very similar to the experiment. Even though there are no mesh distortion problems in such applications, simulating such process using the VOS method highlights the capability of the method in dealing with large rotations and change of free boundaries along the whole surface of a workpiece rather than just one or two sides. The VOS also showed how easily it can deal with sharp corner contacts without the need for any special contact algorithm as it deals with it using the default contact parameters. In addition to estimating the shape of the free surface, different results were validated by comparing VOS to ANSYS, showing that the VOS method can be used effectively to capture the free surface in an Eulerian grid in addition to calculating accurate results and distribution of the stresses and strains. 107 90 o 66 mm 9.5 mm 115 mm 6.5 mm 90 o 25.4 mm 25.4 mm63.5 mm 19 mm 19 mm 6.5 mm 6.5 mm 6 4 m m 3.2 mm 3.2 mm Figure 6.26 Geometrical setup for the V-bending problem 108 (a) 6.35mm punch stroke (b) 12.7mm punch stroke (c) 19.05mm punch stroke (d) 25.4mm punch stroke Figure 6.27 Evolution of the deformed shape from Experiment and VOS in the V- bending problem 0 10 20 30 40 50 60 70 80 0 5 10 15 20 25 30 P u n ch L o a d ( K N ) Punch Stroke (mm) VOS Exp. ANSYS 0 5 10 15 0 5 10 15 20 Figure 6.28 Load-Stroke curve for the V-bending problem 109 (a) 6.35mm punch stroke (b) 12.7mm punch stroke (c) 19.05mm punch stroke (d) 25.4mm punch stroke Figure 6.29 von-Mises stress distribution for the V-bending problem 110 (a) 6.35mm punch stroke (b) 12.7mm punch stroke (c) 19.05mm punch stroke (d) 25.4mm punch stroke Figure 6.30 von-Mises strain distribution for the V-bending problem Table 6.2 von-Mises stress and strain comparison between VOS and ANSYS for the V-bending problem Punch Stroke 6.35mm Punch Stroke 12.7mm Punch Stroke 19.05mm Punch Stroke 25.4mm Stress (VOS) Stress (ANSYS) Strain (VOS) Strain (ANSYS) 208.8 MPa 216.4 MPa 0.133763 0.130869 238. MPa 231.6 MPa 0.23644 0.228923 241.7 MPa 236.6 MPa 0.368352 0.359114 244.5 MPa 242.3 MPa 0.516223 0.462803 111 6.7 Three-Point Bending In this application, the three-point bending (U-bending) is simulated to show the capability of the VOS method in dealing with large curved contact areas. Since the punch in this example is a cylinder and the dies are rollers, we expect the contact regions to be highly curved throughout the simulation in which the VOS can accurately predict the deformed shape of the workpiece while keeping the free surface elements intact. The VOS simulation is compared to experimental results. The process consists of a punch with a cylindrical tip that presses against a sheet metal that is supported by two rollers. The rollers are connected to a supporting frame that allows the rollers to rotate freely. The punch and the rollers are all made from mild steel1018withYoung‟smodulus E = 205 GPa, yield strength Sy = 370 MPa, ultimate strength Su = 440 Mpa and density ρ = 7870 kg/m 3 . The sheet is made from Aluminum 5052-H32withYoung‟smodulusE =70GPaandPoisson‟srationv = 0.33. A tensile test following the ASTM standards was performed using a Tinius Olsen Testing Press to obtain the stress-strain curve of the material. The material behavior was approximated to a bilinear elastic-plastic model for the VOS simulation with yield stress Sy = 195 MPa and a tangent modulus ET = 280 MPa. The geometrical setup of the problem is shown in Figure 6.31. The process was performed for a total punch stroke of 76.2 mm with a constant velocity of 6.35mm/min. The experiment was performed using the Tinius Olsen Testing Press with all the parts being machined for this test. The punch load was recorded during the experiment for different punch strokes. The VOS simulation was performed for half of the model due to symmetry. An 80 x 105 window, meshed into 8766 elements, was used for the Eulerian domain with 539 elements filled by the initial configuration of the sheet. The punch and the roller were modeled as a quarter and a half of a circle using segments of rigid lines. The only boundary conditions assigned to the workpiece were the symmetry boundary conditions at the line of symmetry. The average coefficient of friction for aluminum used in this simulation is 0.1. Similar to the V-bending problem (Section 6.6), a connectivity weighting factor of 1 is given to the bottom boundary of the sheet at the region under the punch. This ensures that the sheet will maintain its thickness since the part in contact with 112 the punch dominates the connectivity process at the upper boundary of the workpiece. All subsequent regions share the same weighting factor at both sides of the workpiece. Figure 6.32 shows the evolution of the deformed shape at different punch strokes during the VOS simulation and the experiment. It can be seen that the deformed shape from the VOS simulation is presented accurately compared to the deformed shape from the experiment. The punch load-stroke curve obtained from the experiment is compared to the one obtained from the VOS simulation in Figure 6.33. The curve from the VOS results matched very well with the experimental curve. At the beginning, the punch load increases with the punch stroke until it reaches a peak value when the workpiece bend flange becomes tangent to the die face. Then the load starts to decrease as the workpiece slides to the side of the rollers and the punch. The development of the von-Mises stresses and strains at different punch strokes using the VOS method are shown in Figure 6.34 and Figure 6.35, respectively. The stresses and strains start to develop at the region under contact with the punch. As the process continues and the contact area with the punch increases, stresses and strains starts to generate and extend towards the end of the workpiece. The values are always higher under the punch land because that is where the workpiece undergoes the maximum tension (bottom boundary) and compression (top boundary). The simulation of the three-point bending problem using the VOS method was com compared to experimental results obtained using a similar setup to the V-Bending problem in Section 6.6 [28]. The evolution of the deformed shape was modeled accurately compared the recoded deformation shapes from the experiment. The VOS method handled the change in contact condition and the large rotations effectively. The change in the free boundary was tracked using the VOS method during the course of the simulation. The VOS method also showed how easily it can deal with sharp corner contacts without the need for any special contact algorithm as it deals with it using the default contact parameters. 113 R19 mm 152.5 mm 6.35 mm R19 mm R10 mm R32 mm 202 mm 115 mm 50 mm69 mm 56.5 mm 63 .5 m m Figure 6.31 Geometrical setup for the three-point bending problem 114 (a) 12.7 mm punch stroke (b) 25.4 mm punch stroke (c) 38.1 mm punch stroke (d) 50.8 mm punch stroke (e) 63.5 mm punch stroke (f) 76.2 mm punch stroke Figure 6.32 Evolution of the three-point bending deformed shape from the Experiment and VOS 115 Figure 6.33 Load-Stroke curve for the three-point bending problem (a) 12.7 mm punch stroke (b) 25.4 mm punch stroke (c) 38.1 mm punch stroke (d) 50.8 mm punch stroke (e) 63.5 mm punch stroke (f) 76.2 mm punch stroke Figure 6.34 von-Mises stress distribution for the three-point bending problem 0 1000 2000 3000 4000 5000 6000 7000 0 20 40 60 80 P u n ch L o ad ( N ) Punch Stroke (mm) VOS EXP. 116 (a) 12.7 mm punch stroke (b) 25.4 mm punch stroke (c) 38.1 mm punch stroke (d) 50.8 mm punch stroke (e) 63.5 mm punch stroke (f) 76.2 mm punch stroke Figure 6.35 von-Mises strain distribution for the three-point bending problem 117 CHAPTER 7 METAL CUTTING APPLICATION Metal cutting is considered to be one of the most complicated manufacturing processes. Even though the basic mechanics are almost the same for all metal cutting operations, each case is considered different due to the complexity that comes from the dependency of the process on many geometrical factors and material properties. Most of the research done in metal cutting used to be based on experimental work but since the advancements of the Finite Element Method (FEM) and numerical simulations in the last few decades there has been more work on studying the machining process by means of numerical analysis. Compared to experiments, presenting accurate results from numerical simulations is still considered a challenging task due to the coupling of high temperatures and high strain rates, and the difficulty in modeling the friction and the contact between the chip and the tool. Due to the large plastic work and high friction, massive amounts of thermal energy are generated causing the temperature to rise in the metal and the tool. This means that an accurate numerical simulation must consider appropriate heat transfer analysis along with the deformation and stress analyses [24, 77, 82]. The natural FE formulation choice in metal forming applications, including machining, is the Lagrangian approach because the mesh is attached to the material points and undergoes the same deformation. In metal cutting this takes care of an important issue which is the unconstrained material flow on the free boundaries during the formation of the chip throughout the machining process. This eliminates the need for a priori assumption for the shape of the chip at steady state stages. This is important since the shape of the chip is highly dependent on the friction coefficient between the chip and thetool,thecuttingspeedandtheshapeofthetool(rakeangle,chamfered,blunt…etc) [77, 83, 84]. The most common way to simulate metal cutting with Lagrangian formulation is to use a nodal separation method in conjugation with a failure criterion [85-89]. In this 118 approach, the line of separation at the cut depth is defined prior to the analysis and the failure criterion is applied to the nodes/elements along that line. As the tool advances through the metal, the nodes starts to detach in front of the tip of the tool based on the failure criterion. Usually the separation aspect of the analysis is modeled in a similar fashion to an elastic-plastic crack propagation problem. One of the drawbacks of this approach is that the material separation happens at a small distance in front of the tool tip rather than at the tip itself [86-88]. Physically, this is an unrealistic approach since such cracks are not observed in actual cutting of ductile materials. Additionally, the process will be modeled with discrete jumps of the crack tip from one node to the next one which requires a very fine mesh in order to minimize the error of the material deformation history results at the vicinity of the tool tip. Although this might be reasonable in some cases where the distance between the tool tip and crack tip is too small, it can only perform well for tools with finite nose radius and fails in applications with chamfered or blunt tools [25, 77]. With chamfered and blunt tools, there exists a dead metal zone under the chamfered length or the curved radius that cannot be modeled properly with nodal separation methods. In addition to that, one must pay attention to the unzipped nodes along the parting line since the line might be pushed out of its position as the tool advances. Another difficulty with this approach is the choice of a failure/damage criterion with the proper values. Most of these criterions are based on parameters that are obtained from experiments and once the experimental parameters or the material used in the process are changed, the criteria will not be valid. Different numerical simulations were carried out using different values and different criterions [88-90]. It is noteworthy to know that Huang and Black [91] concluded after studying different criteria that the type of criterion used in the analysis is not significant but the magnitude of the failure critical value has a great effect on the distribution of the stresses and strains. A viable substitute for the nodal separation method could be achieved by simulating the cutting process as an indentation problem by having the tool pressing against the workpiece causing the material (elements) to flow around the tool [75, 84, 92, 93]. This approach requires extensive remeshing around the tool tip at very small time increments due to the enormous distortion of the elements around the tool. Although this approach overcomes most of the difficulties and problems related to the nodal separation method, it 119 does have its own drawbacks. The constant remeshing of the workpiece in addition to the fine mesh required around the tool leads to huge computational cost [93]. Another negative aspect of this approach is the interpolation of the material variables from the old mesh to the new mesh after each remeshing step. Because of the difference between the old and the new mesh configuration, the approximation that is used in the interpolation process can lead to error accumulation at each remeshing step which may cause deterioration of the results [24]. A less common choice for metal cutting simulations is the Eulerian formulation approach. This approach eliminates some of the problems that arise when using any of the Lagrangian approaches. Since the mesh is fixed in space and the material is allowed to flow freely around the tool, the cutting process can be modeled as an indentation problem without the need for a failure criterion. Additionally, since no mesh distortion occurs in an Eulerian grid, no remeshing steps are required. This makes the Eulerian approach computationally efficient and appropriate for modeling the flow around the tool tip. The main downside of this formulation is that it is not easily adaptable for modeling the unconstrained flow of the material at the free boundaries as the chip evolves during the process [25, 94-96]. To overcome this problem, a priori assumption for the shape of the chip at its steady state is required and the process is simulated only for the steady state stages to estimate the values of the material variables around the chip-tool area. A proper assumption of the chip shape is very difficult to obtain since it depends on many factors and it cannot be predicted accurately without performing experiments with the same cutting conditions. An alternative approach is ALE formulation. Some of the work done with ALE in metal cutting focused on improving the nodal separation method by allowing a more natural treatment of the metal cracking [23]. This is achieved by minimizing the “unrealistic” distance between the tool tip and the separation nodes.ALE can also be used in conjunction with remeshing techniques to minimize the number of remeshing steps. The ALE can also be used to take full advantage of the Lagrangian and Eulerian methods and eliminate the need for remeshing and nodal separation. This requires a proper mesh motion scheme that allows the free flow of the material around the tool 120 edge, as in the Eulerian approach, and the unconstrained flow of the material at the free boundaries, as in the Lagrangian approach, to develop the chip shape naturally [24, 25, 77]. The ALE approach proved to be viable once a proper mesh motion scheme has been implemented in the analysis. Establishing a general and reliable mesh motion scheme for ALE applications in metal cutting is not, however, an easy task although some effort has been done in this direction [23-25]. In this chapter we present a study of the metal cutting process with one sharp tool, two chamfered tools and one blunt tool that are used under the same cutting conditions. The study is carried numerically using the VOS method and results are compared to the numerical work done by Movahhedy et al. [77] using ALE as well as experimental results by Ren and Altintas [97] that were performed for three of the four different tools to compare the forces acting in the process. Figure 7.1 shows the schematic of the metal cutting process for a general case with a chamfered tool. The thrust and cutting forces are defined by Ft and Fc, respectively, where Fs and Fn represent the forces acting along and normal to the shear plane, respectively. The rake and chamfer angles are defined as α0 and α1. The length of the chamfered section is defined as bcf. α0 α1 bcf Ft Fc R Vc R Fn Fs Figure 7.1 Schematic view of metal cutting process with chamfered tool 121 7.1 Experimental Work and Numerical Setup In this work we used P20 steel for the workpiece and carbide of ISO S10 class for the tools. Three sets of orthogonal cutting tests were performed using one sharp and two chamfered tools [97]. The P20 mold steel disks were turned in thrust orthogonal mode on a CNC turning center. All of the tests were performed with a 0 o rake angle and the only differences between the three cases were in the chamfer angles and lengths. The cutting speed was set at 240 m/min with a chip load of 0.1 mm. The same three cases were simulated using the VOS method in addition to an extra case with a blunt (worn) tool with a radius of 0.075 mm. Table 7.1 lists the four cases used in the experiment and the simulation. Table 7.1 List of experimental and numerical metal cutting setups Case α1 bcf (mm) Results 1 2 3 4 0 o -10 o -25 o blunt - 0.0902 0.0863 Radius 0.075 Exp. & VOS Exp. & VOS Exp. & VOS VOS The flow stress model used for the P20 steel workpiece is dependent on the strain rate and the temperature in the following form [77, 97]: *1 lnm nA B C D ET (7.1) with the variable T * defined as: * room melting room T T T T T (7.2) 122 where A = 178.5 MPa, B = 462.4 MPa, C = 0.0438, D = 1.51 and E = 0.99 are the coefficient for the flow stress model. The strain sensitivity parameter m is equal to 0.169 and the temperature sensitivity parameter n is equal to 0.666. The melting temperature for P20 is 1480 o C. The rest of the P20 steel material properties are listed in Table 7.2 along with thermal properties of the carbide tool. All of the workpiece properties are listed for the P20 steel except for the convection heat transfer coefficient where we used a value for C45 steel [78] due to the difficulty in obtaining any values for P20 in the literature. The tool was considered rigid for the stress analysis part in the numerical simulation and only the thermal heat transfer analysis was performed for the tool. It is difficult to apply perfect contact conditions in metal cutting simulations and different values and settings have been used in previous work done in the area. The general consensus in the literature is that a mean friction coefficient value can be obtained from the ratio of the experimental values of the thrust and cutting forces with the same conditions. The mean values obtained from the experimental work [97] fall within the range of values presented in the literature [27, 77, 78, 92, 98]. Table 7.2 Material properties for the workpiece and tool used in metal cutting Material Properties P20 Workpiece Carbide Tool Young‟sModulus(GPa) Poisson‟sRatio Thermal Conductivity, k (W m -1 o C -1 ) Specific Heat (J Kg -1 o C -1 ) Thermal Expansion Coefficient ( o C -1 ) Thermal Convection Coefficient, h (W m -2 o C -1 ) 260 0.3 51.5 470 1.4x10 -6 20 - - 120 343.3 5.2x10 -6 - 7.2 VOS Simulation of Metal Cutting A combination of uniform and non-uniform mesh was used for the Eulerian domain. Regions outside of the deformation zone were constructed using coarse mesh to minimize 123 the number of elements that are not used in the analysis (Figure 7.2). The left region of the domain, where the velocity is applied to model the movement of the workpiece towards the tool, was constructed using a uniform grid with fairly sized elements. In a similar way, the right region where the material exits the machining section was also constructed with fairly sized uniform grid. This is because the material flow is steady in the horizontal direction and there is no change in flow movement or in the general shape. The region enclosing the tool was also meshed with relatively medium sized elements except toward the edges where the contact between the tool and the workpiece occurs. A fine mesh was used around the contact region to get better results in the thermal analysis. Finally, the middle part is the most important part of the domain since it is the region where the chip deformation is expected to take place. This part was constructed with very fine mesh to allow more accurate chip tracking results and better chip shape. Since the VOS method is based on constructing linear lines at the interface of the free surface elements, the finer the mesh the more accurate the free surface shape results become. Figure 7.2 shows the evolution of the chip deformation through the cutting process for the sharp tool. In our VOS FE-code, the stress and the thermal analyses were performed in a staggered fashion. At first, the stress analysis is performed with the workpiece and tool being at room temperature. The temperature data calculated from the internal heat generation in stress analysis is then transferred and the thermal analysis is performed. After that, the new stress analysis step is performed using the temperature data calculated from the thermal analysis as an input. Although this approach does not represent the physical dynamics of the process, it is not expected to introduce any other approximations beyond those adapted in a staggered approach especially since the material properties dependency on the temperature is weak and the time step used is very small. In addition to that, the scheme is easy to implement and the results from this work and from [25, 77] proved to be stable and acceptable. 124 Figure 7.2 Evolution of chip deformation in the VOS metal cutting simulation 7.3 Results and Discussion When studying metal cutting process with chamfered and blunt tools, one of the main points of interest is the existence of the dead metal zone under the chamfer or the worn part of the tool. Experimental and analytical studies stated that a dead metal zone is formed in those regions and that it is mostly independent of the cutting conditions and it depends only on the geometrical profile of the chamfered part [25, 97, 99, 100]. The explanation is that this zone is supposed to fill the part under the chamfer which makes the tool work like a sharp tool with dead metal zone acting as the effective cutting edge of the tool. To examine this phenomenon, we study the velocity field of the material around the tool edge. Figure 7.3 shows the velocity profile of the four different tools used in the VOS simulations. For the sharp tool, it is clear that there is no dead metal zone except for a very small region where the material comes in contact with the tool corner causing the material flow to slow down. For the chamfered tools, it is clear that part of the material is 125 trapped under the chamfer, with low velocities, and almost completes the missing tool nose. The zone is clearer for the case where the chamfer angle is -25 o since the area of the missing nose is bigger which shows that it is dependent on the geometry/chamfer angle. As for the blunt tool, a dead metal zone occurs but to a lesser extent since the curved part of the tool allows a smoother flow of the material. Another interesting observation is the size of the region where material velocity starts to increase from relative stagnation to the working velocity of the workpiece under the tool flank. For the sharp tool, the material almost directly picks up the working velocity. For the chamfered cases, the material takes more time to get back to the working velocity again which also increases with increased chamfer angles. The blunt tool has the largest extent of the stagnation region because the curved part of the tool reaches beyond the corners of the chamfered tools. This occurs because of the dead metal zone which develops a tool nose that is not on the same level of the tool flank causing additional straining at that region. The distribution of the effective stresses is demonstrated in Figure 7.4 for the four different tool geometries. The shear zone is clear in these contours which possesses the highest stress values and is concentrated between the beginning of the chip curve at the free surface to the edge of the tool or the dead metal zone. The stress value is almost constant for all cases along the shear plane. It is also clear that for blunt and higher chamfer angle tools; high stresses extend beyond the shear plane towards the flank of the tool because they allow a smoother flow of the material around the tool edges which also drags the high stresses out of the shear plane causing some stress relieving towards the center of the shear zone. The same observations regarding the effect of the dead metal zone and the shear zone can be seen in the distribution of the plastic strain rate in Figure 7.5. The strain rates are concentrated at higher values along the shear plane than, where in the rest of the workpiece strain rates are much smaller and fairly constant. The distribution also has large values towards the flank of the tool due to the extrusion of the material under the dead metal zone. The effect of the chamfer angle or a blunt tool can also be observed here 126 by noting how longer it takes for the material strain rate to drop to lower values after is passes the tool flank area. Although the chip evolution has a weak dependency on the chamfer angle or radius of the blunt, these factors do have a significant effect on the tangential and, to a lesser extent, on the normal forces. Figure 7.6 compares the maximum stresses in the tangential and normal directions around the tool nose obtained from the VOS simulation and the ones from ALE [77]. It is clear that the increase of the chamfer angle causes a slight increase in the normal stresses acting along the cutting direction. On the other hand, the stresses in the thrust direction, that is tangent to the tool rake face, show an evident increase as we go from a sharp tool to a chamfered and a blunt tool. The increase of the chamfer angle has a more dominant effect in this direction compared to the normal direction. This is logical since the presence of a chamfer causes the buildup of the dead metal zone which strains the material to a larger extent compared to the sharp tool. Figure 7.7 and Figure 7.8 compares the cutting and thrust forces obtained from VOS, ALE [1] and the experimental work, respectively. It is obvious that both forces increase with the increase of the chamfer angle. As for the stresses, the force in the cutting direction is less affected by the chamfer angle and the tool radius in contrast to the forces in the thrust direction which display a more pronounced increase in values. The same observation can be seen from the experimental results which exhibit similar behavior to that of the VOS and ALE results. The forces calculated for the blunt tool are higher than the ones for the sharp and chamfered tools. A rational explanation for this is that the higher and lower portions of the worn tool curve may act as a chamfered edge with a small and large chamfer angles, respectively. There is a notable difference between the numerical and the experimental values which is attributed to the simplification of the friction conditions used in the numerical simulations and the uncertainty of the material model and its thermal properties used in the VOS and ALE simulations. The results from the VOS simulation showed good agreement in values and behavior with the ones obtained from the proposed ALE approach by Movahhedy et al. [77]. The effect of the chamfer angle and blunt radius on the cutting and thrust forces in the VOS and ALE simulation possesses the same behavior of the experimental results [97]. Some 127 differences in the values between two approaches and the experiment can be credited to the assumptions and simplifications made in the friction and material model. In general, the VOS shows a great promise in handling continuous metal cutting applications and the differences between the VOS and experimental results were due to the approximation in the friction and material modeling. It is important to note that the VOS method is only suitable for continuous chip formation and it cannot handle segmented chip formation in its current state. V mm/sec Figure 7.3 Velocity fields for sharp tool, -10 o chamfered tool, -25 o chamfered tool and a blunt tool 128 σeff MPa Figure 7.4 Effective stress distribution for sharp tool, -10 o chamfered tool, -25 o chamfered tool and a blunt tool 129 Rate 1/sec Figure 7.5 Plastic strain rate distribution for sharp tool, -10 o chamfered tool, -25 o chamfered tool and a blunt tool 130 Figure 7.6 Comparison of maximum normal and tangential stresses between VOS and ALE simulations Figure 7.7 Comparison between the cutting forces obtained using VOS, ALE and experimental results 500 1000 1500 2000 2500 0 10 20 30 S tr es s (M P a ) Chamfer Angle Sxx ALE Sxx VOS Syy ALE Syy VOS 150 200 250 300 350 0 10 20 30 C u tt in g F o rc e (N ) Chamfer Angle Experiment ALE VOS 131 Figure 7.8 Comparison between the thrust forces obtained using VOS, ALE and experimental results 50 100 150 200 250 0 10 20 30 T h ru st F o rc e (N ) Chamfer Angle Experiment ALE VOS 132 CHAPTER 8 CONCLUSIONS AND FUTURE WORK 8.1 Summary and Conclusions The simulation of metal forming and machining processes is considered to be one of the most challenging topics in the area of solid mechanics. Both traditional, Lagrangian and Eulerian, formulations have many drawbacks that hinders their use in this field of research. Different schemes and formulations have been presented in the literature to tackle these drawbacks and overcome the difficulties that arise when using FE in simulating metal forming applications. Re-meshing proved to be very effective in dealing with problems related to mesh distortion but on the expense of accuracy, computational time and problem dependency. The ALE formulation on the other hand showed that it may be used effectively in overcoming all the drawbacks associated with other formulations and schemes. The main issue in ALE is that devising a general and proper mesh motion scheme that can be used in wide range of applications is very difficult to attain, hence most of the mesh motion schemes are limited to certain class of applications. In the present work we proposed a new method that is based on the Volume of Fluid (VOF) method [37] to be used with Eulerian formulations to calculate the volume of the solid material, hence called the Volume of Solid (VOS) method, and to track the unconstrained flow of the material at the free boundaries. A general review of the current FE formulations used in metal forming numerical simulations was presented in Chapter 1. In Chapter 2, a comprehensive survey of the VOF method and the different schemes used in the literature was presented highlighting the drawback/advantages of each scheme. This critical study is the basic ground for formulating the VOS method since understanding the different aspects of the VOF method helps in recognizing the proper approach to be used in formulating the VOS method for metal forming applications. 133 The Eulerian finite element equations were presented in Chapter 3 for quasi-static and dynamic analyses. The complete formulation and discussion of the proposed VOS method was presented in Chapter 4 for the case of a uniformly constructed mesh and for the general case of a non-uniform mesh. All topics related to the implementation of the VOS method were discussed in Chapter 5. These topics included the time step stability, updating the velocity field, tracking material point related properties and connectivity of the free surface interfaces. A contact module was also presented in addition to a heat transfer analysis procedures. Several metal forming and metal cutting simulations were performed in Chapters 6 and 7. The results were validated with different results obtained from the literature, experimental work and commercial FE softwares. The contributions from this work may be summarized in the following points: A complete derivation of the Eulerian FE equations based on the principle of virtual work was presented for quasi-static and dynamic analyses. A new method was introduced for the treatment of the convective terms in the Eulerian equilibrium equations based on the work by Bayoumi and Gadala [28, 68]. Unlike the commonly used operator split method where the equations are uncoupled into an advection and a non-advection terms, the new approach utilizes the definition of the VOS method in dealing with these terms without uncoupling the solution. The complete formulation of the VOS method was developed for uniform and non-uniform grids. The formulation included the computation of the fractional volume values, calculation of the free surface directional vector and the construction of the free surface interface. For the uniform mesh case, a PLIC- VOS method is used in addition to a general classification of all the possible free surface interface shapes. A new filling pattern scheme is used for the general case of non-uniform mesh. The filling scheme is developed to construct the free surface with no dependency on the shape of the element. A new free surface connectivity process was employed in the VOS formulation to ensure that each free surface interface is connected to the interfaces at 134 neighboring elements. The process also accounts for every connectivity step to satisfy the conservation of mass and ensure no volume loss during the process. A contact algorithm was used to handle the contact between rigid and flexible bodies. The new module allows the modeling of a moving rigid body within the Eulerian domain and tracks the contact condition. Metal forming simulations have been limited mostly to applications that do not require having a moving punch/die. The new algorithm utilizes the VOS method by recording the movement of the free boundaries and the rigid tool which widens the range of applications that can be analyzed. Several metal forming applications were simulated using the developed VOS FE- code. The applications covered a wide range of applications highlighting the power of the VOS method overcoming problems like the mesh distortion and modeling the unconstrained flow of the material. Different shapes of tools were used in these simulations to demonstrate the effectiveness of the proposed VOS formulation in handling the deformed shape under these tools. Dealing with sharp angles was also handled easily with no penetration. A complex simulation for metal cutting analysis was performed using the VOS method. The simulation consisted of a heat transfer model that was coupled with the stress analysis to include the effects of the heat generation and the heat transfer between the chip and the cutting tool. The results from this model were compared to experimental and numerical results using ALE formulation. No material separation criterion was needed and the VOS method simulated the process in a more natural way as an indentation problem. 8.2 Future Work Even though the VOS formulation presented in this work proved to be a reliable FE formulation in simulating applications with large strains and deformations, there is still room for improvement. The VOS method is new and there is a lot of work to be done 135 towards the future progress of the method either by improving the current scheme or by implementing new techniques. Different topics could be explored for the future development of the VOS method in addition to implementation of new schemes and techniques. Such topics can be summarized as follows: Since the first introduction of the VOF method, different schemes were presented in the literature to improve the efficiency of the method in tracking the free surface of the fluid. A diverse range of free surface re-construction methods were used for different application. Since the VOS method presented in this work acts as an introduction of the method into solid mechanics, more exploration of the VOF methods could lead to better and improved VOS formulations. Extending the VOS method to handle two (or more) flexible bodies; rather than one flexible and one rigid; would introduce a new class of applications that could be analyzed with this method. One of the main advantages of the proposed VOS method in this work is that it defines cells that are not occupied by the material as “true” empty cells, rather than assigning them fictitious values.Thismakes the formulation of a modeling two or more flexible bodies much easier. Improvements could be made on the connectivity process which is considered an important topic in the VOF literature. Additional controlling parameters could be introduced to improve the efficiency of the current connectivity process. New connectivity schemes may also be developed to enhance the connectivity of the free surface interfaces. The connectivity could also be incorporated in the actual formulation of the VOS equations to ensure a more natural connectivity rather than having a separate connectivity step performed after each time step. If a formulation for two or more flexible bodies is developed, then a new connectivity process must be derived to account for both bodies if they come into contact and ensure the mass conservation of both. An extension of the current contact module can be done to include the contact between two flexible bodies. This would widen the range of analyses that could be performed including tool design and flex-flex dynamic impact problems. If the 136 contact problem is elastic-plastic contact then modifications on the current module can be made to track the movement of the elastic body. If both bodies are expected to deform within the plastic range, then a new contact algorithm must be employed. The free surface connectivity of both bodies must be coupled when contact occurs to ensure the integrity of the contact is not lost during the connectivity process of each body. Extending the current VOS formulation to 3D analysis is a topic of interest. This would require modifications mainly in the free surface construction step. It would be challenging to classify all possible interfaces at the free surface for a uniform mesh. On the other hand, the filling pattern might be an easier choice to implement in 3D uniform mesh. In the case of non-uniform mesh, it is a difficult task to apply the filling pattern scheme as it would be complicated to choose a reference node and to calculate the free surface directional vector. Implementations of the VOF method in applications that involve fluid-structure interaction have been widely discussed. In these applications the structure is always fixed and the fluid is the main interest in the analysis. The fluid-structure interface can be explored from the VOS method point of view with the solid being the main point of interest. It would be interesting to derive a set of VOS-VOF equations that could be paired together in such applications. 137 REFERENCES [1] H.D. Hibbitt, P.V. Marcal and J.R. Rice, “A Finite Element Formulation for Problems of Large Strain and Large Displacement”, International Journal of Solids and Structures, Vol. 6, pp. 1069-1086, 1970. [2] R.M.McMeeking and J.R.Rice, “FiniteElementFormulations forProblems of Large Elastic-Plastic Deformations”, International Journal of Solids and Structures, Vol. 11, pp. 353-386, 1975. [3] K.J.Bathe,“FiniteElementProcedures”,PrenticeHallInc.,1996. [4] J. Cheng and N. Kikuchi, “A Mesh Rezoning Technique for Finite Element Simulations of Metal Forming Processes”, International Journal for Numerical Methods in Engineering, Vol. 23, pp. 219-228, 1986. [5] M.S. Gadala, G.A‟E. Oravas and M.A. Dokainish, “A Consistent Eulerian Formulation of Large Deformation Problems in Statics and Dynamics”, International Journal of Non-Linear Mechanics, Vol. 18, No. 1, pp. 21-35, 1983. [6] M. Abo-Elkhier, G.A‟E. Oravas and M.A. Dokainish, “A Consistent Eulerian Formulation of Large Deformation Analysis with Reference to Metal-Extrusion Process”,InternationalJournalofNon-linear Mechanics, Vol. 23, No. 1, pp. 37- 52, 1988. [7] A.M. Maniatty, P.R. Dawson and G.G.Weber,“AnEulerianElasto-Viscoplastic Formulation for Steady-State Forming Processes”, International Journal of Mechanical Sciences, Vol. 33, No. 5, pp. 361-377, 1991. [8] D. Demarco and E.N. Dvorkin, “An Eulerian Finite Element Formulation for ModellingStationaryFiniteStrainElasticDeformationProcesses”, International Journal for Numerical Methods in Engineering, Vol. 62, pp. 1038-1063, 2005. [9] E.N. Dvorkin and E.G.Petocz,“AnEffectiveTechniqueforModelling2DMetal Forming ProcessesUsinganEulerianFormulation”,EngineeringComputations, Vol. 10, pp. 323-336, 1993. 138 [10] E.N. Dvorkin, M.B. Goldschmit, M.A. Cavaliere, P.M. Amenta, O. Marini and W. Stroppianam “2D Finite Element Parametric Studies of the Flat-Rolling Process”, Journal of Materials Processing Technology, Vol. 68, pp. 99-107, 1997. [11] M.A. Cavaliere, M.B. Goldschmit and E.N.Dvorkin,“FiniteElementSimulation of the Steel Plates Hot Rolling Process”, International Journal for Numerical MethodsinEngineering”, Vol. 52, pp. 1411-1430, 2001. [12] M.A. Cavaliere, M.B. Goldschmit and E.N.Dvorkin,“FiniteElementAnalysisof Steel Rolling Processes”, Computers and Structures, Vol. 79, pp. 2075-2089, 2001. [13] D.A. Berazategui, M.A. Cavaliere, L. Montelatici and E.N. Dvorkin, “On The Modelling of Complex 3D Bulk Metal Forming Processes via the Pseudo- Concentrations Techniques. Application to the Simulation of the Mannesmann PiercingProcess”, International Journal forNumericalMethods inEngineering, Vol. 65, pp. 1113-1144, 2006. [14] S. Okazawa, K. Kashiyama and Y. Kaneko, “Eulerian Formulation Using Stabilized Finite Element Method for Large Deformation Solid Dynamics”, International Journal for Numerical Methods in Engineering, Vol. 72, pp. 1544- 1559, 2007. [15] B. Mehmandoust and A.-R.Pishevar,“AnEulerianParticleLevelSetMethodfor Compressible Deforming Solids with Arbitrary EOS”, International Journal for Numerical Methods in Engineering, Vol. 79, pp. 1175-1202, 2009. [16] W.F. Noh, “A Time Dependent Two-Space-Dimensional Coupled Eulerian- LagrangianCode”,MethodsinComputationalPhysics,Vol.3,pp.117-179, 1964. [17] C.W. Hirt, A.A. Amsden and J.L. Cook, “An Arbitrary Lagrangian-Eulerian ComputingMethodforAllFlowSpeeds”,JournalofComputationalPhysics, Vol. 14, pp. 227-253, 1974. 139 [18] J.WangandM.S.Gadala,“FormulationandSurveyofALEMethodinNonlinear SolidMechanics”,FiniteElementsinAnalysisandDesign,Vol.24,pp.253-269, 1997. [19] M.S. Gadala, M.R. Movahhedy and J. Wang, “On the Mesh Motion for ALE ModelingofMetalFormingProcesses”,FiniteElementsinAnalysisandDesign, Vol. 38, pp. 435-459, 2002. [20] M.S. Gadala and J. Wang, “ALE Formulation and its Application in Solid Mechanics”, Computer Methods in Applied Mechanics and Engineering, Vol. 167, pp. 33-55, 1998. [21] M.S.Gadala and J.Wang, “SimulationofMetal FormingProcesseswithFinite ElementMethods”,InternationalJournalforNumericalMethodsinEngineering, vol. 44, pp. 1397-1428, 1999. [22] M.S. Gadala and J. Wang, “Elasto-Plastic Finite Element Simulation of Rolling and Compression between Wedge-ShapedDies”,JournalofMaterialsProcessing Technology, Vol. 97, pp. 132-147, 2000. [23] L. Olovsson, L. Nilsson and K. Simonsson, “An ALE Formulation for the Solution of Two-Dimensional Metal Cutting Problems”, Computers and Structures, Vol. 72, pp. 497-507, 1999. [24] M. Movahhedy, M.S. Gadala and Y. Altintas, “Simulation of the Orthogonal Metal Cutting Process Using an Arbitrary Lagrangian-Eulerian Finite-Element Method”, Journal of Material Processing Technology, Vol. 103, pp. 267-275, 2000. [25] M.R. Movahhedy, “ALE Simulation of Chip Formation in Orthogonal Metal CuttingProcess”,PhDDissertation,theUniversityOfBritishColumbia,January 2000. 140 [26] J.L.F. Aymone,E.BittencourtandG.J.Creus,“Simulationof3DMetal-Forming Using an Arbitrary Lagrangian-Eulerian Finite Element Method”, Journal of Materials Processing Technology, Vol. 110, pp. 218-232, 2001. [27] O. Pantale, J.-L. Bacaria, O. Dalverny, R. RakotomalalaandS.Caperaa,“2Dand 3D Numerical Models of Metal Cutting with Damage Effects”, Computer MethodsinAppliedMechanicsandEngineering”,Vol.193,pp.4383-4399, 2004. [28] H.N.Bayoumi,“ArbitraryLagrangian-Eulerian Formulations for Quasi-Static and Dynamic Metal Forming Simulations”, PhD Thesis, Mechanical Engineering Department, University Of British Columbia (UBC), 2001. [29] T. Belytschko, Y.Y. Lu and L. Gu, “Element-Free Galerkin Methods”, International Journal of Numerical Methods in Engineering, Vol. 37, pp. 229-256, 2004. [30] Y.-M. Guo and K. Nakanishi, “A Backward Extrusion Analysis by the Rigid- Plastic Integralless-Meshless Method”, Journal od Material Processing Technology, Vol. 140, pp. 19-24, 2003. [31] Y. Guan, X. Wu, G. Zhao and P. Lu,“ANonlinearNumericalAnalysisforMetal- Forming Process Using the Rigid-(Visco)Plastic Element-FreeGalerkinMethod”, International Journal of Advanced Manufacturing Technology, Vol. 42, pp. 83-92, 2009. [32] I. Alfaro, J. Yvonnet, E. Cueto, F. ChinestaandM.Doblare,“MeshlessMethods withApplication toMetal Forming”,ComputerMethods inAppliedMechanics and Engineering, Vol. 195, pp. 6661-6675, 2006. [33] Y.Guan,G.Zhao,X.WuandP.Lu,“MassiveMetalFormingProcessSimulation Based on Rigid/Visco-Plastic Element-Free Galerkin Method”, Journal of Materials Processing Technology, Vol. 187-188, pp. 412-416, 2007. 141 [34] Y.J. Guan, G.C. Zhao, X.Wu and P. Lu, “Application of Rigid(Visco)-Plastic Element Free Galerkin Method to Simulation of Plane StrainForming”,Materials Science and Technology, Vol. 25, No. 3, pp. 345-350, 2009. [35] T.Belytschko,D.OrganandY.Krongauz,“ACoupledFiniteElement-element- freeGalerkinMethod”,ComputationalMechanics,Vol.17,No.3,pp.186-195, 1995. [36] L.-C Liu, X.-H. Dong and C.-X. Li, “Adaptive Finite Element-Element-Free Galerkin Coupling Method for Bulk Metal Forming Processes”, Journal of Zhejiang University SCIENCE A, Vol. 10, No. 3, pp. 353-360, 2009. [37] C.W. Hirt and B.D.Nichols,“VolumeofFluid (VOF) Method for the Dynamics of free Boundaries”, Journal of Computational Physics, Vol. 39, pp. 201-225, 1981. [38] C.W. Hirt, B.D. Nichols and N.C. Romero, “SOLA-A Numerical Solution Algorithm forTransient FluidFlows”,LosAlamosScientificLaboratory report LA-5852, 1975. [39] B.D.Nichols andC.W.Hirt, “Improved Free SurfaceBoundaryConditions for Numerical Incompressible-FlowCalculations”,JournalofComputationalPhysics, Vol. 8, pp. 434-448, 1971. [40] F.H. Harlow, “Hydrodynamic Problems Involving Large Fluid Distortion”, Journal of the Association for Computing Machinery, Vol. 4, pp. 137-142, 1957. [41] F.H.HarlowandJ.E.Welch,“NumericalCalculationofTime-Dependant Viscous IncompressibleFlowofFluidwithFreeSurface”,PhysicsofFluids, Vol. 8, pp. 2182-2189, 1965. [42] S. Zaleski, “Computation Of Multiphase Flow By Volume Of Fluid and Boltzmann Lattice Gas Methods”, SHORT COURSE in Modeling and Computation of Multiphase Flows, Zurich, Switzerland, 20-24 March 2006. 142 [43] J.M. Martinez, X. Chesneau and B. Zeghmati, “A New Curvature Technique Calculation for Surface Tension Contribution in PLIC-VOF Method”, Computational Mechanics, Vol. 37, pp. 182-193, 2006. [44] M.S. Kim and W.I. Lee, “A New VOF-Based Numerical Scheme for the Simulation of Fluid Flow with Free Surface. Part I: New Free Surface-Tracking Algorithmand ItsVerification”, InternationalJournal forNumericalMethods in Fluids, Vol. 42, pp. 765-790, 2003. [45] M.S. Kim, J.S. Park and W.I. Lee, “ANewVOF-Based Numerical Scheme for the Simulation of Fluid Flow with Free Surface. Part II: Application to the Cavity FillingandSloshingproblems”, International Journal for Numerical Methods in Fluids, Vol. 42, pp. 791-812, 2003. [46] P. Troch and J.D. Rouck, “An Active Wave Generating-Absorbing Boundary Condition for VOF Type Numerical Model”, Coastal Engineering, Vol. 38, pp.223-247, 1999. [47] Y.M. Shen, C.O. Ng and Y.H.Zheng,“SimulationofWavePropagationOvera Submerged Bar Using The VOF Method With a Two-Equation k-ε Turbulence Modeling”,OceanEngineering,Vol.31,pp.87-95, 2004. [48] M.R.H.Nobari,M.J.Ketabdari andM.Moradi, “AModifiedVolumeOfFluid Advection Method for Uniform Cartesian Grids”, Applied Mathematical Modelling, Vol. 33, pp. 2298-2310, 2009. [49] S.AfkhamiandM.Bussmann,“HeightFunctionsforApplyingContactAnglesto 2D VOF Simulations”, International Journal for Numerical Methods in Fluids, Vol. 57, pp. 453-472, 2008. [50] S.AfkhamiandM.Bussmann,“HeightFunctionsforApplyingContact Angles to 3D VOF Simulations”, International Journal for Numerical Methods in Fluids, Vol. 61, pp. 827-847, 2009. 143 [51] M. Malik, E.S.-C. Fan and M. Bussman, “AdaptiveVOFwith curvature-based refinement”,InternationalJournalforNumericalMethodsinFluids, Vol. 55, pp. 693-712, 2007. [52] A. Theodorakakos andG.Bergeles, “Simulation of SharpGas-Liquid Interface UsingVOFMethodandAdaptiveGridLocalRefinementAroundtheInterface”, International Journal for Numerical Methods In Fluids, Vol. 45, pp. 421-439, 2004. [53] J.P. Wang, A.G.L. Borthwick and R.E. Taylor, “Finite-Volume-Type VOF Method on Dynamically Adaptive Quadtree Grids”, International Journal for Numerical Methods in Fluids, Vol. 45, pp. 485-508, 2004. [54] F.XiaoandA.Ikebata,“An Efficient Method for Capturing Free Boundaries In Multi-FluidSimulations”,InternationalJournalforNumericalMethodsInFluids, Vol. 42, pp. 187-210, 2003. [55] F. Mashayek and N. Ashgriz, “A Hybrid Finite-Element-Volume-Of-Fluid Method for Simulating FreeSurfaceFlowsandInterfaces”,InternationalJournal of Numerical Methods in Fluids, Vol. 20, pp. 1363-1380, 1995. [56] N. Ashgriz, T. Barbat and G. Wang, “A Computational Lagrangian-Eulerian AdvectionRemap forFreeSurface Flows”, International Journal for Numerical MethodsinFluids”,Vol. 44, pp.1-32, 2004. [57] D. Gueyffier, A. Nadim, R. Scardovelli and S. Zaleski, “Volume-Of-Fluid Interface Tracking with Smoothed Surface Stress Methods for Three-Dimensional Flows”,JournalofComputationalPhysics, Vol. 152, pp. 423-456, 1999. [58] W. Yang, S.-H. Liu and Y.-L. Wu, “AnUnsplitLagrangianAdvectionScheme forVolumeofFluidMethod”,JournalofHydrodynamics,Vol.22,No.1,pp.73- 80, 2010. 144 [59] C.JanssenandM.Krafczyk,“ALatticeBoltzmannApproach for Free-Surface- Flow Simulations on Non-Uniform Block-Structured Grids”, Computers and Mathematics with Applications, Vol. 59, pp. 2215-2235, 2010. [60] S.OtsherandJ.A.Sethian,“FrontsPropagatingwithCurvatureDependentSpeed: Algorithms Based on Hamilton-JacobiFormulations”, Journal ofComputational Physics, Vol. 79, pp. 12-49, 1988. [61] G. Son, “Efficient Implementation of a Coupled Level-Set Volume-of-Fluid Method for Three-Dimensional Incompressible Two-Phase Flows”, Numerical Heat Transfer B, Vol. 43, pp. 549-565, 2003. [62] M. Sussman and E.G. Puckett, “A Coupled Level Set and Volume of Fluid Method for Computing 3D and Axisymmetric Incompressible Two-PhaseFlows”, Journal of Computational Physics, Vol. 162, pp. 301-337, 2000. [63] S.P. van der Pijl, A. Segal, C. Vuik and P. Wesseling, “A Mass-Conserving Level-Set Method for Modeling of Multi-PhaseFlows”,InternationalJournalfor Numerical Methods in Fluids, Vol. 47, pp. 339-361, 2005. [64] I.R. Park, K.S. Kim, J. Kim and S.H. Van, “A Volume-of-Fluid Method for Incompressible Free Surface Flows”, International Journal for Numerical Methods in Fluids, Vol. 61, pp. 1331-1362, 2009. [65] D.L. Sun and W.Q.Tao,“ACoupledVolume-of-Fluid and Level Set (VOSET) Method for Computing Incompressible Two-PhaseFlows”, InternationalJournal of Heat and Mass Transfer, Vol. 53, pp. 645-655, 2010. [66] N. Pathak, A. Kumar, A. Yadav and P. Dutta, “Effects of Mould Filling on Evolution of the Solid-LiquidInterfaceDuringSolidification”,AppliedThermal Engineering, Vol. 29, pp. 3669-3678, 2009. [67] J.H. Jeong and D.Y. Yang, “Finite Element Analysis of Transient Fluid Flow With Free Surface Using VOF (Volume-of-Fluid)Method andAdaptiveGrid”, 145 International Journal for Numerical Methods in Fluids, Vol. 26, pp. 1127-1154, 1998. [68] H.N. Bayoumi and M.S.Gadala,“ACompleteFiniteElementTreatmentfor the FullyCoupled Implicit ALE Formulation”, ComputationalMechanics, Vol. 33, pp. 435-452, 2004. [69] M.S.Gadala,“RecentTrends inALEFormulationand Its Application in Solid Mechanics”, Computer Methods in Applied Mechanics and Engineering, Vol. 193, pp. 4247-4275, 2004. [70] T. J.R. Hughes,W.K. Liu and T.K. Zimmermann, “Lagrangian-Eulerian Finite ElementFormulationfor IncompressibleViscousFlows”, Computer Methods in Applied Mechanics and Engineering, Vol. 29, pp. 329-349, 1981. [71] M.S.GadalaandJ.Wang,“Computational ImplementationofStress Integration in FE Analysis of Elasto-PlasticLargeDeformationProblems”,FiniteElementin Analysis and Design, Vol. 35, pp. 379-396, 2000. [72] K. Ravindran and R.W. Lewis, Finite Element Modelling of Solidification Effects in Mould Filling, Finite Elements In Analysis and Design, Vol. 31, pp. 99-116, 1998. [73] A.S. Usmani, J.T. Cross and R.W. Lewis, The Analysis of Mould Filling in Castings Using Finite Element Method, Journal of Materials Processing Technology, Vol. 38, Issues 1-2, pp. 291-302, 1993. [74] T. Nakayama, M. Mori, “An Eulerian Finite Element Method for Time- Dependant Free Surface Problems inHydrodynamics”, International Journal for Numerical Methods in Fluids, Vol. 22, pp. 175-194, 1996. [75] R.H. Wagoner and J.-L. Chenot, “Metal Forming Analysis”, Cambridge University Press, 2001. [76] ANSYS User’s Manual, Swanson Analysis System Inc., 1989. 146 [77] M.R. Movahhedy, Y. Altintas and M.S.Gadala, “NumericalAnalysis ofMetal CuttingWithChamferedandBluntTools”,JournalofManufacturingScienceand Engineering, Vol. 124, pp. 178-188, May 2002. [78] W. Grzesik, M. Bartoszuk and P. Nieslony, “Finite Element Modelling of Temperature Distribution in the Cutting Zone in Turning Processes with DifferentlyCoatedTools”,JournalofMaterialProcessingTechnology,Vol.164- 165, pp. 1204-1211, 2005. [79] N.R. Chitkara and W. Johnson, “Plane Strain Compression of Pre-Shaped Material between Wedge-Shaped Dies”, International Journal of Mechanical Sciences, Vol. 4, pp. 151-164, 1974. [80] K. D. Hur, Y. Choi and H. T. Yeo, “A Design for Cold Backward Extrusion Using FEAnalysis”, Finite Element inAnalysis and Design, Vol. 40, pp. 173- 185, 2003. [81] G.G. Corbett and S.R. Reid, “Quasi-static and Dynamic Local Loading of Monolithic Simply-Supported Steel Plate”, International Journal of Impact Engineering, Vol. 13, pp. 423-441, 1993. [82] A.J. Shih, “Finite Element Analysis of Orthogonal Metal Cutting Mechanics”, International Journal of Machine Tools and Manufacture, Vol. 36, No. 2, pp. 255- 273, 1996. [83] Y.K. Potdar and A.T.Zehnder,“MeasurementsandSimulationsofTemperature and Deformation Fields in TransientMetal Cutting”, Journal ofManufacturing Science and Engineering, Vol. 125, pp. 645-655, 2003. [84] H. Bil, S.E. Kilic and A.E.Tekkaya,“AcomparisonofOrthogonalCuttingData from Experiments with Three Different Finite Element Models”, International Journal of Machine Tools and Manufacture, Vol. 44, pp. 933-944, 2004. 147 [85] Z.-C. Lin and S.-P. Lo,“AStudyofDeformationoftheMachinedWorkpieceand Tool under Different Low Cutting Velocities with an Elastic Cutting Tool”, International Journal of Mechanical Sciences, Vol. 40, No. 7, pp. 663-681, 1998. [86] Z.-C. Lin and Y.-Y. Lin, “A Study of an Oblique CuttingModel”, Journal of Materials Processing Technology, Vol. 86, pp. 119-130, 1999. [87] Z.-C. Lin and Y.-Y.Lin,“AStudyofObliqueCutting for Different Low Cutting Speeds”, Journal of Materials Processing Technology, Vol. 115, pp. 313-325, 2001. [88] G.Shi,X.DengandC.Shet,“AFiniteElementStudyoftheEffectofFrictionin OrthogonalMetalCutting”,FiniteElementsinAnalysisandDesign, Vol. 38, pp. 863-883, 2002. [89] P.A.R. Rosa, P.A.F. Martins and A.G.Atkins, “Revisiting the Fundamentals of Metal Cutting by Means of Finite Elements and Ductile Fracture Mechanics”, International Journal of Machine Tools and Manufacture, Vol. 47, pp. 607-617, 2007. [90] E.-G Ng and D.K.Aspinwall, “Modelling ofHardPartMachining”, Journal of Materials Processing Technology, Vol. 127, pp. 222-229, 2002. [91] J.M. Huang and J.TBlack, “AnEvaluation of Chip SeparationCriteria for the FEM Simulation of Machining”, Journal of Manufacturing Science and Engineering, Vol. 118, pp. 545-554, 1996. [92] T.OzelandT.Altan,“DeterminationofWorkpieceFlowStressandFrictionat the Chip-Tool Contact for High-SpeedCutting”,InternationalJournalofMachine Tools and Manufacture, Vol. 40, pp. 133-152, 2000. [93] G.S. Sekhon and J.L. Chenot, “Numerical Simulation of Continuous Chip Formation during Non-Steady Orthogonal Cutting”, Engineering Computations, Vol. 10, pp. 31-48, 1993. 148 [94] J.S. Strenkowski and K.J.Moon, “FiniteElement Prediction ofChipGeometry and Tool-Workpiece Temperature Distribution in Orthogonal Metal Cutting”, Journal of Engineering for Industry, Vol. 112, pp. 313-318, 1990. [95] K.W. Kim and H.-C.Sin,“DevelopmentofaThermo-Viscoplastic Cutting Model using Finite Element Method”, International Journal of Machine Tools and Manufacture, Vol. 36, pp. 379-397, 1996. [96] J.-S. Wu, J.R. Dillon and W.-Y. Lu, “Thermo-Viscoplastic Modeling of Machining Process using a Mixed Finite Element Method, Journal of Manufacturing Science and Engineering, Vol. 118, pp. 470-482, 1996. [97] H.RenandY.Altintas,“MechanicsofmachiningwithChamferedTools”,ASME Journal of Manufacturing Science and Engineering, Vol. 122, pp. 650-659, 2000. [98] E. Ozlu,E.BudakandA.Molinari,“AnalyticalandExperimentalInvestigationof RakeContact and FrictionBehavior inMetalCutting”, International Journal of Machine Tools and Manufacture, Vol. 49, pp. 865-875, 2009. [99] Y. Karpat and T. O zel, “Analytical and Thermal Modeling of High-Speed Machining With Chamfered Tools”, ASME Journal of Manufacturing Science and Engineering, Vol. 130, 011001, 2008. [100] I.A. Choudhury, N.L.SeeandM.Zukhairi,“MachiningWithChamferedTools”, Journal of Materials Processing Technology, Vol. 170, pp. 115-120, 2005.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Eulerian formulation with Volume Of Solid (VOS) application...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Eulerian formulation with Volume Of Solid (VOS) application to metal forming Al-Athel, Khaled S. 2010
pdf
Page Metadata
Item Metadata
Title | Eulerian formulation with Volume Of Solid (VOS) application to metal forming |
Creator |
Al-Athel, Khaled S. |
Publisher | University of British Columbia |
Date Issued | 2010 |
Description | The volume of fluid (VOF) method has been used extensively in computational fluid dynamics (CFD) applications. The method is used to track the unconstrained flow of the fluid at the free surface in an Eulerian formulations framework. In this work we propose the use of the VOF method in nonlinear solid mechanics and metal forming applications, therefore the new method is called the Volume of Solid (VOS) method. The VOS method is to be used with Eulerian finite element (FE) formulations. In Eulerian formulation, the mesh is fixed in space and the material flows freely within the domain. This type of formulation overcomes many difficulties that are associated with other FE formulations such as mesh distortion and the need for re-meshing or a proper mesh motion scheme. The main drawback of the Eulerian formulation is the difficulty in tracking the unconstrained flow of the material at the free boundaries; therefore an additional scheme is needed to overcome this problem. The derivation of the VOS equations and its implementation in an Eulerian formulation along with applications in metal forming and metal cutting are the main goals of this research. In this work, we present the derivation of the Eulerian quasi-static and dynamic equilibrium equations. The complete development of the VOS method is presented for the commonly used uniform mesh using a piecewise linear interface construction method, and for a general case of a non-uniform mesh using a new filling pattern scheme. The formulation of the VOS method includes the calculation of the fractional volume value of the solid in each element and finding a directional vector to locate the interface of the free surface. The implementation of the VOS method is discussed by covering important topics that include the connectivity of the free surface interfaces in addition to introducing and modeling a new algorithm for the rigid-flexible contact problem. The formulation presented in this work is implemented in a VOS FE-code and used to simulate several metal forming and metal cutting applications. The VOS simulations are validated by comparing the results with other experimental and numerical work presented in the literature. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-09-28 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0071321 |
URI | http://hdl.handle.net/2429/28759 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2010-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2010_fall_al-athel_khaled.pdf [ 2.1MB ]
- Metadata
- JSON: 24-1.0071321.json
- JSON-LD: 24-1.0071321-ld.json
- RDF/XML (Pretty): 24-1.0071321-rdf.xml
- RDF/JSON: 24-1.0071321-rdf.json
- Turtle: 24-1.0071321-turtle.txt
- N-Triples: 24-1.0071321-rdf-ntriples.txt
- Original Record: 24-1.0071321-source.json
- Full Text
- 24-1.0071321-fulltext.txt
- Citation
- 24-1.0071321.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0071321/manifest