Numerical Evaluation of Surface Roughness Influences on Cold Formability of Dual‐Phase Steel

DP1000 steel sheets with different surface treatments are taken to three‐point bending tests to evaluate their damage and fracture behavior. It is found that different surface treatments yield different fracture behaviors. Conventional finite‐element (FE) simulation with an uncoupled material model is able to capture the material's response only for the sheets with the perfectly smooth surface condition. A multiscale strategy is introduced to quantitatively include the surface information into the material model. The localization at the microlevel simulation allows adjusting the fracture criterion of the surface elements in the macrolevel simulation. Therefore, a good agreement between FE simulation and experiment is yielded.


Introduction
The automotive industry has been in a highly competitive situation to find the best material to fulfill its requirements in terms of light weight, strength, crashworthiness. Advanced highstrength steel (AHSS) has been developed and widely used to satisfy the car producers for these purposes. Dual-phase (DP) steel is one example of this steel type. It is one of the most popular choices for multiple reasons. Its unique microstructure, i.e., ferritic matrix with martensite islands, plays a crucial role to enhance the performance of this steel class. In general, the plastic behavior as well as fracture behavior of DP steels depends not only on the ferrite and martensite phase fraction but also on grain refinement and carbon content.
Cold formability is another key factor that manufacturers consider. The forming limit diagram (FLD) has been one of the most classic tools to evaluate the steel's inhomogeneity. [1] This concerns the "global" formability of the steel. In other words, it includes the deformation of the material in large regions simultaneously. The local formability, on the other hand, reveals ability of the material to withstand the concentrated deformation. [2] Examples of standard for the local formability test include the hole expansion test [3] and three-point bending test. [4] The DP1000 steel yields a very good global formability but lacks its local formability. A small amount of extrinsic influences such as residual stress from sheet manufacturing or surface inhomogeneity could lead to local failure.
Dual-phase (DP) steels generally fail in a ductile manner when subjected to mechanical loads at room temperature. The involved mechanisms include void nucleation, void growth, and void coalescence. These mechanisms accumulate as strain develops. There are two main approaches to model such mechanisms and predict the material's failure in the macroscopic scale-uncoupled and coupled models. The uncoupled models generally define a ductile failure criterion based on plastic strain value where the failure is observed. The criterion designates a boundary between failed and not-failed material points. Nevertheless, they take no consideration of the fact that the internal mechanisms degrade the not-failed material point's response before it fails. The Johnson-Cook failure model is one of the most famous models of this type. [5] The coupled models, on the other hand, consider the void relevant mechanisms during the material's plasticity by taking an effective void volume term as well as void evolution term into account. They soften the material's response from yielding until fracture. A renowned model of this type is the Gurson-Tvergaard-Needleman (GTN). [6,7] Both modeling approaches have their own pros and cons. The uncoupled model can recognize the material's fracture instant while acknowledging no plasticity degradation. In contrast, the coupled model realizes the internal damage while disregarding the fracture criterion. To combine these ideas, a hybrid approach called the "modified Bai-Wierzbicki (MBW) model was introduced. [8][9][10] This model allows us to include the damage softening effect before fracture as well as define a strain limit as a function of stress state to identify its fracture criterion. Therefore, after yielding, instead of modeling the material in two stages-failed and not-failed-one can classify it into three stages: plasticity, damage softening, and failed. For material such as DP1000, whose damage softening stage is relatively short, one can apply MBW and a specific parameter set to consider it as an uncoupled model. This study presents experimental results of DP steels with two extreme surface treatments under a cold forming process. It is discovered that not only internal phenomena such as void relevant mechanisms contribute to the DP steel's failure behavior but also external factors such as surface condition. [11][12][13][14] DOI: 10.1002/srin.202000141 DP1000 steel sheets with different surface treatments are taken to three-point bending tests to evaluate their damage and fracture behavior. It is found that different surface treatments yield different fracture behaviors. Conventional finiteelement (FE) simulation with an uncoupled material model is able to capture the material's response only for the sheets with the perfectly smooth surface condition. A multiscale strategy is introduced to quantitatively include the surface information into the material model. The localization at the microlevel simulation allows adjusting the fracture criterion of the surface elements in the macrolevel simulation. Therefore, a good agreement between FE simulation and experiment is yielded. Until now, there has not been a model that quantitatively takes the surface condition into account as a failure criterion or damage-induced softening function. Modification factor approaches have been developed to fulfill the missing information from surface imperfection. [15,16] Because the aforementioned macroscopic material models are generally implemented via finite-element method (FEM) software in which scaling limitation exists, creating a simulation model in millimeter scale while preserving micrometer-scale geometry from the surface condition would take the simulation to last essentially forever. This study introduces a multiscale strategy to extend the prediction capability of the MBW model.

Material
A dual-phase steel DP1000 was chosen for this study. It is a widely used steel grade in the automotive industry. The chemical composition is shown in Table 1. For detailed microstructure information of the steel, readers are referred to a previous study of the investigating team. [17] The sample was prepared by grinding, polishing, and etching with HNO 3 . A previous in-depth comparative study between DP1000 and DP600 can be found in a previous study. [18] It reveals the differences between DP steel grades As DP1000 contains higher martensitic phase fraction-to be specific, 62% ferrite and 38% martensite-the strength of DP1000 is much higher than that of DP600.
As the rolled material is expected to achieve anisotropic property, this study considers only the rolling direction because only a slight difference was found from our preliminary study. A series of quasi-static uniaxial tensile tests in its rolling direction at room temperature were performed to investigate its basic mechanical properties as well as its flow curve. The steel's yield strength of 770 MPa, ultimate tensile strength of 983 MPa, uniform elongation of 5.2%, and A 80 fracture elongation of 11.0% were discovered. A Swift-Voce hardening law was selected, as shown in the following equation [19] where σ is extrapolated equivalent stress, ε 0 is yield strain, and ε p is the equivalent plastic strain (PEEQ). According to the previous study, the damage initiation strain and the ductile fracture strain of DP1000 are close to each other. [18] It confirms the fact that applying an uncoupled model allows us to sufficiently capture the material damage and fracture behaviors.

MBW Material Model with Surface Factor
The MBW model has proved its capability to predict the ductile damage behaviors of various types of steels. The model defines damage and fracture on a macroscopic scale by local strain and local stress states. In other words, it considers the local PEEQ ε p , which is triaxiality η and lode angle parameter θ dependent, as the main variable to drive damage and fracture phenomena. Given principal stresses σ 1 , σ 2 , and σ 3 , the three invariants p, q, and r of the stress tensor read The stress triaxiality η and lode angle θ are defined by As the lode angle is in symmetry, it shall be characterized in the range of 0 ≤ θ ≤ π=3. Therefore, the normalized lode angle or lode angle parameter is calculated by While loading is applied but the damage criterion is not fulfilled, the conventional yield potential is used to describe the plastic deformation of the material. Once the damage condition is reached, the damage softening effect governed by the damage variable D is taken into account onto the yield stress σ y . The yield potential Φ is expressed by To indicate the initiation of ductile damage, a so-called damage initiation locus (DIL) is defined as a relation between PEEQ and stress states, which is formulated by where c i 1 , c i 2 , c i 3 , and c i 4 are DIL parameters to be calibrated from macroscopic experiments. The surface factor c s is an additional parameter to represent how surface imperfection softens the local material. This factor is in the range of 0 ≤ c s ≤ 1 for material points on the surface, and c s ¼ 1 otherwise. To calibrate it, a multiscale analysis and surface modeling shall be conducted. It is determined by taking the ratio between the damage initiation plastic strain at the meso level ε i meso and the damage initiation of plastic strain at the microlevel ε i micro . More details on surface factor calibration are given in the upcoming paragraphs. Upon damage initiation, the damage variable D evolves as plastic strain develops. The relation between damage D and the PEEQ ε p is assumed to be linear, which reads where σ y;i is yield stress at damage initiation, G f determines energy dissipation between damage initiation and fracture, and ε f defines fracture strain, which is a function of stress states, in which its parameter c f 1 À c f 4 are calibrated from macroscopic experiments. This equation provides hybrid behaviors between the coupled and uncoupled model. Nevertheless, it is observed that DP1000 yields a very short damage accumulation phase during the macroscopic experiments. In other words, the damage initiation strain is close to the fracture strain ε i % ε f . It can be described by setting G f to a small value as well as setting c f 1 À c f 4 to be slightly higher than the corresponding c i 1 À c i 4 . By applying this assumption, the MBW model behaves in an uncoupled way. From the experimental results, the calibrated MBW parameters of c i

3-Point Bending Test
The bending test was performed according to Verband der Automobilindustrie (VDA-German Association of the Automotive Industry) 238-100. [4] DP1000 sheet of size 20 Â 60 mm with 1.0 mm thickness was placed on two 18 mm diameter support rollers. The prepared surface, whether grinding or polishing, faced downward. The support rollers had a distance between surface and surface of 6 mm. A punch with a radius of 0.5 mm pushed the sheet from the top surface with 12 mm travel distance at the speed of 1 mm min À1 and at room temperature. The crack on the steel sheet was expected to appear on the bottom surface due to plane strain tension loading. The test schematic is shown in Figure 1.

Surface Preparation
This study aims to investigate and model the surface condition of the steel during forming as the main factor for its fracture behavior. Therefore, two surface conditions were prepared from an identical steel sheet. The sheet was cut into samples of 60 Â 20 mm where the 60 mm length was in the rolling direction. Both surface conditions were ground by 80-grit sandpaper. One set, called "smooth", was fine-grinded and polished with the grain size of 1 μm afterward. Another set, namely, "rough", was kept as is after the first grinding process. Considering the thickness of the sheets, the smooth set was processed by grinding and polishing, while the rough set was processed by only grinding. Therefore, the thickness of the smooth set is supposed to be slightly thinner than that of the rough set. In terms of visual appearance, the surface conditions of smooth samples and rough samples are shown in Figure 2a,b respectively. The difference between the surface of both sample sets can be noted by naked eyes.
To quantify the difference, they were taken to a white light confocal microscope analysis to investigate their surface topography. The magnification was 20Â, which covered a 1.475 Â 1.475 mm area. Each sheet was measured at the positions (x, y) as follows-(30, 5), (30, 10), (30, 15), where the direction of x, y and the reference point (0, 0) are indicated in Figure 2a,b. The measurement result of "smooth" and "rough" surface samples are shown in Figure 2c,d respectively. Please note that small jitters in Figure 2c are due to light reflection during the measurement process, not the real surface's height. They were eliminated in the noise reduction process after the measurement.
The measured data were decomposed into waviness and roughness components by using a Gaussian filter. [21] The transfer function from the raw surface profile a 0 and roughness surface profile a r is defined by where λ is the wavelength of raw surface profile and λ c is the predefined cut-off wavelength. The cut-off wavelength λ c plays a role to identify what is the biggest wavelength component that shall be included in the analysis. It should be set according to the mesoscale simulation base length, which corresponds to the finest mesh size of the component level simulation FE model. In this study, the cut-off wavelength λ c was 50 μm. The filtered height distributions of both surfaces for each sample set were accumulated and assumed to be normally distributed. Therefore, the roughness component was fitted to the normal distribution equation as follows f ðxjμ, σÞ ¼ 1 ffiffiffiffiffiffiffiffiffi ffi 2σ 2 π p exp ðx À μÞ 2 2σ 2 (14) where μ is the mean value of the roughness and σ is the standard deviation of the roughness. The fitting takes only the roughness component into account and finds these two parameters.
The fitted values are μ ¼ À0.1244 μm and σ ¼ 0.0429 μm for the smooth set and μ ¼ À1.2970 μm and σ ¼ 0.6293 μm for the rough set. It can be noted that the roughness deviation of the smooth set is close to zero. Therefore, this surface set can be considered as a perfectly smooth surface representation.

Simulation Setup
A multiscale simulation is applied in this study to transfer the surface roughness information to component simulation. First of all, the three-point bending test is considered as a component level simulation. The first trial assumes that the surface of the steel sheet is perfectly smooth. It allows finding the location, time instant, and stress state of the critical region as the bending proceeds. This information is transferred to the mesolevel simulation, which takes only a 50 Â 50 μm area. The measured surface roughness is considered at this level. The simulation result of the local region in the meso level is taken into account as microlevel. As a result, the surface factor c s is found by bridging the information between component level and micro level by the meso level.

Component Level Model
The overall geometries such as sheet dimensions, support diameter, support distance, and punch geometry are set according to the bending test mentioned in the previous paragraph. The isometric view of a finite-element (FE) model created in ABAQUS is shown in Figure 3. The support rollers and punch are regarded as rigid bodies. The only deformable body in this simulation is the steel sheet. It is meshed by C3D8R elements. The supports are fixed concerning their displacement and rotation. The punch moves in the vertical direction downward by 12 mm travel distance. The global mesh size is 0.30 mm. The regions that contact the support rollers are meshed by 0.15 mm elements. The region that contacts the punch and the opposite side is meshed by a 0.05 mm element. Figure 3 shows how it is meshed on the right half of the sheet. The friction coefficient between the sheet and support and that between the sheet and punch are set to 0.2. It can be noted that the FE model in this level cannot include roughness into its geometry as the finest mesh here is 0.15 mm or 150 μm, while the roughness has its length unit %4 μm. The sheet material's constitutive model is defined by the MBW material model implemented in the VUMAT user subroutine in Fortran. The parameters are set according to calibrated DP1000 values in Figure 3. However, the c s is set to 1.0 to all elements in the first simulation to obtain the stress states, PEEQ, and fracture behavior of the sheet at its perfectly smooth surface condition. After obtaining information from the other scale's simulation, the c s of all elements on the bottom surface will be altered to represent the effect of surface imperfection.

Submodel Generation
The simulation in the meso level acts as a bridge between the component level and the micro level. As the mesh at the critical point on the sheet is of size 50 μm, the mesolevel model is defined accordingly by a rectangle with a base size of 50 Â 50 μm. On the top side of the rectangle, points are randomly placed according to normal distribution with the mean value and standard deviation value from Figure 4. These points are connected to each other by polynomial lines to construct geometrical roughness features. The representative submodels with smooth sample surface as well as rough sample surface information are shown in Figure 4a,b respectively. The submodel is meshed, in general, with 1 μm elements, while the top 20% of the height is meshed with 0.6 μm elements.
It is noted that the implementation of the MBW user subroutine is not applicable to the 2D element. Therefore, a small thickness is included in the geometry. As the loading case in the component level is three-point bending, it is expected that the stress state of the critical element is plane-strain tension. Therefore, nodal displacement in the thickness direction is constrained such that the plane strain condition is satisfied. Boundary conditions on the left side as well as the bottom side of the submodel are set as symmetrically in xand y-direction, respectively. The displacement boundary condition in the x-direction is applied on the right edge of the submodel. In this particular case, the displacement is 20% of the submodel size, which is 10 μm.
As the current application of interest is VDA three-point bending, the critical element loading is plane-strain tension loading. This allows the 2D simulation scheme to capture the material behavior accurately. However, in case the critical elements experience a complex loading or a nonproportional loading, it is better to create a sub\model in a 3D fashion. Further development on creating a 3D submodeling scheme that integrates the surface information is ongoing.

Multiscale Simulations
The FE simulations in the component level model are conducted using the MBW material model. The resulting reaction force and punch displacement are extracted from the simulation output database. In addition, the local responses, i.e., stress triaxiality η, lode angle parameter θ, and PEEQ, of the critical element are stored. Please note that the stress state over PEEQ evolution is called the strain path. This information indicates that the stress state of the critical element during the three-point bending process, which is η % 0.56 and θ % 0, is plane-strain tension. Therefore, the boundary conditions of the mesolevel simulation are defined as plane-strain tension, as mentioned in the previous paragraph.
The submodel for both smooth and rough surfaces are simulated. The volume-averaged strain path of all elements in the submodel is referred to as the mesolevel response. The responses from both surface types yield a good agreement with the response of the critical element in the component-level simulation. This indicates a scale-transferability of the proposed strategy. Therefore, from this point on, the mesolevel response in the submodel can represent the critical element's response in the component-level simulation. Figure 5 shows localization in the submodels. The volume-average response of elements within a 5 Â 5 μm rectangle around the localized region is defined as the microlevel response. The 5 μm region is selected as it is the physical detectable microscopic damage. [9,22] The strain path plot comparison and the procedure to calibrate the surface factor c s of the rough surface model are shown in Figure 6. First of all, the strain path of the meso level, in other words, the strain path of the critical element of the component level, and the DIL derived from Equation (9) with c s ¼ 1.0 are plotted in Figure 6a. Without considering the submodel, the component-level simulation predicts the damage initiation at the time instant t i . Once the submodel is taken into account, the strain path of the microlevel response is plotted in Figure 6b. It is found that, at the time instant t i , the microlevel strain path exceeds the DIL significantly. In fact, it reaches the DIL at the time instant t k where t k < t i . The corresponding time instant t k of the mesolevel response is identified in Figure 6c. The surface factor c s is calculated by the ratio between the PEEQ of the mesolevel response at the time instant t k and the PEEQ of the microlevel response at the time instant t k . Therefore, the c s for the rough surface sample yields c s ¼ 0.344. As this surface factor shifts the DIL of the material, the modified DIL is shown in Figure 6d. For the smooth-surface model, the strain paths of the meso level and micro level do not deviate from each other significantly. Only slight localization is found. Thus, the surface factor c s ¼ 1.0 is applied.
The modified DIL is applied onto the constitutive material model of the top surface elements in the component-level simulation model to represent its surface condition. The material model of the bulk material below the surface remains original.

Validations
The DP1000 sheets of the "smooth" and "rough" set are bent under the VDA three-point bending test. After bending, it can be noted that no crack is observed from the smooth sheet, while major cracks are found by naked eyes on the rough sheet. Figure 7a,b shows the bottom surface after bending of the smooth sample and rough sample, respectively. It can be concluded from the experiment point of view that surface condition plays a role in the material's mechanical behaviors.
Not only can qualitative difference be observed from these two surface conditions, but quantitative difference can also be observed in Figure 8. The force-displacement responses of the rough surface sheet show an earlier force drop during the bending. This is due to ductile damage softening and crack propagation. On the other hand, these phenomena are not presented during the loading.
By applying FEM simulation at the component level using the MBW constitutive model without considering the surface  condition, in other words, setting the surface factor c s ¼ 1.0, the force-displacement of the simulation agree with the smooth surface sample's result nicely. Once the surface condition is considered, the calibrated surface factor for the rough sample c s ¼ 0.344 is applied on the surface elements. The force-displacement response of the simulation demonstrates an early drop in the same manner as the rough surface experimental result.

Conclusion
1) This study identifies influences of the surface condition in damage and fracture behavior for DP1000 during a VDA three-point bending process by comparing experimental results from sheets that are prepared by different processes. The rough ones were ground while the smooth ones were ground and polished. 2) The rough ones pronounced earlier damage and fracture mechanisms than the smooth ones during the three-point bending test. This indicates a significant contribution of the surface condition in the damage and fracture behavior of the material.
3) A hybrid constitutive material model, MBW, with an additional surface factor c s is introduced together with its calibration scheme by a multiscale simulation strategy. 4) The finite-element (FE) simulation in the component level with the MBW material model can capture the material's mechanical behavior of the smooth surface sheet precisely. 5) A submodel is generated by including the surface information from white light confocal microscopy. This submodel represents the critical element in the component-level model. 6) Localization from the surface geometry is captured and taken into account for surface factor c s calibration. This factor modifies the damage initiation criteria of the material point. 7) The calibrated surface factor is applied to every element on the sheet's surface to reflect its softened response due to surface condition. 8) With the proposed strategy, the componentlevel simulation can capture the early damage and fracture of the rough surface sheet accurately. 9) The current approach is limited to applications in which the critical element loading condition can be considered as 2D loading with a proportional strain path. In case the critical element requires a complex loading condition, a 3D submodeling approach should be applied instead.