Advanced constitutive modeling and application to industrial forming processes

Continuum constitutive descriptions of plasticity suitable for finite element simulations of sheet forming processes are succinctly discussed in this presentation. Although multi-scale approaches allow for a more explicit representation of the physical deformation mechanisms occurring at microscopic scales, they are usually not suitable for industrial applications because of the quick turnaround time needed for process design simulations. Therefore, advances in classical concepts such as plastic anisotropy and strain hardening are still very much in demand. Actual forming simulation examples conducted for the steel industry are presented for illustration purposes.


Introduction
Numerical simulations for industrial forming processes have been very useful in the last decades, since the start of the NUMIFORM conference series.They have strongly contributed to the decrease of time for process optimization.Progress in deformation mechanics, heat transfer and microstructure evolution and their time and space discretization have allowed a number of commercial codes to provide robust numerical tools for this purpose [1,2].There are still a number of challenging aspects regarding the code themselves not only on the mechanics of finite elements but also on boundary and interface conditions and on material constitutive descriptions.This articles only deals with the latter aspect but if the goal of numerical simulations is to achieve accurate predictions of processes above a certain scale, much effort is still needed in other areas as well.
The constitutive description must not only account for the physical aspects to describe the proper mechanisms but also for the input data needed and the identification of the material coefficients.
It must translate in mathematical terms the relationship between a number of variables with their evolutions and the stress-strain response within a framework set by mechanics and thermodynamic principles.A pure numerical description, which specifies the state of a material with a set of data is more empirical.It leads to poor extrapolation of states, which are not included in the database.However, sometimes, it is not practically possible to exclude this type of input as, for instance, for the case of hot press forming (HPF), such as schematically shown in Fig. 1.HPF simulations using thermo-mechanical-metallurgical numerical models require input information such as temperature and strain rate-dependent stress-strain relationship, which are not the only material characteristics needed.Thermal properties such as conductivity, latent heat, phase or temperaturetransformation-cooling diagram, and metallurgical states such as composition, grain size, phase or atomic diffusivity.It is also worth mentioning the coupling, sometime strong, between all these variables as well as with temperature and strain rate.Therefore, for such a problem, it is not wise to reject a priori any form of empirical input.
Another consideration necessary to choose a constitutive model is the scale of the material features needed, the size of the problem in terms of node or element number and the associated computation time.The scale is dictated by the microstructural information needed to characterize a product as, for instance, the grain size, which requires simulations at the meso-scale.It is possible to envision cases were details at a finer scale are necessary, for instance for the manufacturing of microdevices.However, in general, the description of lower scales conflicts with computation time.In sheet metal forming processes for automotive applications, the purpose of the simulations is to optimize the process geometry, boundary and interface conditions, load and, possibly, the material itself.Therefore, FE simulations must be run a number of times before an acceptable optimized solution is achieved.This prohibits the consideration of any material microstructural features explicitly.Thus, even if the idea of defining the constitutive description from the atomic scale is scientifically and philosophically very attractive, it is usually not practical.
The size of the simulation and the product itself also influence the nature of the constitutive description for the same reason, that is, computation time.An automotive panel of the size of a car door might require a large number of elements if the goal of the simulation is to detect plastic localization or fracture zones.In contrast, the stretch-forming of a panel for the fuselage of an airplane, several meters in length, may not be such a large size simulation because plastic deformations are usually small.Thus, the size of the FE element may be relatively large.Nevertheless, large product sizes imply generally that the number of elements or nodes is large, that is, the simulation of only one process condition long.
Other issues for practical industrial applications are the ability of a constitutive model to be well understood by metal forming analysts, to be identified at the lowest possible cost and to require a minimum number of experiments.Ideally, the best constitutive model for industry is a description that provides good accuracy with low computation time and necessitates only a few experiments for the coefficient identification.

Constitutive description example
In this example, some of the authors' own models are presented but there are many ways to describe a flow theory of plasticity.More emphasis is put on the reasons for the authors' assumptions since the detailed formulation can be found in the relevant publications [4][5][6][7].Room temperature deformation for rate-independent plasticity is considered in this section.Since the addition of a viscous effect can be done relatively simply, rate independent plasticity only is discussed.The goal is to define a yield condition, a flow rule and a hardening model.
The von Mises yield condition is the most popular for isotropy because it is one of the simplest.It contains some physics because it describes yielding as the limit of the elastic shear distortion energy, is in reasonable agreement with many experimental results and is independent of the hydrostatic stress as observed in metals.In this sense, it is as physical as the Schmid law, except that it is not expressed at the same scale.Similarly, Tresca expresses yielding based on a critical macroscopic shear stress, which, as a consequence, is independent of the hydrostatic stress.However, von Mises is more popular in numerical simulations because it leads to a smooth yield surface.Another isotropic yield condition is that of Hershey [8] where σ σ is the Cauchy stress tensor and ′ σ the principal values of the deviator.This yield condition is similar to that of von Mises but non-quadratic and is again based on an average of the principal shear stresses.Again, it contains as much physics as the Schmid law but no microstructure information.In addition, the yield surface shapes predicted with this criterion are in excellent agreement with those calculated with crystal plasticity model when the exponent a is set appropriately, usually significantly larger than 2 [9].The point of this paragraph is to say that von Mises, Tresca, Hershey and other suitable yield criteria are equally valid to represent isotropy and contain some physics but no microstructural information.
In order to introduce anisotropic formulations, a simple method is to use linear transformations of the type and to substitute the principal deviatoric stresses in Eq.
(1) by the principal values of the linearly transformed deviator [4,5].The advantages of this method are numerous.The yield condition reduces to the isotropic formulation if the 4 th order tensor C (or L ) reduces to the corresponding identity tensor.Moreover, the yield surface remains convex if the anisotropic yield surface associated with the isotropic yield function is convex.Usually, the convexity condition is more easily verified with principal stresses.C and L have a structure, which is consistent with the macroscopic symmetry of the material.Finally, when yield is defined conventionally with an offset strain of a fraction of a percent, it leads to yield surface shapes in agreement with experiments [10,11].
Under the isotropic hardening assumption, the yield surface only expands because it is controlled by one parameter, the reference stress σ R or, equivalently, the effective strain, ϕ .The former is defined by the plastic work equivalence, i.e., ξ dϕ = σ pq dε pq (3) Note that when ξ is von Mises, it is possible to retrieve the corresponding effective strain from Eq. ( 3).However, if the effective stress is defined by Eq. ( 1), Hershey's approach, or any isotropic or anisotropic yield function, the von Mises effective strain is not proper.The effective (or accumulated) strain always increases, unless specific recovery mechanisms are active, because plastic strains accumulate irrespective of the current strain increment.This corresponds to the physical observation that the density of dislocations, which produces slip before being trapped, increases.The reference curve σ R ϕ ( ) is mathematically expressed by an equation, which can take numerous forms and may or may not contain microstructural information.However, very often the stress is represented as a monotonic function of the effective strain because the material hardens as the accumulated strain increases.Thus, the effective stress and strain are physically motivated although they do not necessarily contain microstructural information.Experimental of iso-plastic work contours for DP590 steel approximated by Yld2000-2d (reprinted from [11] with permission from Springer).
Isotropic hardening is particularly convenient for industrial applications because it requires only a limited number of tests to identify the constitutive coefficients.In addition, the use of shell elements for sheet forming problem reduces the computation time drastically.As a result, a state of plane stress is usually assumed.Thus, a form of Eq. ( 1) modified with two linear transformations to describe plastic anisotropy, the so-called plane stress Yld2000-2d yield function, is more convenient.It contains eight anisotropy coefficients, which require a minimum of four mechanical tests to identify.Three uniaxial tension tests in different directions provide six data since each test allows the determination of a flow stress and r-value (width-to-thickness strain ratio in uniaxial tension).Since the bulge test has recently become a standard, it can be used in addition to provide input data for balanced biaxial tension.However, other tests could be used as well.In addition, new techniques, such as the virtual fields method, are now emerging to identify a full set of constitutive coefficients [12].These methods are based on the measurement of full inhomogeneous displacement fields with digital image correlation.Researches are conducted to extract all the anisotropy coefficients from one single test, which would reduce considerably the number of experiments to conduct.
However, this approach, although well established for elasticity, requires much more research for plasticity.
Finally, the question of the flow rule has been debated for several decades and a resurgence of the discussion has occurred recently [13,14].The associated flow rule (AFR) assumes that the strain increment is normal to the yield surface while the non-associated flow rule (NAFR) assumes it is not.The constraint of the AFR is so high that it does not appear to be reasonable.However, general considerations of crystal plasticity with only the Schmid law for the crystal behavior as the main assumption indicates that metals, which deforms by microscopic shearing, should deform according to the AFR.Thus, even if the AFR is not strictly respected, the physics of plastic deformation indicates that the strain increment should be reasonably close to the normal of the yield surface.Because it is difficult to precisely define AFR and NAFR, the low scale physics for metals tends to justify the use of the AFR.This is why in this work, the AFR is considered more suitable, that is The isotropic hardening assumption leads to reasonable simulation results when loading is roughly proportional.However, when sharp changes in strain path occur, as often the case in forming, the accuracy of the predictions is usually not as good.This is particularly true when load reversals occur such as bending and unbending around a die radius, leading to the so-called Baushinger effect.The springback, that is, the elastic shape distortion of a part due to the removal of forming loads, strongly depends on the Bauschinger effect and the stress-strain behavior after reversal.Since springback is a critical problem for high strength sheet materials such as modern steels, constitutive models that account for the Bauschinger and associated effects are important for industrial applications.This effect can be captured by non-linear kinematic hardening rules [15,16], in which the yield surface translates in stress space.Kinematic hardening has been used very efficiently during a half of a century to account for anisotropic hardening.However, a different approach was recently proposed in an effort to make industrial applications easier.More than the model itself, the arguments to justify the use of this model are developed below.The so-called homogeneous anisotropic hardening (HAH) approach makes use of the isotropic hardening concept discussed above.
The yield condition is described by q ĥ : ′ σ − ĥ : ′ σ q { + f 2 q ĥ : ′ σ + ĥ : ′ σ where f 1 , f 2 are two scalar state variables that evolve rapidly as a function of strain and distort the yield surface.When f 1 = f 2 = 0 , the standard isotropic yield condition evolving under the isotropic hardening assumption as seen above is recovered, that is, ξ ′ σ σ ( )= σ R (ε ) .ĥ , a tensorial state variable called microstructure deviator, controls the location of the distortion on the yield surface.This variable evolves slowly as a function of the strain because it is assumed to represent the dislocation structure that develops in a material.Although this model is macroscopic, it is consistent with a number of observations made about dislocation patterning under monotonic and nonmonotonic loading conditions, which indicate that dislocation structures tend to evolve slowly towards the structure corresponding to a particular loading direction.Thus, the material is physical, although it does not explicitly include microstructural features.The material behavior predicted with this approach when the strain path is linear is strictly the same as that obtained using the anisotropic yield function considered under isotropic hardening.Although, there is no back-stress or kinematic hardening component, this approach successfully describes the Bauschinger effect as illustrated in Fig. 3 for a balanced biaxial stress state followed by load reversal.
Figure 3. Yield surface evolution for DP490 steel calculated for as received material, after balanced biaxial tension and after subsequent balanced biaxial compression (reprinted from [17] with permission from Elsevier).

Applications
A few examples are given in this section that underline the progress and current limits of macroscopic material models.

Indentation
The first example is an application to indentation [17].
The process is schematically explained in Fig. 4a.The process consists of two steps: 1) A flat sheet is first prestretched in balanced biaxial tension with a large punch to a given strain and subsequently deformed by an indenter operating in a quasi-static mode.Preliminary simulations indicated that elastic deformation is preponderant in the whole specimen since the maximum plastic deformation near the indented area is less than 10%, loading is highly non-linear with biaxial reversals and the elasto-plastic behavior of the material near balanced biaxial tension requires special attention.Thus, in order to accurately describe the behavior of the material, the biaxial elastic modulus as well as its evolution as a function of the accumulated strain has been accounted for.Plastic anisotropy, strain hardening and its transient stage during reversal require a good description based on biaxial measurements, as well as simple shear testing with load reversal.The former allows the characterization of the biaxial elasto-plasticity in proportional loading and the latter, the transient hardening effects.Additional uniaxial tension tests are necessary to characterize plastic anisotropy.In order to compare the effect of the constitutive description, two models were employed in the simulations.The first, referred to as conventional, consisted of a constant uniaxial elastic modulus, von Mises yield condition with the associated flow rule and Swift isotropic hardening description.
The other, referred to as advanced constitutive model, consisted of the biaxial elastic modulus and its decrease as a function of the accumulated plastic strain, the Yld2000-2d non quadratic anisotropic yield function, and the HAH distortional framework to account for the Bauschinger effect and other transient hardening phenomena.The simulations were conducted with Abaqus/Standard (implicit).Figure 4b compares the experimental indentation load for a DDQ steel as a function of the indenter displacement with the predicted results computed with the conventional and advanced models.This figure clearly indicates the superiority of the advanced model.Of course, the better material description leads to the better results.This implies that a careful material characterization is necessary and requires tests, which are not necessarily standard.This implies an additional significant effort to gain accuracy.For practical purpose, the permanent indent depth is about 30% higher with the advanced model compared with the standard.However, the computation time is significantly higher with the advanced model.If this is not an issue for this example, this could be prohibitively different for much larger problems.As a whole, accuracy comes with a cost that industry may or may not be willing to assume.The important question is what maximum error/budget ratio is allowed for a particular process and product.

Fuel cell separator
The next example is about forming of ultra-thin sheets (less than 0.15 mm thick) for PEM fuel cells (PEMFC) [18].Fig. 5 shows the geometry of the mesh of the cathode, which has a size of 400 by 220 mm.The blank, initially flat, becomes corrugated with micro-channels after forming.A channel is less than 1 mm wide, approximately 0.5 mm deep and features corner radii of 0.3 mm.The forming consists of three steps, namely, forming, trimming and springback.Experimentally, the PEMFC cathode was formed on a HSK direct-drive digital servo-press, which, because of its high accuracy compared to traditional presses, is necessary to successfully make product from ferritic stainless steels.
This problem is very challenging from many viewpoints.The material itself is so thin that it is very difficult to measure the mechanical properties.In fact, the standards for sheet metal characterization are not valid for thickness values lower than 0.5 mm.Tests for extra-thin sheet materials require extreme care and are sometimes impossible to conduct as, for instance, reversals for the characterization of the Bauschinger effect.In addition, since the grain size is not negligible compared to the thickness, this feature might be needed in the characterization input data, as well as a model to transfer this information to the continuum scale.In passing, another issue with many thin materials is that they are often heavily deformed and do not provide sufficient elongation in uniaxial tension for characterization.Thus, in spite of all the advances made with material property measurements, it is expected that additional efforts in this area will be necessary in the future.
For the simulations, even though the size of the product is modest, this problem is already very large.The simulation were run with about 1.5 millions shell elements although many more are required because the mesh size of 0.25 mm is barely smaller than the radius of the channel corners.In spite of that, the computation is prohibitively too large for process optimization in industry.As an aside, sheet forming process optimization is almost always conducted with shell elements and explicit solvers in industry.One of the reasons is computation time while the other is numerical convergence.Due to the many tool/workpiece contacts that occur in forming, it is very difficult to solve a complete problem, even a small one, using an implicit solver.The use of solid elements is generally not considered because of computation time.
In this example, in order to evaluate buckling and twist springback, and to optimize drawbead positions, a simple constitutive description only, similar to the conventional model of the previous section (von Mises, etc.), was considered in order to keep the effort at a reasonable level.Nevertheless, a realistic solution of the drawbead positions was obtained.The plastic flow localization, which occurs during forming (Fig. 5b), was analyzed using a plane strain analysis throughout the channel section.The search of an optimum process condition was reached for this particular example but real products are expected to be commercialized only after down-gauging.
Therefore, a number of serious challenges still remain in numerical and experimental methods in industrial forming processes for such products.

Hot press forming
The hot press technology (HPF) is very important for the automotive industry because it allows the manufacturing of products such a B-pillar from boron-added steel, which leads to very high strengths because of the martensitic structure.However, the material must be formed at high temperature in the austenite domain and die-quenched at a high cooling rate to avoid diffusive and softer phases to transform.Also of importance is to tailor the properties of the product by producing zones with full martensitic microstructure and others interspersed with softer phases.The use of tailored-welded blank is another option [19,20].For instance, for the B-pillar example, the upper end should be very strong for, during a crash, resistance to intrusion in the passenger cabin and the bottom part more ductile for energy absorption.The partial cooling technology is among a number of technical solutions to this end.
In this process, forming takes place approximately from about 800 ºC down to 600 ºC.Then, the critical die cooling stage occurs from the temperature at the end of forming stage down to about 200 ºC.Segments of the die are cooled using water circulating in channels at various flow rates while other are not (Fig. 6a).This produces cooling at different rates for different locations in the blank.This process is difficult to model because it is time and temperature dependent.This requires a fully coupled thermo-mechanical-metallurgical solution [21,23].The plastic behavior itself, stress-strain and localization, must be temperature and strain rate dependent.Then, phase transformation kinetics must be known accurately, which demand the knowledge of a heat transfer description.All these parameters are coupled since heat transfer coefficients, for instance, are at least pressure and temperature dependent.In addition, the material itself is subjected to drastic microstructural changes associated with the transformation kinetics.Of course, all this requires a very large quantity of input data, some of it very empirical, based on composition, or extracted from thermodynamics databases and software.
An example of such simulation is given in Fig. 6a for a non-straight channel with a U-shape cross-section.The die was segmented in order to produce different cooling rates.The cold die zone leads to full martensite in the blank, while the heated zone can produce softer phases, depending on the die segment temperature.Fig. 6b compares the cooling curves from the start of forming to a time after which, the microstructure remains stable.The predicted cooling curves are in good agreement with the experiments but the preliminary results were not.As a result, it was necessary to measure the heat transfer coefficients between a particular tool steel and the blank as a function of the temperature and the contact pressure.
The Vickers hardness can be calculated, yet from an empirical formula based on phase volume fraction and composition.It leads to predicted values that are in good agreement with experiments (Fig. 6c) but at the expense of additional costs and delays.This example indicates that empirical methods are necessary to completely solve such a problem.A more scientific approach to this problem with FE and more physically-based models, namely, constitutive, heat transfer, metallurgical states and reaction kinetics, is not readily available yet.This example suggests that there are still huge efforts to conduct just from the material behavior viewpoint.In spite of this, the results of thermo-mechanicalmetallurgical finite element simulations can be a tremendous help for process designers and practitioners.

Discussion
Advances in continuum constitutive modeling have improved the accuracy of numerical simulation methods in industrial forming processes.In the sheet metal forming area, a number of obstacles still remain in order to produce reasonably accurate results in a time and cost efficient manner, even with continuum constitutive models.The requirement of high turnaround time and low cost requires the use of compromises that prevent the achievement of higher accuracy.
Lower-scales approaches (Fig. 7) are not suitable yet for direct applications in sheet metal forming codes because huge amounts of microscopic information should be preprocessed first, one way or another.Nevertheless, mesoand micro-scopic material models can be used to guide and refine the development of higher scale constitutive descriptions (Fig. 8).
In addition, new experimental techniques and procedures are still required in order to provide sufficiently accurate input data to the numerical forming process codes.Again, limitations arise because of the complexity of the database to prepare, such as for HPF, or because of the nature of the product itself, such as extra thin sheet materials, which are hurdles to overcome.At the same time, all these obstacles constitute new challenges for science and engineering.Alternate routes for the description of the macroscopic behavior are showing up such as full field measurements and virtual field methods [12], which provide a wealth of information about the material behavior.However, this information has to be processed first before it can be used optimally for forming application.Finally, other forms of input data can be generated from virtual experiments conducted with lower scale microstructure-based models.In any case, even though the field of simulation for metal forming and advanced constitutive modeling, particularly in industrial forming processes, is not an emerging technology, this field is still vibrant and full of innovations yet to discover.This is also "fuelled" by the imperative necessity of sustainable development, which drives the emergence of lightweight materials with enhanced attributes.This exciting time is coming with challenges for the development of advanced constitutive descriptions as well.

Figure 4 .
Figure 4. Indentation of panel for drawing quality steel (reprinted from [17] with permission from Elsevier); a) Testing schematics and; b) Load-displacement curve.

Figure 7 .
Figure 7. Multi-scale description showing the huge amount of input data to encapsulate in simple concepts for application to manufacturing.

Figure 8 .
Figure 8. Application of multi-scale approaches through transfer towards continuum constitutive plasticity models.