A Sequential FEM-SPH Model of the Heating-Remelting-Cooling of Steel Samples in the Gleeble 3800 Thermo-Mechanical Simulator System

This article presents a sequential model of the heating-remelting-cooling of steel samples based on the finite element method (FEM) and the smoothed particle hydrodynamics (SPH). The numerical implementation of the developed solution was completed as part of the original DEFFEM 3D package, being developed for over ten years, and is a dedicated tool to aid physical simulations performed with modern Gleeble thermo-mechanical simulators. Using the developed DEFFEM 3D software to aid physical simulations allows the number of costly tests to be minimized, and additional process information to be obtained, e


Introduction
In engineering new processes new methodological approaches can be developed thanks to the use of modern Gleeble thermomechanical simulators combined with commercial computer simulation systems [1].On the other hand, original dedicated simulation systems are alternative to commercial simulation systems [2].Already at the engineering stage, or when adapting the existing numerical codes, they allow a solution to be optimally matched to the specific requirements of the user.Also a lower implementation cost, as well as the ease of use of a tool like this, are worth mentioning [2].In solutions of this type, usually the number of parameters controlling the numerical model is limited to the minimum, and the user's interface itself is intuitive and simple to use, even for beginning engineers.It has inspired the author of this study, leading to the development of a concept of the methodology of integrated modelling, combining the advantages of physical simulation and numerical modelling, to aid designing new processes.
The key component of the developed methodological solution are spatial numerical models based on the finite element method (FEM) [2] and the smoothed particle hydrodynamics (SPH) [3,4].
The intensive development of computer technologies, accompanied by the continuous growth in the computing capacity, enables more and more complex physical processes to be effectively modelled.It allows new computing methods to be applied, which could not be used earlier due to extremely long computing times.
A fast development can be observed in the area of the socalled meshfree methods, allowing many physical phenomena to be modelled.They include deformation, solidification [5], or both turbulent [6], and multiphase [7] flows.The SPH method was originally developed to simulate astrophysical phenomena [8].At present, it is still used to describe effects of this type, e.g. a simulation of meteorite fragmentation after entering the atmosphere [9].Thanks to numerous tests, results featuring satisfactory accuracy for complicated problems have been obtained, and it has been found that this method can be extended onto other complex physical issues.This method is characterised by the application of particles, where a computing mesh is not needed to calculate spatial derivatives.Momentum and energy equations are written in the form of ordinary differential equations, which facilitates computing.For instance, the pressure gradient determines the force acting between a couple of particles.In this method, the particle state is described in accordance with the Lagrangian approach, where co-ordinates move together with particles, contrary to the Eulerian approach, where co-ordinates do not change their positions.Modelling flow-related effects is widely used in e.g. the automotive industry, for issues related to designing fuel cells [10], or the high-pressure casting process [11].Also developmental research is carried out on the application of the SPH method in commercial software LS-DYNA for the simulation of flow processes in casting moulds [12].The SPH method has also been applied to simulate a fluid flow in a porous material, where particles constituting the matrix were randomly assigned fixed positions [13].In paper [14] the authors presented the application of the SPH method to the investigation of the diffusion effect in spatial porous materials.The effect of convection in the steady and in the non-stationary state was examined.The findings were compared to other available solutions, and the model developed in this studies was used to compute the diffusion coefficient.Due to problems occurring when using mesh methods (e.g. the finite element method), the SPH method is often applied to model problems, in which large deformations and loss of material cohesion occur.Paper [15] presents the application of the SPH method to simulate the deformation of retaining walls and to describe the landslide shift, where the use of traditional modelling methods turned out to be useless because of large deformations and displacements.The correctness of the proposed approach was verified with tests.Paper [16] presents a material cutting model with the SPH method.The research conducted for a two-dimensional model has confirmed the suitability of this method for correct estimation of shearing forces, which was presented with a few examples of transverse shearing.A similar issue was raised in paper [17].In this case, the coupling of the SPH with the FEM method to model the process of aluminium sheet trimming was presented.Implementing the SPH method in the LS-DYNA program allowed a smooth transition between two domains to be determined without referring to the contact definition.Numerical computing was performed for the spatial strain and stress state.The authors of paper [18] presented the application of the SPH method to simulate the solidification process.The tests were conducted considering the issue of the multi-component mixture solidification.After verifying the findings for a unary case, an attempt of a binary alloy solidification simulation was made.Results of the simulation were verified experimentally.Taking into account the phase transformation and the so-called latent heat released or absorbed during the solidification and melting is a very important issue for the melting and solidification processes.Article [19] introduced and tested various variants of the SPH method concerning simulation of the solidification and heat conduction.The assessment of the new method for various analytical and numerical results confirmed its accuracy, and primarily the ease of numerical implementation.The next area of the SPH method application is the modelling of metal flow and solidification in the high-pressure casting process described in paper [20].In this case, the geometrical complexity and high velocities of the liquid metal, which creates a strongly 3dimensional flow with a substantial free area, fragmentation and splashing, is particularly important.Knowledge of the flow combined with heat transfer and solidification is an important area for this kind of a process.Another example of application of the SPH method for modelling the metal solidification process was presented in paper [21].The solidification of liquid metal in the fast cooling conditions was simulated.The research of the authors confirmed that the SPH method could be successfully used to model the solidification process, and to predict the areas within the sample tested volume where defects in the form of voids or the loss of material continuity caused by shrinkage can occur.The mathematical model presented in the paper will be the basis for the development of DEFFEM 3D software aimed at developing a comprehensive numerical model allows the simulation of deformation of steel in semi solid state.Compared to model based on finite element method (FEM) [22], it allows to eliminate problems with a large mesh distortion during deformation of steel in the semi-solid state (=lack of convergence of optimization algorithms) as well as simultaneous simulation of solidification with local flows of solidifying steel in the mezoscale.

The FEM model
The concept of the proposed sequential solution uses the thermal 3D model based on the FEM [22,23], developed during previous studies.It is based on the non-stationary solution of the Fourier-Kirchhoff differential equation, with an internal voluminal heat source [22,23]: where: τtime.ρdensity, c pspecific heat, λheat transfer coefficient, Qrate of heat generation (the current flow in the sample).
The equation ( 1) was coupled with a functional model of generating heat resulting from the flow of electric current through the sample [22]: where: change in current intensity versus time, change in resistance versus temperature, heating intensity versus time.
The boundary and initial conditions were adapted to the conditions of the Gleeble 3800 simulator module used for the experiments [22].The temperature field within the steel sample volume, obtained as a result of resistance heating is the initial condition in the form of a known temperature field for the solution obtained with the SPH method (simulation of solidification).

The SPH method
The SPH method is a simulation method using interacting particles.A particle in this method is defined as a certain fluid volume, where its physical parameters such as e.g.mass, position in the 3D space, velocity, density or temperature are associated with its centre of mass.In the numerical simulation, the parameters of a system evolving from a certain initial state through a pre-set number of steps (in time) are observed.At each simulation step, the interaction forces between particles are calculated, and equations of motion are solved, therefore, particle positions and velocities are determined at subsequent steps of the simulation.When the model is coupled with the thermal solution also the effects related to heat transfer and conduction can be taken into account.The main idea of the SPH method is the approximation of the field of any physical quantity at selected points in the space.The basic approximating formula is [22,24,25]: where:  is the particle mass,  density, is the kernel function computed for particles with indexes  and , index  corresponds to each neighbouring particle .The quantity ℎ is the so-called smoothing length.Equation (3) describes, how any scalar quantity at point   is calculated on the basis of the values of   from points   .The sum in equation ( 3) passes all particles of the modelled system.The kernel function  has a limited domain, so in equation ( 3), the sum includes only the particles for which the value of function  is non-zero.Many kernel function forms are available, each of them features appropriate properties.On the implementation side, most often splines are used, which enable the trimming radius to be introduced to the simulation.In the presented solution, the cubic spline kernel function was chosen [25]: The gradient of the function  is given by [25]:

Continuity and momentum equation
The SPH method is a discrete version of Navier-Stokes equations.By the digitization of the Navier-Stokes equations, equations for the value of pressure gradient are obtained for the particles, and forces impacting the particles are computed from these equations.Quantities obtained with this method allow the equations of motion to be solved.In the Lagrangian description, the basic equations describing the flow of solidifying metal are [22,24,25]: where:  means time,  velocity,  pressure, viscosity, whereas  external forces acting on the medium.Equation ( 6) is a continuity equation describing changes in density versus time.On the other hand, equation ( 7) is an equation of motion describing the acceleration of the medium analysed.
Using an approximating formula (5), equations ( 6) and (7) in the SPH formalism assume the form [22,24,25]: Viscosity in the SPH method is introduced by an additional term Π ij in the equation of motion (9).In the presented solution, a formulation proposed by Monaghan and Morris was applied [24]: where:  is a parameter introduced to the equation to avoid a singularity in equation ( 10) when (  −   ) approaches zero,  is dynamic viscosity, whereas  = 4.96333 is a constant independent of  [14].Equation ( 9) presents a change in the particle motion as a result of action of pressure, viscosity and forces acting on the fluid.The time step ∆ for integration is limited by the Courant condition [26]: where:  is the sound velocity.

Equation of state
To solve equation (9), it is necessary to compute the pressure on the basis of the adequate equation of state [25]: where: ρ ref is a set reference value of density, while for liquid metal it is assumed that parameter γ = 7.The coefficient β is so selected that relative changes in density are small enough (below 1%).This condition is met for [25,26]: where:   is the maximum velocity.It is assumed that the sound velocity  = 10  , which allows the maximum value to be determined immediately.

Energy equation
The thermal model of heat conduction is based upon the enthalpy method [22,26,27]: where:  is enthalpy,  thermal conductivity,  temperature.Equation ( 14) in the SPH formalism is approximated with the relationship [27]: Equation (15) ensures that heat transfer is automatically continuous across material interfaces.It allows various materials with essentially different specific heats and conductivities to be precisely simulated.
If we assume a constant value of specific heat for the liquid and the solid phase, the relationship between enthalpy and temperature can be presented with the following dependences [26,27]: The solidification heat  describes changes in the amount of energy discharged during particle solidification,   and   are the liquidus and solidus temperatures, respectively,  _ and  _ the specific heat of the liquid and the solid phase, respectively.The temperature of each particle is computed on the basis of the relationships [26,27]: The phase transition (from the liquid to the solid state) will occur when the particle temperature achieves a value lower than the solidus temperature   .In the proposed solution, particles representing the solid phase are in the model treated as a pseudoviscous fluid and move as a result of very high viscosity.The proposed solution allows the mutual forces acting on solid particles and liquid particles to be maintained.The maximum time step ∆ for integration of the energy equation is [26]:

Physical and thermal boundaries
The physical boundary condition is defined as a special class of stationary particles, so-called "dynamic particles" [28], for which equations of motion are not solved in the solution.An approach like this allows a constant position of particles in the space, which does not change over time, to be maintained.Effectively, any shape boundary can be easily simulated.The heat exchange between the liquid/solid and boundary particles occurs by considering the heat conduction.The heat transfer process occurs by cooling down the solution domain through the boundaries.The boundary condition of heat exchange with the environment was defined in the form of heat fluxes, assuming a constant value of substitute heat transfer coefficient  [26].

Test material and thermo-physical data
The material selected for the tests was steel C45 with the chemical composition presented in Table 1 The necessary thermo-physical data for the needs of the numerical simulations were determined with commercial software JMatPro on the basis of chemical composition of the selected steel.Table 2 summarizes the computed thermo-physical properties of steel C45.

Heat conduction and transfer in the adiabatic conditions
The first test verifying the correctness of the numerical implementation of the mathematical model was the simulation of heat conduction and transfer in adiabatic conditions (no heat exchange with the environment occurs).The solution domain with dimensions 1x1x1 m consisted of 80,000 particles.In the simulation, the initial temperature for randomly selected 40,000 particles was assumed 1000°C, for the other 40,000 particles it was 400°C.Fig. 1A presents the generated solution domain with the initial distribution of the temperature field.The following was assumed in the computing: The analysis of the obtained results (Fig. 1B -Fig.1D) indicates that the system approaches the steady-state and the stabilization of the temperature field at a level of 700°C (Fig. 1D).Temperature equalization in the volume of the computational domain was obtained after 57 seconds.Therefore, one can conclude that in terms of quality the developed numerical tool correctly handles mechanisms of heat transfer and conductions.

Computer aid of experiment
For the purpose of preliminary experimental verification of the implemented numerical model, Gleeble 3800 thermomechanical simulator was used.Figure 2 presents a diagram of the adopted configuration of the measurement system.In the experimental tests, cylindrical samples with dimensions Ø10x125mm were used, and they were secured on both sides in copper grips.For the so selected configuration, the sample free zone was about 67 mm.The temperature was measured with two thermocouples.The first thermocouple TC4, at the same time being the control thermocouple, recorded the temperature at the sample surface.On the other hand, the thermocouple marked as TC3 recorded the temperature at the sample core.Based on the determined solidus and liquidus temperatures of the test steel (see Table 2), and bearing in mind that as a result of resistance heating of steel samples within the Gleeble 3800 simulator system the highest temperature values occur within the sample core [22], the following experiment schedule was adopted.The sample was heated to a temperature of 1350°C at a constant heating rate of 20°C/s, followed by heating to a temperature of 1430°C at a rate of 1°C/s.At the subsequent stage, the sample was held at a temperature of 1450°C for 30 seconds to stabilize the temperature within the sample volume.Finally, it was cooled down to the nominal temperature of 1380°C, and after holding for 10 s, cooled down to the nominal temperature of 1368°C.The experiment was made twice in order to obtain a reliable temperature measurement within the sample core.In addition, during the experiments, a special quartz shield, with the primary purpose of protecting the simulator's interior against a leak of liquid metal at the melting process execution stage, was used.The modelling procedure of numerical simulation execution for experiment was performed at three main stages: Stage 1: At the first stage the sample was heated.This was followed by melting and controlled cooling to the test nominal temperature of 1380°C, in accordance with the adopted experiment schedule.A special thermal solution based upon the finite element method was used to assess the temperature field within the sample volume at this stage.The obtained difference between the sample surface temperature of 1380°C, and the temperature of its core of 1420°C (for the first experiment), and 1418°C (for the second experiment) was 40°C and 38°C, respectively.Therefore, referring to the solidus and liquidus temperatures as the reference temperatures, a zone comprising a mixture of the solid and liquid phases formed within the sample volume [22].More details concerning the numerical model of resistance heating, experimental verification, boundary conditions and the heating/remelting procedure itself within the Gleeble 3800 thermo-mechanical simulator system are presented in papers [22,23].
Stage 2: At the second stage of the model procedure, SPH solution domain was generated.Particles are generated within the volume of the sample together with their temperature initialization based on the interpolation temperature from finite element mesh nodal results.
Stage 3: At the third stage, the SPH model developed as part of the study was used to simulate the sample solidification stage in the Gleeble 3800 thermo-mechanical simulator.The adopted simulation time of 0.6 second was equal to the maximum time achieved during the experiment (Fig. 5).Fig. 3 presents the evolution of the mushy zone shape changes within the sample volume during its solidification.The visualization was limited to particles with their temperature higher than the solidus temperature.Fig. 3. Evolution of the mushy zone shape changes within the sample volume during the solidification As time went by, the size of the mushy zone decreased (Fig. 3-A to Fig. 3-K) to achieve the smallest size at the final stage of the solidification process simulation (Fig. 3-L).Analysing the final distribution of the temperature field presented in Fig. 4, one can observe that the obtained maximum temperature of the sample core of 1414.27°C is still higher than the solidus temperature (1412.42°C) of the test steel.Bear in mind that the numerical simulation was made with some simplifications (e.g.constant heat effective coefficients).It did not consider potential local flows of the solidifying steel at the stage of controlled heating/melting and cooling to the test nominal temperature (it stems from current limitations of the developed software DEFFEM 3D).In addition, the numerical algorithm did not include potential grouping of particles followed by their local flow as a group during the solidification, or the forces that occurred as a result of flow of the solidifying steel within the solidifying skeleton of the solid phase.The development of the numerical model with particle grouping algorithms or taking into account forces of interaction between the liquid phase and the solid phase in the mezo-scale should allow results featuring a better compliance with the experiment to be obtained.

Conclusions
The paper presents the main assumptions of a spatial numerical model being a sequential FEM-SPH solution as regards the simulation of heating combined with local remelting of a steel sample, followed by its free cooling in the Gleeble 3800 simulator tool system.The numerical implementation was carried out as part of the DEFFEM 3D software being developed.The developed solution is the basis for the development of an advanced mathematical model based on SPH allowing the simulation of steel deformation in the semi-solid state.The simple numerical tests were performed in order to verify the correctness of numerical implementation (simulation in adiabatic condition).The temperature equalization in the volume of the computational domain was reached after 57 seconds.Using the Gleeble 3800 thermo-mechanical simulator, preliminary experiments allowing additional verification of the mathematical model were carried out.The obtained relative error between the measured and calculated temperature values in the sample core was below 1%.It can therefore be concluded that examples of results presented in this paper show that the adopted model assumptions are correct.

Fig. 4 .
Fig. 4. Final distribution of the temperature field of the sample free zone ( = 0.6 .)

Figure 5
Figure5presents measured and computed changes of temperature in the sample core.During experiments, for test number 1, the sample core solidified in 0.6 seconds.For test number 2, this time was about 0.42 seconds.Analysing the course of temperature change obtained as a result of numerical computing, one can observe that the sample core solidification time is longer than the times obtained during the experiments.Despite these differences, one can find that the obtained results feature good compliance and a similar characteristic of the course.The total relative error between the computed values of temperatures and the values measured during the experiment was below 1% (0.144% and 0.36%, for test numbers 1 and 2, respectively).