Detailed Technical Report by MEFOS

AuthorM Karlberg

Page 275

A4 1 Material of interest

In this project two different materials has been considered, C-Mn-Nb and Cr steels. Where the most relevant measurements are yield stresses versus temperatures and strain rates in the range of interests. One should though point out that if microstructure modelling is of interest, data for these models as well could be important but these data has not been used in this project. However, the first stage in this project was to use the data for modelling of the yield stresses based on dislocation and vacancy theories which means that a lot of work can be done by deformation data since magnitudes of the microstructure parameters generally are known and does not differ so much between different materials. Yield stress data for C-Mn-Nb steel has been collected at the pilot plant in Freiberg where it has been measured for a number of different strain rates at different temperatures. The collected data is especially important for the optimisation of parameters in the existing constitutive model. When parameters have been optimised, validation of the constitutive model for the C-Mn-Nb steel can be performed for compression tests as well as for the benchmark example.

A4 1.1 C-Mn-Nb

Measurements has been sampled at the Swedish institute of material research (SIMR) by compression tests at the strain rates 0.8, 3.5, 10 and 30/s and at temperatures from 800°C up to 1200°C, at every even hundred degrees. All samples have been exposed to a prior austenitisation for 10 hours at 1250°C. This data originally was produced for a similar project at SIMR and therefore was very convenient to use for almost the same objectives in this project.

Data has been processed at the partner TU Freibergs mill facility in Germany where samples are measured at a several temperatures between 600°C up to 1200°C but only temperatures of interests are 800, 850, 900, 1000, 1100 and 1200°C. Strain rates were sampled at 0.1, 1 and 10/s. This kind of data is corresponding to quite low strain rates which typically for the benchmark case is around 20/s. Anyhow the strongest dependence of the strain rate on the yield stress is usually in the lower strain rate regions, lets say between 0.1 up to 10/s whereas this data still gives a lot of valuable information when optimising the parameters in the model.

A4 1.2 Cr-Mo-V steel

Material of interest in the Nordic committee at MEFOS and is categorized as tooling steel. Where the characteristics are high toughness against cyclic temperature changes, high strength at elevated temperatures, easy to harden and machine tooling. Typical yield strengths for this material are about 1400 MPa at room temperatures.

A4 2 Constitutive models

Some of the models presented in this report are the same that were used in a former project, Ref. [A4.1]. The models used in that project were mainly basic equations for description of the evolution of the dislocation density and the vacancy concentration, recovery processes was the same as in this project further explained in Section 2.1. In order to develop andPage 276further increase the complexity of the models a cooperation between SIMR, MIKRAB and MEFOS has been put together. The reason for this was that SIMR together with MIKRAB has great experience and knowledge in microstructure modelling and material science.

A4 2.1 Dislocation density and vacancy formulations

Dislocation density evolution is generally described as the difference between generated mobile dislocations which acts like carrier of the plastic strain and immobilised dislocations. The generation of dislocations is represented by the following relationship:


Where m denotes the Taylor factor, a generalised correcting factor for the critical stress in a slip system correlated to the external applied stress, this parameter is slightly dependent of the strain as also the Burgers vector denoted b is, L is the mean free slip distance for mobile dislocations. Burgers vector usually is considered to be linearly dependent of the temperature by e.g. the following equation:


T is the temperature, C a slope constant and finally b0 the initial burgers vector.

The mean free slip distance is described by Equation (A4.3) as follows.


Where the slip distance L, is inverted proportional against the square root of the dislocation density, !. L0 is the relative level of the slip distance, R0 and kL are proportionality constants needed to be optimised. The mean free slip distance can best be described as the mean distance where a dislocation can move without meeting any obstacles like other dislocations, precipitations or grain boundaries. The recovery process is dependent of several factors where dislocation climbing is one of the dominant ones especially when hot rolling is considered. Recovery of the dislocation density:


Where M is the mobility and c is the vacancy concentration. Since recovery influenced by the climbing process implies diffusion of vacancies or interstitials to the area of climb, this process is a thermal dependent process. This model is only considering point defects like vacancies and not interstitials since the activation energy for interstitials usually are 2-4 times higher than for vacancies, Ref. [A4.4]. The mobility is described as: T


Dm, the vacancy migration times the vacancy concentration is representing the diffusion coefficient. G is the temperature dependent shear modulus, k is the Boltzmann’s constant.

Page 277


By combining Eqs. A4.1-A4.6 this gives us an equation for the derivative of the dislocation density.


A formulation of the vacancy concentration in Eq. (A4.7) is needed since it plays an active role in the recovering part. The vacancy concentration could as well as the dislocation density equations be divided into two different parts; generation and recovery as follows: dt


Deformation of the material produces a considerable amount of vacancies due to among others interactions between dislocations.


C is a constant that needs to be optimised and the rest of the parameters are the same as described in the text above. Recovery is represented by diffusion of vacancies towards areas with high concentration of dislocations.


C is a constant and Dm times c is once again representing the diffusion coefficient for diffusion of vacancies and c0 is the formulation of the equilibrium number of vacancies present in a material, see Eq. (A4.11). This equation is well known and widely used to predict the number of point defects in a material.


Qv is the activation energy for generation of vacancies; R is the universal gas constant and T the temperature. Dm, the vacancy migration used in Eqs. (A4.5), (A4.7) and (A4.10) is represented by Eq. (A4.11).


Dm0 is the initial migration and Qm is the activation energy for migration. Migration could be described as the rate at which a point defect moves from site to site in the lattice.

Page 278

A4 2.2 Recrystallisation

Recovery is a collection of different processes whereas the climbing and cross slip phenomena’s are two basic ones. The maybe most important one or at least a very powerful one is the recrystallisation which has been in the scope of researching for several decades. Since MEFOS by itself not have the deep material insight a cooperation group where put together in the spot of the project. The main objective was to develop the presented constitutive model to cover recrystallisation in multi pass situations and dynamic recrystallisation where the material has a mixed structure of deformed and recrystallised grains. It was started with a recrystallisation model developed at SIMR, Ref. [A4.2]. The model had been used for one-step deformation and static recrystallisation and was programmed in another computer language. It needed to be transformed into FORTRAN in order to fit the code format of the FEM program LS-DYNA. This work was performed but unfortunately it turned out that the model not was applicable at all for dynamic situations. When the situation were discovered the first model had to be reformulated and developed in order to cover situations with mixed grain sizes etc. The present model is based on the previous static recrystallisation model and was developed by SIMR, MIKRAB and MEFOS together. Worth to mention is that the equations in the recrystallisation model is prepared for effects influenced by precipitation.

A4 2.2.1 Subgrain

The whole concept was based on the subgrain structure development in deformed materials. Before we have achieved fully developed grains this stage is proceeded by formation of subgrains and the mean subgrain size can be approximated by Eq. (A4.13):


ksubgrain denotes a constant which has to be optimised and L is the mean free distance of dislocation migration. The growth rate for subgrains is described by:

[VER FORMULA EN PDF ADJUNTO] where the two formulas is used weather there is an ongoing deformation or not. Basically during deformation the growth rate is dependent of the dislocation density and its generation rate in deformed areas. KL is a constant and L is the mean free slip distance as denoted before but in this case it would maybe be represented as the mean free distance between obstacles in the matrix. When no deformation is present the subgrain growth is driven by reduction of subgrain boundary energy. The driving force is represented by Eq. (A4.15), where Fp is the pinning force of particles.

Page 279

A4 2.2.2 Nucleation

The classical nucleation theory can not be adopted for recrystallisation. Instead the description of nucleation in the present model is based on the abnormal subgrain growth mechanism. The critical subgrain size rc, is represented by eq. A4.16 and can be described as the transition size from subgrain embryo into a growing new grain. rc is far larger than the mean subgrain size. This would mean that only abnormal growing subgrains have the possibility to reach the critical size and the nucleation process becomes an abnormal process.


[VER SIMBOLO EN PDF ADJUNTO] is the grain boundary energy [VER SIMBOLO EN PDF ADJUNTO] is the dislocation density in deformed respectively the initial material state, cd and kpg are constants.

The abnormal growing subgrains must satisfy the size and disorientation conditions described further in Ref. [A4.3]. A critical subgrain size, rac is defined as the threshold for abnormal growth:


Abnormal growth rate:




During deformation, subgrains generally not having time for abnormal growth, the new grains are only coming from the subgrains that are larger than the critical size, rc, following relationship hold for the number of grains:


Where [VER SIMBOLO EN PDF ADJUNTO] is total site number and p(r) is the size distribution function of the subgrains.

After the deformation, the nucleation is dependent on the number of abnormal growing subgrains and their growth rate.

[VER FORMULA EN PDF ADJUNTO] rt is the initial size ( at t=0) of the subgrains which reach rc at time t:

Page 280


A4 2.2.3 Growth and coarsening

The final grain size of recrystallised structure depends on both grain growth and coarsening. The growth and coarsening happen concurrently during recrystallisation. Here the growth means the grain growing in the deformed part, and the coarsening means grain coarsening inside the recrystallised part. In order to calculate the growth and coarsening simultaneously, recrystallised grains are treated as two classes – surface grains and interior grains. The surface grains grow to the deformed area. Reduction of storage energy is driving force and grain boundary plays as restrainer, described in Eqs. (A4.23) and (A4.24). The interior grains follow coarsening mechanism and the reduction of grain boundary is the driving force, see Eqs. (A4.25) and (A4.26). The grain growing rate combining growth and coarsening is given in Eq. (A4.27).


Where Ninte and Nsurf are grain numbers of interior and surface grains.

A4 2.2.4 Dynamic and Multi-step deformation

After the onset of recrystallisation, the material is separated into two parts: a recrystallised and non-recrystallised matrix. In case of dynamic recrystallisation, the deformation continues in both parts, but at different strain rates. The recrystallised (soft) part has a higher strain rate and the already deformed part has a lower strain rate. The grain sizes in these two parts are given as Rrec and Rdef. The global strain rate is given as [VER SIMBOLO EN PDF ADJUNTO]. The strain rate in these two parts can be calculated from Eqs. (A4.28) and (A4.29). In case of multi- step deformation, if the recrystallisation is not complete during the pause time, the situation is similar as dynamic recrystallisation case. The recrystallised and non-recrystallised materials will be deformed according to Eqs. (A4.28) and (A4.29). If the recrystallisation isPage 281complete during pause time, then the deformation can be treated simply in the same way as if it was a one step deformation.




A4 2.2.5 Flow stress

The influence of microstructure at the flow stress can be well described by:


[VER SIMBOLO EN PD ADJUNTO] is the initial frictional stress generally dependent of the temperature and the fraction of precipitations in the matrix, [VER SIMBOLO EN PDF ADJUNTO] is the mean dislocation density in one cell for partly recrystallised structures and the dislocation density for fully deformed or re-crystallised structures. [VER SIMBOLO EN PDF ADJUNTO] is a parameter used to describe the strength of obstacles, this parameter is treated as a constant independent of e.g. temperature and strain rate etc. The parameters m, G and b is the Taylor factor, shear modulus and Burgers vector respectively calculated and described as in the proceeding text. These parameters are slightly dependent of temperature and the strain level.

Page 282

A4 3 Optimisation

Optimising in order to find extreme points of a function or systems of functions is a topic that by itself can be an objective of research. Optimisation generally becomes very complicated when the number of degree of freedoms are increasing and when the functions to be optimised are coupled and strongly dependent of each others. When optimising, maybe the hardest thing is to know weather or not the global minima is reached or if there is a local minima that has been found in the error function. The basic idea however here is to create a function F which for each set of the vector x returns a scalar which can be minimized regarding to a representative norm.


F is a function dependent of sub-functions, F is continuous and has continuous derivatives in the n dimensional space. This is essential since optimisation algorithms usually need calculations of the first and second order derivative. The functions f1 to fn are the set of first order ordinary differential equations representing the constitutive model presented in the proceeding section. The vector x are forming a multidimensional space where all elements are representing one axis and where the minima should be found. A lot of efforts have been put on the optimisation part since the parameters are important for the ability to predict the final flow stress. It is impossible to present all the attempts made within the project except for some essential steps used in the optimisation process.

A4 3.1 Screening

Screening means that the measured data which should be used for the optimisation has to be checked and/or cleaned from irrelevant samples which always are present. Almost in all cases there is a need for screening the data in order to extract as much information as it is possible from it. There are always scattering in measurement devices due to vibrations, human factor and samples are performed in another reason etc. where all factors, called secondary factors are influencing the result in an order that not are represented by physical reasons rather then by randomised factors.

A4 3.1.1 Oscillations

In this project the C-Mn-Nb steel data were gathered in the frame of this project but needed filtering algorithms in order to stabilize the strain rate which was oscillating heavily. When using the data in models which are directly dependent of the strain rate like the constitutive models in chapter 0 one is obtaining troubles if the amplitude in the strain rate is too high. Once again data is data and this amplitude might be realistic but for modelling purposes it can be more favourable e.g. to use the local mean value instead.

For this reason a method called the “Loess” method, Ref. [A4.5] were used in the smoothing process of the data. "Loess" are derived from the term "locally weighted scatter plot smooth" and is locally approximating data by linear regression, samples in the neighbourhood to the one to be approximated is weighted which is the width of the filter. However in Fig. A4.1, smoothed strain rate data is presented which clearly are pointing out the effectiveness of the filtering algorithm. Usually ten up to about fifty samples are used in the local regression model.

The problem with an oscillating strain rate is that when integrating the constitutive model with a specific time increment different from the sampled one. The constitutive models arePage 283sensitive at the strain rate level whereas strong oscillation can influence the optimisation result.

A4 3.1.2 Interpolation

Data for the Cr steel was as mentioned in the Section 0 not produced for this project from the first time. Hence the number of samples was too few and it is not possible to see what history the data has, Fig. A4.2, if they have been smoothed, if the strain rate is corrected etc. However the easiest and most cost efficient way to expand the dataset is to interpolate new samples based on the original samples with the drawback that information in the extra samples are constructed and not real. The shape of the flow stress curve is though almost continuous and smooth so there is no major risk by using interpolation. In Fig. A4.2, the marked curve is representing the interpolated data which seems to behave in a reasonable way except for the end of the curve where the curvature seems to be slightly too high.

A4 3.2 Normalizing

When minimising F in Eq. (A4.55) by an algorithms that uses the gradient calculations it is logical to realise that if two directions or elements of the input vector x have magnitudes of 1015 and 10-5 like the dislocation density respectively the mean free slip distance, there is a risk that the error in the derivative calculations for one parameter is larger than the absolute value of another one. Therefore it could be necessary to normalise the input vector. An easy way to achieve this is to initially use an internal parameter vector usually equal to the unit vector outside the algorithm, inside the algorithm the elements of the parameter vector x, is transformed to accurate magnitudes by element wise multiplication with vector m. This means that:


Where n is the number of parameters to be optimised and xi is the “scaled” parameter vector.

A4 3.3 Norm and optimisation criteria

As important as the rest of the optimisation is the level where one could say that the optimisation is successful and should be terminated. A convenient and very common error estimation to use is the variant of the L2 norm called the RSM (Root Mean Squared) error or simply RSME:


[VER SIMBOLO EN PDF ADJUNTO] denotes the measured respective the predicted sample and n the number of samples. An advantage of using the RMSE is that it is so widely used that it easily can be related to similar work or that one is familiar with it.

Typical RMSE for optimisation of one strain rate at one temperature use to be in-between the range 0.5-1.5 MPa. For optimisation of one set of strain rates at one temperature the error has increased to about 4-10 MPa. When optimising a set of yield stress curves with various strain rates at one temperature the error for each of the strain rates is use to be approximately the same. When optimising both temperatures and strain rates the error is generally increasing since it becomes more difficult to find the true global minima. UsuallyPage 284the optimisation error is dominated and very sensible by the predictions at lower strains basically because of two reasons. First is of a direct kind which means that the prediction error is large when the derivative is high and the second one is that if the initial derivative are wrongly predicted this influents the rest of the integration of the constitutive models since there are differential equations.

A4 3.4 Algorithms, strategies and results

There are a numerous number of different optimising algorithms widely used and available in both literature and as computer code. Which one who will produce the best optimising result is generally a very hard question to answer. And in many cases it will need a great knowledge about functional behaviour and algorithms in order to choose the optimal one. Generally a trail and error procedure is the only way to find out the best method for the specific problem. A broad categorisation of methods would be, methods that only use function evaluations like the Nelder-Mead’s Simplex method, gradient based methods like “Steepest descent” and higher order derivative methods such as the Newton’s method. Newton based algorithms usually need the Hessian matrix to be calculated which is very time- and memory consuming process. Quasi-Newton methods are a modification of these and are avoiding this drawback by updating the Hessian matrix and are generally much more efficient, example of algorithms are the well known Levenberg-Marquardts method and BFGS, Ref. [A4.6]. What can be said though is that the RAM memory usage differs a lot and should be a criterion itself when choosing an algorithm. For optimisation of the parameters in these constitutive models all kind of optimisation algorithms has been investigated but the Levenberg-Marquardts method with a line search method has generally produced the best results.

What was discussed above is one aspect of the complete optimising problem, which however not is fully highlighted yet since the way in applying the algorithms has its own important role and can be varied in numerous ways.

The stress level in the flow stress data is dependent of three major variables, the strain rate, the temperature and the strain level. In Section 2 there were discussed about recovery mechanisms such as climbing and cross slip. Climbing is one of the processes that mainly are dependent of the temperature instead of the strain rate since the way a dislocation is climbing is due to thermal activated diffusion of a vacancies or interstitials. This is pointing out that the optimisation itself could be separated to cover the temperature or/and the strain rate dependentness each by each or all together. In the project a lot of different strategies have been applied where among the most successful ones has been to first optimise for a medium temperature let say 900°C. This is tuning the constants into a reasonable level whereas the optimisation can be extended with strain rates at another temperature or temperatures, assumed of course that the first optimisation were successful.

Since a material response are dependent of many primary and secondary effects which has their strongest significance at different temperatures, strains or strain rates it is generally difficult to optimise the parameters for all situations. When the optimised parameters have been used in FEM simulations usually the parameters have been optimised for two or three temperatures and for all strain rates. The important thing then is that the temperature span in the optimisation is covering the span that will be present in the simulations. One should though remember that two primary effects that have a great significance are recrystallisation and precipitation which not have been included in the optimisations.

Page 285

A4 3.5 Optimisation of Cr steels
A4 3.5.1 Polynomial prediction of the initial frictional stress

For all flow stress data, the most difficult part to optimise except for the recrystallisation part of the curve is the one at strains lower then 0.05. One simple reason is that there was no information available for this particular Cr data, see Section 0. Other indirect reasons are that there are high derivatives in this area and that the dislocation density is increasing so fast that even small deviations in optimised parameters have great effects at the shape of the rest of the curve.

To overcome the missing of samples in the low strain regions, optimisation were just applied to the data that existed at strains greater then 0.05. A polynomial function, Eq. (A4.58), was fitted to the sampled stress levels at strains equal to 0.05 for all temperatures and strain rates:


For stress calculations at strains lower then 0.05 a linearly interpolated stress level based on the predicted one could be used, which is shown in Fig. A4.3. The polynomial predictions are very accurate relatively and coefficients in the polynomial function were easily fitted to the measured stresses at strains equal to 0.05. In Fig. A4.4 this method is shown whereas predictions at different temperatures and strain rates are plotted in the same figure. The black curve is showing a prediction at 750 C which is an extrapolation of the constitutive model, the red curves are interpolations and the blue ones are predictions at data that the model already has been exposed to. The figure is showing that the curves not have continuous derivatives which can have effect on the convergence of FEM solution. The prediction at 750°C are actually not relevant since the temperature is to low but it is very important to expose the model for extrapolations in order to get a feeling for the robustness of it. In Fig. A4.5, the extrapolation of the strain rate is plotted. Also in this figure there is no strange behaviour of the curve outside of the strain rate span used in the optimisation.

A4 3.5.2 Temperature dependent initial frictional stress

The method applied in the proceeding section is possible to implement, it gives reasonable results and the major advantage is that it easy to find a global minima that cover all temperatures and strain rates with the same parameters. The drawback though is that it is not accurate in a physical meaning to assume that all curves should have the same dislocation density at strains equal to 0.05 but different stresses. In fact the only or almost the only variable of greater significance in the flow stress equation, Section 0 is the dislocation density. A more physically accurate representation then is to assume that the initial frictional stress at strains equal to zero is a thermally dependent function like:


Since data for strains lover then 0.05 not was available for the Cr steel, two alternatives were possible. First, one could approximate this part quite reasonable by just construct samples that are maintaining the shape of the curves. Secondly one could produce a sampling sequence where these parts are measured. The first alternative was adopted by financial and logical reasons. However, by extending the curves to zero strains and approximating a reasonable initial frictional stress the optimisation becomes more difficult. In Fig. A4.6 predictions for strain rates at 800° is made where the RMSE is approximately 1.63 MPa. When optimising for several temperatures the RMSE generally increases, oncePage 286again because the models has to cover a wider range of effects. In Fig. A4.7 and Fig. A4.8 optimising results for two and three different temperatures are presented. RMSE’s in the figure are 4.2 MPa respectively 6.78 MPa.

A4 3.6 Optimisation of C-Mn-Nb steels

The flow stress measurements of the C-Mn-Nb steels were performed by TU Freiberg. For optimisation, flow stresses at three strain rates, 0.1, 1 and 10/s was logged and has been used in the optimisation together with five different temperatures 800, 850, 900, 1000, 1100 and 1200 °C. The strain rates does not cover the full range of strain rates present in the benchmark rolling case but it cover the most important region where the flow stress behaviour changing at most. If the optimised parameters are physically reasonable the models have showed a stable behaviour when extrapolating.

Notable for optimisation of these data is that no use of the polynomial prediction of the initial frictional stress has been used. This is due to that data for lower strains are present in the measurements and that it is not a physically accurate interpretation.

The optimisation for the C-Mn-Nb data has been much more difficult than optimisation of Cr steel data. There is unclear what the reason for this is but one reason might be that the highest strain rates sometimes at several temperatures are showing a peculiar shape at low strains which could come from troubles with the measurement device. Another reason could be the presence of particles but this would be case for the Cr steel for the Cr steels as well. However, different optimisation techniques and algorithms have been used where the result has been approximately the same. Figures A4.9 to A4.13 are showing results for optimisation at temperatures from 850°C up to 1200°C where the typical RMSE lies around 4-8 MPa. When looking at Figs. A4.9 to A4.13 it is obvious that the dynamic recrystallisation not is represented by any equations in the constitutive model. Since the recrystallisation model not has been fully implemented into the rest of the constitutive models the optimisation should just concern fitting of parameters in the vacancy concentration- and dislocation generation equations. Figure A4.14 is showing the result for optimisation of flow stresses at two temperatures and at three strain rates.

A4 4 Implementation

In the former project, Ref. [A4.1] a subroutine were also implemented into LS-DYNA for applying the constitutive model to FEM simulations. A lot of things has though been changed since then e.g. that the FEM program has changed, the constitutive model has been developed in some areas, and optimisation is improved.

The implementation has been compiled only with object files for the single processor version which means that the solution times could be decreased if they were a problem, for small simulations like the one made in this project this however does not matters.

A4 4.1 Explanation of subroutine

In Section 0, where the constitutive model is presented it can be seen that the governing microstructure equations are a set of first order ordinary differential equations. To solve these they must be integrated in respect of time. This can be achieved by some advanced algorithm like e.g. the well known Runge-Kutta method or if the time increment is very short by a direct approximation like the Euler method where the derivative is multiplied by the time increment as in Eq. (A4.41).

Page 287


Since LS-DYNA is based on an explicit FEM method which uses very short time increments it is very convenient to use fast Euler integration schemes for the subroutine without reduced accuracy in the calculations.

Flowchart A4.1: schematic figure for implementation of the UMAT (User Defined Material Subroutine)


In Flowchart A4.1, a schematic description of the implemented UMAT subroutine is described. It works like the subroutine is called at every time increment with the time, strain,Page 288old stresses, time increment history variables and material specific constants as arguments. In the beginning of the program a subroutine to the main program is initialising history variables and material specific parameters and at the first time increment the initial values also has to be calculated, marked as sequence 1 in Flowchart A4.1.

Based on calculated incremental deformations and strains in the main program first step in the UMAT is to update the stress tensor for an assumed fully elastic deformation.

A4 4.1.1 Basic equations

Dependent of what kind of FEM formulations that are used the incompressibility is treated in different ways. In LS-DYNA the incremental stress is represented by Eq. (A4.42).


Where [VER SIMBOLO EN PDF ADJUNTO] is the Kronecker delta and [VER SIMBOLO EN PDF ADJUNTO] is the bulk modulus


G is the shear modulus:

[VER FORMULA EN PDF ADJUNTO] and [VER SIMBOLO EN PDF ADJUNTO] in Eq. (A4.42) is the deviatoric strain components used in order to neglect effects influenced by the hydrostatic stress.


And finally the [VER SIMBOLO EN PDF ADJUNTO] is the volumetric strain which is anticipated to be small in comparison with the deviatoric strains since the material in scope are almost incompressible.


By Eqs. (A4.42) to (A4.46) and the stress at tn-1 the updated stress is represented by following relationship:


The updated stress in Eq. (A4.47) is a first elastic trial where it is assumed that all deformations are producing only elastic responses.

Updating the incremental effective strain rate is performed by Eq. (A4.48):


Page 289

A4 Plasticity

In order to be able to describe the plastic behaviour of an elastoplastic material when it is deformed, it is essential to have a representative yield function of the form.


Where the dots do denotes material specific state variables like the dislocation density or any other variable describing the microstructure which influences the yield surface. The material response will be pure elastic as long as the yield function is lower then zero.

[VER FORMULA EN PDF ADJUNTO] and elastic or plastic when the function is equal to zero.


Notable is that positive values of the yield function are inadmissible. The function f is usually implemented as


Where [VER SIMBOLO EN PDF ADJUNTO] is the effective elastic stress according to von Mises plasticity theory and [VER SIMBOLO EN PDF ADJUNTO] is the yield stress calculated by the constitutive model.


[VER SIMBOLO EN PDF ADJUNTO] is the deviatoric stress tensor calculated by


In order to check if yielding occurs or not the only step that last in the subroutine is to calculate the yield stress which is done by the constitutive models and Eq. (A4.35), the modified Hirsch equation. This is marked as sequence 2 in Flowchart A4.1. If the statement in Eq. (A4.50) is true the program will adopt the elastic trial calculations and save the incremental stresses as history variables. If the statement in Eq. (A4.50) not is fulfilled the material response will be plastic and/or elastic. The yield function, Eq. (A4.51) still is the governing equation and the overshooting stress has to be transformed back to the yield surface by a penalty algorithm. There are several ones to use with a different level of complexity, in this subroutine however the well known Radial return method, Ref. [A4.7], were implemented because of its simplicity and reasonable accuracy.

Page 290

A4 Radial return

The purpose with the algorithm is to radially scale the overshooting deviatoric stresses back to the flow surface defined by Eqs. (A4.51) and (A4.52). Routines like the Radial return method is not automatically implemented in LS-DYNA’s source code, at least there is no easy way to call existing routines and this has to be manually implemented by the user. First step is to convert the overshooting stress to a plastic strain increment which is done by Eq. (A4.54).


The plastic strain contribution has to be added and stored into a specific history variable devoted for the effective plastic strain. Next step is to calculate a penalty factor lambda to scale the deviatoric stress tensor back to flow surface.


The last step is to subtract the converted plastic strain from the elastic stress case


A4 4.2 Troubleshooting of compiling

In LS-DYNA the common way of extending the code with your own subroutines, like a new material subroutine is to write your code in FORTRAN and compile/link it together with the compiled source code in form of object files delivered by the software supplier. It could be a lot of troubleshooting with this since there is no or at least a very poor explanation of what parameters that should be used or what parameter that are public etc. Another thing that can make this hard to work with is that there are no or limited possibilities to debug the code.

In later part of 2004 the computer used for simulations at MEFOS were replaced by a PC cluster, this implied that LS-DYNA had to be updated and reinstalled in order to run on the new platform which was different from the former one. Since there not exists any options to compile and link subroutine and object files in a running mode, this caused big troubles to find a working compiler. The problem was to find a compiler that could run on the desired platform compile and link the subroutine with the object files libraries delivered by the supplier. Usually the suppliers use to have a number of versions of object files compiled with different compilers in order to be suitable for the customer’s system requirements. In this case MEFOS could not find a compiler that could run on MEFOS platform and in the same time be compatible with the delivered object libraries. Finally, by contact with the Livermore Software Technology Corporation in USA one solution were find and could be installed at our system. This did though delay the planned work and unfortunately brought several weeks of extra work at problems not significant for the project.

Page 291

A4 5 FEM Simulation
A4 5.1 Introduction

Two different simulation models have been set up during the project, the first one is a model for simulating compressions of circular test specimens and the second one is a model for simulating hot rolling. In the beginning of the project a benchmark example was worked out by the partners in order to have one reference point when discussing models etc. The benchmark model is based on geometries in the pilot plant at TU Freibergs workshop. Since MEFOS did not have any mill data available in the beginning of the project, the rolling model were just implemented but not used whereas the compression model was created for validation of the constitutive model against the measured data.

Major options that the FEM models were going to handle were:

Thermo mechanical simulations: this includes radiation, convection and conduction

Pause time simulation

Reasonable short simulation times

Contact, both thermal and mechanical

Both models, compression and benchmark, has been used for both thermo mechanical and only structural simulations. There is though an increase in simulation time when solving the relevant heat equations. Pause time simulations were not fully developed since it showed to be a quite time consuming work. The reason for this was that LS-DYNA is based on explicit FEM and explicit solutions is very seldom longer then a couple of seconds since the time steps are very short which means that the simulation times becomes very long if the pause times considered is about 5-10 minutes. This is though dependent of the number of elements, single or multi processor mode, or in what extent the simulation could be mass scaled.

Generally the thermo mechanical simulation possibilities have been improved in LS-DYNA but are still not flexible enough. E.g. when simulating pause times it would be very convenient to manually extend the time step and “turn off” the structural solver and only use the thermal solver which have longer time increments. This is impossible today when a user defined subroutine is compiled to the code.

When the rolling model were created, thermal conduction in use with the contact algorithm was very inefficient in searching for elements in contact. Due to this, the conduction option was put to a side and the most of the thermal modelling was performed with radiative and convective boundary conditions. This has effects on the solutions since it usually is a temperature difference between the roll and work piece of 8-900°C but it was chosen that the increased simulation time not was worth the small improvement in the thermal solution that is present since the simulation times are so short.

Page 292

A4 5.2 Compression model
A4 5.2.1 Geometry, Boundary conditions and specification of simulation

The compression model was set up since it was convenient to be able to make fast simulations for validation of relatively simple cases and geometries. The compression models also are reflecting the way the data has been sampled whereas this is a good tool to make sure that models are implemented correctly and that the optimisation is accurate. For a long time there were only data available for compression tests and therefore it was needed to create a tool for validation of these particular yield stress data.

The model has been used for both mechanically and thermo-mechanically coupled simulations and is simulating the constant or varying strain rate compressions of circular specimens up to strains equal to 1, see Fig. A4.15.

In order to reduce the number of nodes and elements symmetry is used whereas only one eighth of the original specimen is modelled. Nodes on the bottom surface have been constrained to move in the axial or y-direction and are free to move in the x and z directions. Nodes on the inward surfaces are constrained in the x and y directions respectively in order to maintain the position of the cut. Boundary conditions for the top surfaces have been modelled in several different ways depending if contact should be included or not.

With contact included the stress state will act in a slightly different way since the friction is “holding” the upper surface. There is no easy work to predict an accurate magnitude of the friction coefficient which not has been of priority within this work where a coulomb friction model has been used where the friction coefficient typically has been put to 0.3. The forging dies have been modelled as rigid walls which mean that the dies not are deformable which also are affecting the results.

When contact not has been used, the set of nodes forming the upper surface has instead been given a prescribed motion and are free to move in the x-y plane. This is more easily to implement and there is less parameters that has to be prescribed for e.g. friction and penetration parameters. Without contact no un-symmetric thickening effect takes place which has influence at the effective stress state during the compression. This is however considered as minor effects and could be considered in order to improve predictions.

A4 5.2.2 Results and discussion

Simulations with the compression model are very fast even though there is a physically based material subroutine that has to be integrated at every time increment. For compression simulations up to a strain equal to 1 and at 1/s in strain rate, the simulation time is about ten minutes and for lower or higher strain rates the simulation time increases respectively decreases. The simulation time is dependent of the mass scaling, which in all simulation has been about ten to thirty times the original density. Validation of the model has been done at several temperatures and in Figs. A4.17 to A4.21 predictions of the effective yield stresses for the Cr steel at 900°C and at different strain rates are presented. The figures are presenting three different curves, the blue curve is showing the measured data, the red dashed one is the resulting yield stress from the FEM simulation and finally the black curve is representing the resultant curve when optimising the constitutive model. Figures A4.21 and A4.22 are showing the same kind of results but for the C-Mn-Nb. It can be seen that the results are accurate but with some deviances where the yield stress curves are flattening out or where the change in derivative is largest. In Figs. A4.16, A4.17 and A4.18 recrystallisation are present as can be seen at the dramatic softening of the curvesPage 293around strains equal to 0.4. In this figures it becomes clear that no constitutive model are dealing with this phenomena since the predictions are still maintaining an almost constant stress level. The predictions of the yield stresses of the Cr-Mo-V steel generally looks better but that might be effects from the construction of the initial frictional stress that might be too high whereas the changes in derivative are easier to predict.

In Figs. A4.22 to A4.25 results for the history variables for C-Mn-Nb steel are plotted. These figures are showing a simulation where contact is included and the resultant effective stress and effective plastic strain has both their maximum at the centre of the specimen. In this figures the effective plastic strain level is approximately 0.5 which is half of the strain maximum in the simulation. The dislocation density, plotted in Fig. A4.22 has at this stress level increased from an initial level of 1013 /m2 up to a level in the magnitude of 1015/m2. If this are reasonable or not is hard to answer since it depends on different factors and must be validated by itself in order to ensure that this is the case. This is of course a risk when validation of the prediction ability of the constitutive model is based on just one parameter like the yield stress as in this case. The prediction of the yield stress can be very good whereas the microstructure parameters like the dislocation density not are physical at all. One way to get around this problem is to put constraints at some parameters when optimising the model and relate the constraints on measurements like for instance measurements of the dislocation density. However, Fig. A4.22 is showing that the dislocation goes hand in hand with the effective stress level as it should do with Eq. (A4.34) in mind. In Fig. A4.25 the vacancy concentration is plotted.

In Section 0 it was mentioned that one attempt to get around the missing of initial data was to predict the initial frictional stress and it was mentioned that a problem was that small deviations in the strain rate is affecting the prediction of the [VER SIMBOLO EN PDF ADJUNTO] 0.05 stress level. In Fig. A4.26 the drawback with this formulation is obvious for a FEM simulation for the Cr-Mo-V steel.

A4 5.3 Rolling model
A4 5.3.1 Geometry, Boundary conditions and specification of simulation

In order to be able to compare the different partners a unified model were set up as a benchmark example as mentioned earlier. This benchmark example is a hot rolling case of a work piece with geometrical and technical prescriptions. The same case was also set up in TU Freiberg pilot plant and milling was performed in pilot scale at that environment.

This model has as the compression model also the thermo mechanical coupling condition available where simulations have been performed with both options. The constitutive model always needs a temperature but this can be prescribed without any solution of the temperature. When thermo mechanically coupled simulations has been performed the heat fluxes from the material to the surrounding takes place by radiation and convection. The conductive heat loss from the material to the work roll is the major surface heat flux together with the radiation in rolling situations. There is un- fortunately no possibility to use the conductive contact option when one is using contact entities as models of the work rolls which means that one has to model it as rigid body. The problem formulation that occurs is that if the extended simulation time due to the increased number of elements is worth the slightly more accurate temperature calculation. For surface conditions the answer would be yes but for the overall result for the material behaviour it is discussible and was discarded in this work.

As was clarified in the discussion above was that the work roll has been modelled as a contact entity which means that it is described as a mathematical surface used for definingPage 294contact, they are not considered as deformable bodies whereas no flattening of the rolls are considered in the model.

For simplicity reasons the model is split up into four quarters where just one are modelled because of the symmetry. Figure A4.27 is showing the geometric model where the small triad is indicating the positive directions of the model. All surfaces with a negative normal at the work piece are modelled as symmetry surfaces. The upper work roll which is the one that has been included due to the symmetry conditions are showed where it is obvious that it actually does not consists of any elements. The major advantage by using the entity option is that it decreases the simulation times due to the lower number of elements.

The friction force between the roll and the work piece is modelled by a coulomb friction model and no deformation of the roll takes place.

Figure A4.29 is showing the calculated yield stress or the effective stress which is seen to have its maximum somewhere in the deformation zone produced by the work roll. Relaxation becomes present in the positive x direction which is indicated as blue coloured areas and effects of the higher heat fluxes at the edges and corners could clearly be seen. In Fig. A4.32 the flow stress is plotted for two nodes where one is located in the x wise midsection and the other one is located in the front section of the work piece. The red curve is showing the stress level in the mid section and the stress is increasing from the start time even when this element not has been exposed to any deformation yet. The reason why the dislocation density and flow stress is increasing is that the material is experiencing an increased strain rate level because of the deformation introduced shock waves that are travelling forth and back through the material, see Fig. A4.34. The waves which are a result of the first contact between the work piece and the roll are only producing elastic distortions but they still have effects on the stress state whereas the flow stress are increasing onto it reaches the roll gap when time is approximately 0.17 s. When the work piece has entered the roll gap, flow stress is dramatically increasing up to its maximum level at approximately 190 MPa and the following relaxation is reducing the stress by 15% to approximately 160 MPa. In Fig. A4.34 the effective strain rate for one element in the mid section of the plate is plotted, it is showing that this mill operation are producing strain rates up to approximately 18/s.

A4 5.3.2 Results and discussion

A number of different simulations have been performed within the project. But for the benchmark rolling simulations only the C-Mn-Nb steel has been of interest since this is the actual connection to the other partners and the natural point for discussing and comparing the models.

In Figs. A4.27 to A4.31 results from one of the simulations are shown. The simulation times in the figures are 0.3 seconds which means that the pass not are completed yet and in the rear part of the work piece there are still an ongoing deformation.

In Fig. A4.28 the effective plastic strain are plotted, and there are some deviations in the work piece. Since the effective plastic strain is considered, no direct direction dependentness could be detected. The corner and the edge zones are experiencing a slightly lower effective strain because it not is increasing as much in length as the centre zone. In Fig. A4.33 the effective plastic strains are presented at the two elements as the flow stress presented for. Notable are that the strain level are slightly higher in the mid then in the front section.

Page 295

Figure A4.31 are showing the calculated temperature profile and since the rolling time is only 0.35 seconds there are not enough time for any larger temperatures changes. The small temperature changing that is occurring is that all outward surfaces are slightly cooled and the corners and edges are as expected the zones with the highest heat loss showing a decreased temperature which influences the microstructure and resulting in a higher flow stress then the warmer kernel of the plate. In this figure heat conduction between the work roll and the test piece to the roll has not been included in the simulation.

Figure A4.30 is showing the vacancy concentration which remains almost constant when material has passed the deformation zone.

A4 6 Discussions

Physical models for describing the microstructure connected to material properties is an expanding part of the modelling area. Reason for this is that results can be used for a more accurate prediction of materials mechanical properties like the yield stress, fracture stress and elongation which can be used in process control situations. The underlying physical theories describing material behaviour like dislocation, recrystallisation and recovery is outlined for a several decades ago but of practical reasons not applied to reality before the computer science made this possible during the last decade. Today it is possible to run extensive and complicated models describing at lots of phenomena’s without any increase in computational times compared to a situation for 10-15 years ago.

The perfect constitutive model is a model that can predict the material response under environmental influences like the temperature, deformation, chemical components etc and do perfect predictions of the microstructure which in turn gives the relevant mechanical property like the yield stress. This will of course never be the situation but one can continuously refine the prediction ability by adding reasonable physically based equation to the constitutive model. The model is today covering microstructure properties like the dislocation density and the vacancy concentration and their active role in the recovery and generation processes. Recoveries in the model are based on classical theories for thermal/diffusion controlled climbing and glide processes. Two major parts are still missing in the model which is the description of precipitations and recrystallisation. Since precipitation models generally need information about the chemical equilibrium in the material it will imply that these are calculated or predicted. During the project it were decided that it not was reasonable to put efforts in this work within this project. The recrystallisation model is developed in cooperation between MEFOS, SIMR and MIKRAB and are working for dynamic recrystallisation situations and when the matrix consists of mixed grain size structures; deformed and undeformed. The recrystallisation model are integrated with the original constitutive model and can be used for predictions but has not been implemented into LS-DYNA because a considerable time had to spent on compilation trouble and that the recrystallisation model had to be reformulated.

The constitutive model is working and seems to handle the most relevant microstructure properties and their evolutions but this necessarily does not mean that it predicts all microstructure variables as good as it predicts the flow stress. To be sure that this is the case one has to verify e.g. the calculated dislocation densities with measured ones, this is not made within this project but by experience from literature and earlier projects reasonable validation still can be made.

One very important part of the project is the optimising of the parameters in the constitutive model which also has been very time consuming. If we consider that the equations in the models are fully describing the physics of the microstructure and that there are noPage 296measurement errors in the data present it would be easier to find a global minimum of the optimisation. That is not the situation and there is very hard to know if the global minimum is reached or if optimisation could be further improved since there will always be secondary effects that have influence in the final solution. Another aspect of the optimisation is that there are hard to know if the constitutive model is holding enough degrees of freedoms. If the constitutive model has too many degrees of freedoms the risk is that the model will be extremely dangerous to extrapolate outside the “input space” with respect to strain rate and temperatures. On the other hand, if the degrees of freedoms are too less the predictions will not be able to follow curvatures as well, which can be thought of a straight line which has a low degree of freedom. Dependent on in what situation the constitutive model will be used the input space in temperature for the optimisation should be in proportion to the problem temperature. This is a balance between the prediction ability which increases for a narrow temperature span and the temperature span needed in the simulations.

The implementation of the UMAT was made with help from the old project but there were however some changes made since the old UMAT were implemented in another version and on another computer platform. There is though not a straight forward work to implement a UMAT into LS-DYNA and without support there is impossible to know the name of different variables needed for calculations of the yield stresses. One could also expect that in order to handle plasticity, radial return would be an already programmed subroutine that could be called which is probably is but unfortunately not is public. However, the implementation of recrystallisation model was delayed partly because the change of computer system at MEFOS where the compiler and operative system had to be changed and caused problems in finding a combination that were working our tasks. Another reason for delay was that the recrystallisation model had to be developed and reprogrammed since the original one could not handle the mixed grain size structures in a reasonable way.

When the complete constitutive model was implemented, two kinds of deformation test cases were implemented; one for simulating a compression case and the other one for a benchmark rolling case. Simulations have been performed in different complexities with contact, thermal coupling options enabled in some cases. The main objective though is to check that the constitutive model is working reasonable and not to put a lot of efforts in improving for e.g. the friction or some other simulation specific parameter. The objectives of the simulations have been that they should be less detailed but detailed enough to catch the overall problem formulation. Results from simulation are showing a good agreement with measured data and when formulations of the initial frictional stress were reformulated the convergence in the solution procedure became even more stable.

A4 7 References

A4.1 Ingham, P M: 'The effect of strain reversal and strain-time path on constitutive relationships for metal rolling/forming processes', Final Report, ECSC Project, 7210.EC/811, 2000, ISBN-92-894-1546-0, pp99-139.

A4.2 A physical model for prediction of microstructure evolution during thermomechanical processing, Wang, X; Siwecki, T; Engberg, G, Materials Science Forum. 2003, Vol. 426-432, pp3801-3806.

A4.3 Humphreys F.J., Acta mater., 45, (1997), 4231-4240.

A4.4 Theory of dislocations, Hirth, J. P; Lothe J, second edition, ISBN 0-471-09125-1.

Page 297

A4.5 Robust locally-weighted regression and smoothing scatterplots, Cleveland, W. S. (1979), Journal of the American Statistical Association, 74, 829-836.

A4.6 Optimisation toolbox, User’s guide, version 2, Coleman, T; Branch, M.A; Grace, A.

A4.7 Calculation of Elastic-Plastic Flow, Methods in computational physics, Wilkins,M. L; Alder, B; Fernbach, S; Rotenberg, M. pp211-263, Academic press, New York, 1964.

A4.8 LS-DYNA, Keyword user’s manual, version 960, Volume 1, 2001.

Page 298


Fig. A4.1: Oscillations in strain rate


Fig. A4.2: Typical Cr steel data sampled at SIMR


Fig. A4.3: Descriptive figure to show the use of initial frictional stress polynomial predictions


Fig. A4.4: Optimised flow stress curves with interpolated initial parts


Fig. A4.5: Extrapolation of the strain rate


Fig. A4.6: Optimised flow stress curves with temperature dependent initial frictional stress


Page 299

Fig. A4.7: Optimisation for Cr steel data at three different temperatures and four strain rates at each temperature


Fig. A4.8: Optimisation for Cr steel data a two different temperatures and four strain rates at each temperature


Fig. A4.9: Optimisation for C-Mn-Nb data at 850°C and at three strain rates


Fig. A4.10: Optimisation for C-Mn-Nb data at 900°C and at three strain rates


Fig. A4.11: Optimisation for C-Mn-Nb data at 1000°C and at three strain rates


Fig. A4.12: Optimisation for C-Mn-Nb data at 1100°C and at three strain rates

Page 300

Fig. A4.13: Optimisation for C-Mn-Nb data at 1200°C and at three strain rates.


Fig. A4.14: Optimisation of two different temperatures and at three strain rates mean RMSE 5.79 MPa


Fig. A4.15: FE model for simulating compression tests


Fig. A4.16: Compression of Cr steel at 900°C and 0.8/s


Fig. A4.17: Predictions at 900 and 3.5/s for Cr steel


Page 301

Fig. A4.18: Predictions at 900°C and 10/s for Cr steel


Fig. A4.19: Predictions at 900°C and 30/s for Cr steel


Fig. A4.20: Predictions at 900°C and at 10/s for C-Mn-Nb steel


Fig. A4.21: Predictions at 900°C and at 1/s for C-Mn-Nb steel


Fig. A4.22: Dislocation density at 900°C and at 10/s for C-Mn-Nb steel


Fig. A4.23: Effective stress at 900°C and at 10/s for C-Mn-Nb steel


Page 302

Fig. A4.24: Affective plastic strain at 900°C and at 10/s for C-Mn-Nb steel


Fig. A4.25: Vacancy concentration at 900°C and at 10/s for C-Mn-Nb steel


Fig. A4.26: Validation of simulations made with predicted initial frictional stress


Fig. A4.27: benchmark testing case, colour bar indicates the dislocation density level


Fig. A4.28: Effective plastic strain for C-Mn-Nb steel


Fig. A4.29: Effective stress for C-Mn-Nb steel


Page 303

Fig. A4.30: Vacancy concentration for C-Mn-Nb steel


Fig. A4.31: Temperature for C-Mn-Nb steel without thermal conduction through the rolls


Fig. A4.32: Effective plastic stress in mid and front positions


Fig. A4.33: Effective plastic strain in mid and front positions


Fig. A4.34: Effective strain rate in mid positions


VLEX uses login cookies to provide you with a better browsing experience. If you click on 'Accept' or continue browsing this site we consider that you accept our cookie policy. ACCEPT