A multi-barrier model assisted CAFE method for predicting ductile-to-brittle transition with application to a low-carbon ultrahigh-strength steel

The conventional micromechanical approaches today are still not able to properly predict the ductile-to-brittle transition (DBT) of steels because of their inability to consider the co-operating ductile fracture and cleavage mechanisms in the transition region, and simultaneously to incorporate the inherent complexity of microstructures. In this study, a complete methodology with coupled cellular automata finite element method (CAFE) and multi-barrier microcrack propagation models is presented to advance the prediction of DBT. The methodology contains three key elements: (i) a multiscale CAFE modelling approach to realize the competition between ductile damage and cleavage fracture and embrace the probabilistic nature of microstructures, (ii) a continuum approach to estimate the effective surface energy for a microcrack to penetrate over particle/matrix interface, and (iii) a method to calculate the effective surface energy for the microcrack to propagate across grain boundaries. The prediction of DBT therefore needs only (1) the stress-strain curves tested at different temperatures, (2) the activation energy for DBT, (3) the ratio between the size of cleavage facets and cleavage-initiating defects, and (4) key statistical distributions of the given microstructures. The proposed methodology can accurately reproduce the experimental DBT curve of a low-carbon ultrahigh-strength steel.


Introduction
Comprehending the factors affecting the ductile-to-brittle transition (DBT) of structural materials, e.g. steels, is of paramount importance to structural integrity assessment. The prediction of the DBT of ferritic steels has been widely carried out (Needleman and Tvergaard, 2000;Rossoll et al., 2002;Samal et al., 2011;Samal et al., 2008;Tanguy et al., 2005a, b;Tvergaard and Needleman, 1986;Needleman, 1988, 1993) by using approaches based on Gurson model (Gurson, 1977;Tvergaard and Needleman, 1984;Zhang et al., 2000) or Rousselier model (Rousselier, 1987) coupled with the RKR model (Ritchie et al., 1973), or the local approach, e.g. Beremin model (Beremin et al., 1983). However, these predictive approaches have their limitations. Firstly, these approaches represent in general a type of post-processing solutions, which evaluate cleavage fracture by the stress ahead of a crack tip calculated from a ductile constitutive model. The competition and interaction between ductile damage and cleavage in the transition region are absent. Secondly, it is difficult for finite element method to process two types of fractures within a single mesh since ductile fracture and cleavage occur within two independent material length scales, i.e. the ligament between dominant voids and the cleavage fracture unit, respectively. Thirdly, fracture toughness is strongly dependent on the given microstructures, which are often inadequately described (Shibanuma et al., 2018) due to their true complex probabilistic nature. Although the distribution of defects, e.g. particles in a material can be considered in the local approaches, other intrinsic statistical features, i. e. the distribution of initial voids, grain size, and misorientation of grain boundaries have also to be involved, so that the scatter of fracture toughness can be accurately captured. Only by including all these factors, the correlation between microstructure and toughness can be revealed.
One multiscale approach, i.e. the CAFE method (Das et al., 2006;Shterenlikht, 2003;Shterenlikht andHoward, 2004, 2006;Shterenlikht and Margetts, 2015;Wu et al., 2005), which couples cellular automata (CA) and finite element (FE), offers a viable alternative framework for the prediction of DBT, and overcomes the above-mentioned shortcomings. In order to reproduce DBT of a ferritic TMCR steel, a temperature-dependent threshold of grain boundary misorientation θ th has been utilized previously Wu et al., 2005). However, the proposed temperature-dependent θ th is heuristic and lacks of physical significance. Therefore, the key for success of the CAFE method lies in the proper identification of physically-based variables to capture the dependence of fracture toughness on temperature in the transition regime (Li et al., 2019).
A material constant σ u in the Beremin model correlates with the effective surface energy γ eff for brittle fracture (Beremin et al., 1983;Lambert-Perlade et al., 2004). Recently, it has been found that a temperature-dependent σ u can reflect the effect of temperature on fracture toughness (Cao et al., 2011;Gao et al., 2005;Moattari et al., 2016;Petti and Dodds, 2005;Qian et al., 2015;Tanguy et al., 2005b;Wasiluk et al., 2006). The temperature-dependence of σ u implies that the effective surface energy varies with temperature. Determination of temperature-dependence of γ eff relies on the implementation of appropriate cleavage fracture mechanisms. In this study, the multi-barrier micromechanics model (Lin et al., 1987;Martín-Meizoso et al., 1994;Pineau et al., 2016), which can well describe the explicit cleavage crack propagation process, has been chosen for estimating the effective surface energy. According to the multi-barrier model (Lin et al., 1987;Martín-Meizoso et al., 1994;Pineau et al., 2016), cleavage fracture occurs following three elementary stages: microcrack initiation within particles, microcrack penetration over the particle/matrix interface, and microcrack extension through the grain boundary. The first barrier and the second barrier, i.e. the particle/matrix interface and grain boundary, are governed by two respective effective surface energies γ pm and γ mm (Pineau et al., 2016).
The effective surface energy γ pm is usually regarded as temperatureindependent (Kawata et al., 2018) and can be measured from the relationship between local fracture stress and the reciprocal square root of particle size in terms of the Griffith theory (Bowen et al., 1986;San Martin and Rodriguez-Ibabe, 1999). It varies normally in a range of 7-20J/m 2 for steels, depending on the microstructure and method of measurement (Kawata et al., 2018). These measurements were normally performed with specimens tested at very low temperatures, e.g. 77 K-143 K, where pure cleavage occurs accompanied with very limited plastic deformation. However, the effective surface energy γ pm may increase with the increase of temperature as more plastic work dissipates prior to the penetration of a microcrack over the particle/matrix interface. Enlightened by the shielding effect of dislocation dynamics on crack tip (Lin and Thomson, 1986;P. B. Hirsch, 1989), a method has been established in (Li et al., 2019) to estimate the temperature-dependence of γ pm in the transition region for a ferritic steel.
The effective surface energy γ mm is known to be strongly temperature-dependent (Lin et al., 1987;Linaza et al., 1995Linaza et al., , 1997Wallin et al., 1984). Following the Griffith theory, a temperature-dependent γ mm has been established for a Ti/V steel by San Martin et al. (San Martin and Rodriguez-Ibabe, 1999) through calculating the size of cleavage islands surrounded by ductile fracture and fracture stress required for the microcrack propagation across a grain boundary. However, the γ mm was overestimated since the sizes of individual cleavage islands are likely to be larger than the statistically descriptive average cleavage facet size (d CFS ), which governs the cleavage crack propagation. This is the reason for the extremely rapid increase of the estimated γ mm , e.g. from 50J/m 2 at 173 K to 500J/ m 2 at 223 K (San Martin and Rodriguez-Ibabe, 1999). A less temperature-dependent γ mm has been evaluated according to the Griffith theory by Kawata et al. (2018). In their work, the mixed-mode stress intensity factor near the grain boundary is calculated by measuring the local stress and the orientation of a neighboring cleavage facet compared with the cleavage nucleation facet. Since cleavage fracture occurs abruptly, the measured local fracture stress can hardly be identified as the critical stress for microcrack propagating over the particle/matrix interface or that of a grain boundary. This could lead to an overestimation of γ mm at low temperatures because the critical stress for a microcrack penetration into a grain is expected to be higher than that for the crack trespassing grain boundaries (Li et al., 2019;Lin et al., 1987;Linaza et al., 1997).
In a preliminary study (Li et al., 2019), the effective surface energy concept was incorporated into the CAFE method to predict the DBT of a ferritic steel reported in the literature. Due to the lack of necessary microstructural data, various nonphysical approximations were made. In this study, an attempt is made to formulate a complete methodology for predicting the DBT through coupling the CAFE method and the multi-barrier model. The predicting power of the proposed methodology is verified by the application to a low-carbon ultrahigh-strength steel. Within the proposed methodology, the competition between ductile damage and cleavage fracture, and the probabilistic nature of microstructures are intrinsically represented by the CAFE method. With the multi-barrier model, physically-based variables can be determined to describe the variation of fracture toughness in the transition region. A continuum method is formulated to quantify the temperature-dependent γ pm for a microcrack penetration over the particle/matrix interface. Two scenarios to identify the intermediate parameters for the estimation of γ pm are adopted and discussed. A temperature-dependent γ mm is obtained by a method established to represent the competition between particle-size controlled and grain-size controlled unstable cleavage crack growth. The effective surface energy γ mm is estimated from γ pm and measured microstructure features related to the cleavage crack propagation. Characterization of microstructure, mechanical tests, and fractographic analysis of the low-carbon ultrahigh-strength steel are performed to capture the essential parameters for the prediction of DBT.

The modelling methodology
The methodology explored in the present work that combines the CAFE method and the multi-barrier model to predict the DBT of a steel is schematically illustrated in Fig. 1. It contains three key elements: (i) the competition between ductile damage and cleavage processed by the CAFE method, where the probabilistic nature of microstructures is incorporated; (ii) a physically-based variable γ pm that describes the second stage of cleavage propagation is estimated by using a continuum model based on the shielding effect of dislocation mobility at crack tip; (iii) a correlation between γ mm and γ pm established based on the likelihood between particle-size controlled and grain-size controlled cleavage fracture in transition region to depict the occurrence of the third stage of cleavage propagation.

The CAFE method
In the CAFE method (Fig. 1), a set of cellular automata (CA) arrays are designed to reflect the mechanical essentials of a microstructure and handle the corresponding evolving failure of the material, while the stresses and strains are processed by finite elements (FE) (Shterenlikht, 2003;Shterenlikht and Howard, 2006). To reproduce DBT with two micro-mechanisms operating simultaneously, two distinct cell-confining cubic arrays, i.e. ductile and brittle, are mathematically assembled into the material integration points to represent relevant microstructural features in a finite element. The damage variable and strain increment tensor are distributed to arrays in the material integration point, by which the ductile damage or cleavage is treated within ductile or brittle array. At the structural level, the failure of an element is decided via a so-called integrity indicator Σ, which first gathers the cells' evolution, i. e. manifesting the progress of materials' failure, and then is fed back to the finite element. An explicit dynamic algorithm using VUMAT is chosen for the implementation of CAFE into the finite element program ABAQUS, on which crack growth following a natural failure path is accomplished in the numerical simulation of the Charpy test through element removal (Shterenlikht, 2003;Shterenlikht andHoward, 2004, 2006). The Rousselier model (Rousselier, 1987) is used to describe the constitutive response in ductile cell array. Its yield function is expressed by where σ 1 and D are material constants; H(ε eq ) is the material's hardening property; ε eq is the equivalent strain; σ eq and σ m are equivalent stress and mean stress; B(β) is the function of damage variable β that is relevant to plastic deformation and damage; and ρ is relative density, a function of the initial void volume fraction f 0 , which can be expressed by The spacing of primary voids which denotes the damage unit of a material, is taken as the ductile cell size L D . In the Rousselier model, the ductile failure is assumed to occur when the damage variable β reaches its critical value β F .
With an assumed penny-shape defect ahead of a macroscopic crack tip in a brittle cell array, the fracture strength of a brittle cell is calculated by where γ eff is the effective surface energy for cleavage fracture; E and ν are Young's modulus and Poisson's ratio, respectively; and d g is the grain size. In this study, the cleavage facet size d CFS , measured from fractographic analysis, is adopted as the size of brittle cells L B . To consider microcrack nucleation, a fraction of brittle cells η is assumed to describe the grains with particles where a microcrack has already nucleated (Shterenlikht, 2003;Shterenlikht and Howard, 2006). Since high-angle grain boundaries can naturally resist cleavage crack propagation (Stec and Faleskog, 2009), an arbitrary orientation is allocated to every brittle cell, where a maximum value is assumed to be 70 • . Furthermore, a threshold θ th needs to be set for crack propagation from one cell to the other. A cleavage crack propagates when the local stress in a brittle cell exceeds σ F and the misorientation between a cell and its neighboring cell is lower than θ th .
In both cell arrays, the state of a cell at next time step depends on the To calculate the fracture strength of a brittle cell, γ mm is calculated through a correlation with γ pm (g), in which γ pm is estimated by using a continuum model (h). By applying stress-strain behavior at different temperatures (d) into the finite element model (e), the load-displacement response and DBT of the Charpy test can be produced using the CAFE method implemented with the estimated γ mm (i).
state of itself and the neighboring cells' state at the current time step (Shterenlikht, 2003). The simulation of the damage development in ductile arrays or fracture propagation in brittle arrays is achieved by changing the state of failed cell from the initial state of 'alive' to 'dead'. Since 'dead' cells lose their load-bearing capacity, the cells surrounding a 'dead' cell will carry an elevated local load, which induces a local strain/stress concentration. The CAFE method solves this by using concentration factors C D for ductile cells and C B for brittle cells. In addition, a solution has been proposed to let CAFE method automatically locate such a close-neighboring cell in the vicinity of the 'dead' one (Shterenlikht, 2003;Shterenlikht and Margetts, 2015). The concentrated cells' state, Ω ti+1 m , are determined by Eq. (4-1) and Eq. (4-2) for a ductile and brittle cell at the next time increment, respectively.
where β ti m is the damage variable of a strain-concentrated cell at current time increment; σ ti m is the maximum principal stress on the stressconcentrated cell at current time increment; and θ m is the orientation at the stress-concentrated cell. Due to the cumulation of damage, the integrity indicator Σ with an initial value 1.0 reduces successively before N D or N B attains its maximum value N max D or N max B , N D and N B being the numbers of dead cells in ductile and brittle arrays, respectively. At this time, Σ becomes zero, and the material represented by the integration point is assumed to be failed. The associated element will be removed once it receives the message of zero Σ. The integrity indicator Σ is calculated by where N D and N B represents the quantity of 'dead' cells in ductile and brittle array, respectively. Because two arrays take the position of the same physical space, once a cell in an array is failed, the cell (cells) in another array is assumed to be failed too. This enables the material inside the integration point to fail in a single manner, i.e. ductile or brittle. Therefore, the competition between ductile damage and cleavage dominates the fracture at structural level because the two failure mechanisms are completely independent. A mapping rule is used to synchronize the evaluation of cells within two arrays, which allows the cells failed in a ductile array to be projected into the relevant brittle array, and vice versa (Shterenlikht, 2003). The CAFE method has been presented in more detail in (Shterenlikht, 2003;Shterenlikht and Howard, 2006).

Determining the surface energy for microcrack penetration over particle/matrix interface γ pm
Failure of material is assumed to depend on the result of the competition between atomic debonding (cleavage) and the dislocation nucleation and emission (plasticity) (Rice and Thomson, 1974). Dislocations can cause a shielding effect on crack tip within a finite region, which is the origin of material's toughness against crack propagation (Thomson and Sinclair, 1982). The shielding effect of dislocation dynamics, i.e. dislocation nucleation and emission, on crack tip has been evaluated through a continuum model (Hartmaier and Gumbsch, 2005;Li et al., 2018;Nitzsche and Hsia, 1994). The DBT of iron and steels, i.e. the body centered cubic metals, is induced by the thermally activated dislocation dynamics (Argon, 2000). Based on the same mechanism, the continuum model designed for studying DBT of single-crystal iron (Li et al., 2018) has been extended to ferritic steels (Li et al., 2019). Under the same circumstances, the continuum model is further extended to a martensitic steel in this study.
An upper half of the symmetric boundary layer model is presented in Fig. 2. It is assumed that the elastic particle with a radius a occupies the center of the model, e.g. a small circle around the crack tip. The small particle is completely enclosed by the matrix, e.g. martensite, within the model. In the matrix, the material is assumed to be time-dependent plastic, which is adopted for describing plasticity produced by dislocations. Inspired by the work of (Hartmaier and Gumbsch, 2005), the radius R of model is designed with a size 20 times larger than that of a particle. A pre-crack with a tip-radius of 1.15 × 10 − 4 R is initially introduced into the model to represent the cleavage initiation within the particle. The applied stress intensity factor K a І is implemented into the modified boundary layer model with a loading rate K by controlling the outer-layer nodal displacement.
In this continuum model, the time-dependent plasticity due to the dislocation mobility is described by the equivalent plastic strain rate ε (Li et al., 2018(Li et al., , 2019 where k B is Boltzmann constant; ε 0 is the reference strain rate; Q is the activation energy for dislocation velocity; m is a temperature-dependent stress exponent for a wide range of stress level (D. Hull, 2011); T is temperature in Kelvin; σ e is the von Mises equivalent stress; and σ 0 is a normalization stress. It is assumed that failure of material occurs once the local stress intensity factor at crack tip K t І reaches its critical valueK ІC , which is determined by the material's surface energy γ s for cleavage fracture. K t I is normally lower than K a I , particularly at high temperature, since stresses within the particle could be relaxed by plastic deformation around the crack tip (Li et al., 2018). When the failure occurs, K a І can be considered as the material's fracture toughness. Following the modified Griffith theory G c = 2(γ s + γ p ) = 2γ eff , the effective surface energy is calculated as (Irwin, 1948) The effective surface energy γ pm for microcrack penetration over the particle/matrix interface can be estimated according to Eq. (7).

Determining the effective surface energy for microcrack propagation across the grain boundary γ mm
For the formation of unstable cleavage fracture, the second stage of cleavage crack propagation in the multi-barrier model is a critical step (Kroon and Faleskog, 2005). The critical fracture stress for microcrack penetration into the matrix, namely the particle cleavage strength σ pm (Lin et al., 1987), is expressed by where c crit. is the particle size of cleavage initiation. When unstable cleavage fracture is initiated, the grain boundary it encounters first can be the obstacle for cleavage fracture propagation. At the moment when the microcrack propagates across the grain boundary, the critical fractures stress, so-called grain strength σ mm (Lin et al., 1987), can be written Equations (8) and (9) indicate that cleavage fracture propagation is essentially governed by the competition between σ pm and σ mm (Li et al., 2019), i.e. particle-size controlled or grain-size controlled. A criterion has also been deduced for the microcrack either being arrested or propagating across the first grain boundary (Li et al., 2019), upon which the lower limit of the γ mm can be estimated by This suggests that a minimum value of γ mm can be obtained from γ pm and the ratio of d g /c crit. , in which the grain size d g should be replaced by d CFS in Eq. (10) since the latter is better suited to characterize the crystallographic unit of cleavage (Kim et al., 2003). d CFS is assumed to be temperature independent (Li et al., 2019). In addition, an artificial linear relation between temperature and c crit. has been employed to estimate γ mm in Eq. (10) (Li et al., 2019), which was inspired by the study by Lee et al. (2002). However, to accurately estimate γ mm , the variation of both d CFS and c crit. with temperature should be measured for a specific material. Recall Eq. (3), γ mm estimated here is exactly the one required to assess the critical cleavage fracture stress. In the following, γ mm will be integrated into the CAFE method to predict the DBT.

Material and experiments
The studied material is a S960 grade laboratory hot-rolled asquenched steel. Its chemical composition is 0.13C-1.1Mn-0.7Cr-0.16Mo-0.03Ti (wt.%). The cut cast pieces were first homogenized at 1050 • C for 2 h, hot-rolled at temperatures from 1050 • C to 850 • C with a total reduction of 80% to a final thickness of 12 mm, and directquenched in water to room temperature. These plates were then further reaustenitized at 870 • C for 20 min, and quenched in water to obtain a homogeneous and thoroughly hardened as-quenched lathmartensitic microstructure with an equiaxed prior-austenite grain structure (Pallaspuro et al., 2019), Fig. 3. The average grain size is ~1 μm, and the grain size distribution shown in Fig. 4 is described by a three-parameter Weibull distribution. Its scale, shape, and location parameter are 0.00044, 0.125, and 0.391, respectively.
Tensile tests were conducted on 10 mm-thick rectangular specimens at room temperature using a MTS 250 kN test machine and following EN ISO 6892-1 (ISO, 2009b). The room-temperature yield strength is 963 ± 18 MPa, tensile strength 1252 ± 14 MPa, and uniform elongation 3.6 ± 0.2% (Pallaspuro et al., 2019). To obtain the stress-strain behavior of other temperatures based on that of the room-temperature, the yield strength at each temperature is first calculated by using the relation reported in FITNET (FITNET, 2008). Next, the hardening laws of the material at different temperatures are acquired by power-law fitting of the smoothed tensile curves of a similar A514 AQ steel tested at corresponding temperatures (Partin et al., 2010). Then, the true stress-strain curves at different temperatures can be acquired as shown in Fig. 5(a). For the study material, at room temperature the strain hardening exponent n = 0.086 and the strength coefficient K = 1644 MPa are obtained, and at − 100 • C the fitted n = 0.081 and K = 1753 MPa are  obtained as well. Charpy V-notch (CVN) impact tests were conducted with 10 mm × 10 mm × 55 mm L-T oriented specimens at temperatures between 193 K-313 K according to EN ISO 148-1 (ISO, 2009a). The experimental impact toughness at different temperatures is presented in Fig. 5(b). The DBT curve is obtained by fitting the test results following a method described in (Wallin, 2011a), Fig. 5(b). The upper shelf energy (USE) is 205 J, and the DBT temperature T 50 is 254 K. Here, T 50 is the temperature where the absorbed energy equals to the half of USE reduced with the lower shelf energy (LSE) (Ghosh et al., 2013).
Further microstructural and fractographic characterization is carried out with ZEISS Sigma field-emission scanning electron microscope equipped with both EDAX Hikari XP electron backscatter diffraction (EBSD) camera and TSL OIM software and EDAX Apollo X energy dispersive spectroscopy (EDS) detector and Team software. All the measurements of the characteristic microstructural and fractographical units are approximated with either a circular, square, or triangular shape, and the resulting sizes are given as equivalent circle diameters.
The fracture surfaces of CVN specimens tested at higher temperatures, i.e. ~293 K and 273 K, are analysed to measure the relevant microstructural features of ductile damage like void spacing and initial void volume fraction f 0 . A typical dimple morphology is shown in Fig. 6 (a) and a close-up in Fig. 6(b), in which the size of the large dimples (≥~10 μm) and their average distance between two or three closest large dimples are measured. Observable inclusions with sizes of few micrometres are still presented at the bottom of the large dimples, Fig. 6 (b). These are identified mainly as Al and Ca based oxides and sulphides, i.e. AlO, CaS, CaAlO, and CaOS. The dimples on the fracture surface of the specimen failed by ductile damage exhibit a uniform feature on not only the size but also the distance, Fig. 7. Through measuring in total 81 large dimples, an average size of larger dimples D avg.
dim. of 17.1 μm is obtained, and their average distance L avg.
dim. of 34.3 μm is approximately double of that, which is used as the void spacing and ductile cell size L D , as mentioned before. The f 0 is determined as the area of these inclusions that have nucleated the large dimples divided by the total area in the images of ductile fracture. The resulting void volume fraction is 0.0025, as measured from a total area of ~114 600 μm 2 and with an average inclusion diameter of 3.3 ± 1.1 μm of 35 recorded inclusions.
Fracture surface analysis indicates that normally only one broken    6. (a) The appearance of a ductile fracture surface on a specimen tested at room temperature, (b) a close-up to the dimples with an illustration how the void size and spacing are determined.
coarse particle located at a certain distance from the pre-crack or notch root dominates the origination of cleavage fracture, which accords with the weakest link model (Linaza et al., 1997;San Martin and Rodriguez-Ibabe, 1999). Thus, through tracing the river patterns, the cleavage fracture initiation sites, i.e. the critical broken particles adjacent to sufficiently large cleavage facets, can be identified on the fracture surfaces in regions narrowed down with the aid of an optical stereo light microscope. On the fracture surface of a CVN specimen tested at 193 K, a small region beneath the notch root is analysed, Fig. 8(a). Fig. 8 (b) shows a larger magnification of the area where brittle failure has originated, in which a typical river-pattern of cleavage is visible, and a broken CaAlO particle is located as a cleavage fracture initiation site. In the whole data set of 18 CVN specimens, CaOS or CaAlO based inclusions or clusters of them and TiN inclusions/clusters were the most frequent ones to nucleate the failure. The sizes of these inclusions are measured on the fracture surfaces of 14 identified cases on samples tested between 193 K and 253 K. The measured results are presented in Fig. 9(a), in which the size of critical broken particles c crit. increases with the temperature.
Cleavage facets develop within grains during the propagation of a brittle microcrack, which is originated from a critical particle. The identification and measurement of the first cleavage facets are performed to determine the cleavage facet size d CFS , Fig. 8(c). Due to the complexity of the grain-refined martensitic microstructure and its very fine effective grain size, it is in most cases ambiguous or impossible to determine the sequence and/or only a single cleavage facet that would have caused cleavage crack propagation to an extent of a global failure.
In these cases, the primary cleavage fracture size is determined as the size of the larger facets immediately surrounding the cracked particles that have caused the specimen failure, Fig. 8(c). These primary cleavage facets are critical with respect to the third stage of cleavage propagation, since they have a small misorientation with the cleavage initiation site (Linaza et al., 1997). The primary cleavage facets are identified on the fracture surface of the CVN specimen tested at 233 K and 213 K as shown in Fig. 9(b). The primary d CFS shows a temperature-independent feature, in which the average primary cleavage facet size d avg.
Since cleavage fracture occurs along the {100} crystallographic planes, the deflection or arrest of a crack connects only to the misorientation of a grain boundary, which is usually temperature-independent. Then, we can plot the cumulative probability of d CFS to reveal the statistical characteristic of the cleavage fracture units in the material, Fig. 9(b). The size of larger cleavage facets in the cumulative distribution at the 90% probability level d 90% CFS is approximately 13.8 μm, Fig. 9(b). As the cleavage initiation is based on sampling a critical weakest link in the stressed microstructural volume, such an approximation of an 80% or 90% probability level is justifiable, and has sufficient experimental validation in the literature (Isasti et al., 2014;Kaijalainen et al., 2013;Pallaspuro et al., 2018).
To identify the high-angle grain boundaries able to stop crack propagation in a Charpy test in the given material, the misorientation of the grain boundaries where secondary cracks arrest was measured with EBSD from immediate vicinity of the main crack front from the centreline crosscuts of CVN samples tested at 193 K, Fig. 10. An average of  arrest of 48 • is obtained from the measured 34 independent secondary cracks, which will be adopted as the threshold of the misorientation for cleavage propagation as mentioned in Section 2.1. All results obtained from the fractographic analysis on the CVN specimens are concluded in Table 1.

Prediction of the DBT
To reproduce the transient response of the Charpy test specimens, an explicit dynamic process is modelled using Abaqus 14.1, where a user subroutine VUMAT is used to implement the CAFE method. Since void spacing and d CFS measured on the material are so small (measured in Section 3), the memory required to compile the code exceeds the upperlimit of the available super computer. To mitigate this problem, a subsize 5 mm-thick CVN specimen with a notch radius 0.25 mm and a notch depth 2.0 mm is adopted following EN ISO 148-1 (ISO, 2009a). To model the Charpy-V notch test, the elastic body of striker and two anvils is meshed with C3D8R and C3D6 type elements, Fig. 11(a). The elastic-plastic body of specimen is meshed with C3D8R type elements. Since fracture is anticipated to happen in a narrow region in the middle of a specimen, i.e. in the damage zone, cells that handle damage are assigned only into the elements of this region with a mesh size 0.25 mm, Fig. 11 (b). The model for the CVN specimen consists of total 30600 elements while there are 2320 elements located in the damage zone. The initial velocity of the striker is 5.5 m/s. A friction coefficient 0.15 is utilized to model the contact between the CVN specimen and the    striker/anvils. The predicted absorbed energy of the 5 mm-thick specimen is converted into that of a full-size specimen according to the methodology proposed by Wallin et al. (Wallin, 2011a;Wallin et al., 2016).

The effective surface energy for microcrack penetration over particle/ matrix interface γ pm
It should be noted that the plasticity induced by dislocation dynamics around the crack tip depends on loading rate, Eq. (6). An approximate loading rate K = 5.14 × 10 4 MPam 0.5 s − 1 is obtained from J *, which is fitted from the linear section of the relationship between the loading time and J-integral. Here, the J-integral is calculated by simulating the Charpy test on the 5 mm-thick CVN specimen. With the loading rate, γ pm in the transition regime can be calculated according to Eq. (7). *The loading rate K can be determined by dK Here, dJ dt is a constant through fitting the linear section of the relationship between the loading time and J-integral. J is a time variable. Since J-integral in its linear section varies slightly with loading time, an average of J-integral in this linear section is used to estimate the K .
A correlation between the loading rate K and the critical DBT temperature T c has been found through experiments (Hirsch and Roberts, 1991), which can be described by where E a is the activation energy for DBT, and equivalent to Q in Eq. (6). T c is defined as the temperature where brittle fracture shifts to pure ductile fracture (Hartmaier and Gumbsch, 2005;Li et al., 2018;Nitzsche and Hsia, 1994). It indicates that K t I will never exceed K IC at T c . Here, K IC = 2.012 MPam 0.5 is obtained with the effective surface energy of 9J/ m 2 , which was measured for martensitic steels at 77 K (Bowen et al., 1986). Eq. (11) implies that parameters used in the continuum model introduced in Section 2.2 can be identified through comparing the predicted T c to the experimental one under the same loading rates. A minor variation of E a within a range of 0.2-0.5 ev has been found for single-crystal iron, polycrystal iron, and steels (Giannattasio et al., 2007;Tanaka et al., 2009). This indicates that a small difference among parameters identified with E a of distinct steels can be anticipated. Tanaka et al. (2009) have also adopted Eq. (11) to describe the correlation between K and T c of a low-carbon steel based on four point bend tests. E a of a low-carbon steel is adopted to identify the parameters, with which the γ pm of the low-carbon ultrahigh-strength steel can be calculated in present study. For the exploration of a solution to calculate γ pm , the gap between two steels can be possibly ignored. Recall Eq. (6), the plastic deformation due to the dislocation dynamics under specific loading rate depends on the parameters, i.e. ε 0 , Q, T, m, and σ 0 , in which Q is the one measured by Tanaka et al. (2009), the temperature-dependent m can be acquired from (Turner and Vreeland, 1970) and σ 0 is assumed to be a constant. The only parameter to be identified is. ε 0 .
There exist two scenarios to identify ε 0 . In the first scenario S.1, it is assumed that DBT temperature T 50 measured from the Charpy tests can be adopted to characterize the transition from brittle to ductile in a theoretical manner, namely T 50 equivalent to the critical DBT temperature T c . T 50 of Charpy test with the 5 mm-thick specimen, i.e. T 5mm 50 , can be acquired by subtracting 20 K from the transition temperature T 50 = 254 K measured for the Charpy test of 10 mm-thick specimens in Section 3. Note that the 20 K offset is applied when estimating the transition curve of 5 mm-thick specimen from that of 10 mm-thick specimen according to the methodology proposed by Wallin et al. (Wallin, 2011b;Wallin et al., 2016). Accordingly, a target critical DBT temperature for the Charpy test with the 5 mm-thick specimens, i.e. T 5mm c , can be set initially as 234 K. The continuum model is applied to predict the T 5mm c with the loading rate of 5.14 × 10 4 MPam 0.5 s − 1 . If it equals to the target value 234K, parameters including ε 0 can be identified. Further, these identified parameters are adopted to predict the brittle-to-ductile transition curves for different loading rates so as to build the correlation between K and temperature according to Eq. (11). Then, parameters, particularly ε 0 , can be verified if the correlation between K and temperature predicted by the continuum model is similar to the one reported by Tanaka et al. (2009). In the second scenario S.2, it is assumed that the DBT of studied material, i.e. obtained by Charpy test on 5 mm-thick specimens, can be characterized by the T 5mm c predicted by using continuum model with identified parameters and the same loading rate. The parameters including ε 0 are optimized and identified through directly reproducing the experimental correlation between K and temperature of a low-carbon steel measured by Tanaka et al. (2009). It can be seen that S.1 is basically a semi-experimental parameter identification approach since T 5mm 50 , i.e. the target T 5mm c , is obtained from the experimental results of the identical material. On the other hand, S.2 is an approximate approach on the assumption that the predicted T 5mm c can characterize the transition of the low-carbon ultrahigh-strength steel when the 5 mm-thick specimens are used.
To compare the predicted correlation between K and T c with the experimental one, following the procedure introduced in (Li et al., 2018(Li et al., , 2019, the strain rates applied on the four-point bend test are converted into the loading rate K . Under these specific loading rates, T c is predicted through adopting the continuum model (Section 2.2). It is found in Fig. 12(a) that the slopes of the fitted lines produced by the two scenarios are comparable to the experimental result. This implies that the parameters listed in Table 2 are capable of producing a similar relationship described by Eq. (11). In addition, different intercepts of the two fitting lines based on the two scenarios can be found in Fig. 12(a). In (Li et al., 2019), a constant m has been used to calculate γ pm of a TMCR steel, which led to an under-predicted fracture toughness. To improve this, the exponent m varied with temperature is employed in this study. It has been shown that the predicted T c and fracture toughness are slightly influenced by the elastic zone or particle size of the continuum model ( Fig. 2) (Li et al., 2019). In this study, particle size a = 1 μm has been adopted in the later simulations.
Fracture toughness of the low-carbon ultrahigh-strength steel in DBT region is predicted for Charpy impact test with the 5 mm-thick specimen by applying the parameters identified (see Table 2). Following Eq. (7) and the two scenarios, two different effective surface energy γ pm are achieved, Fig. 12(b). In addition, two fitted curves for the temperaturedependent γ pm are presented in Fig. 12(b).

The lower limit of effective surface energy for microcrack propagation across grain boundary γ mm
According to the method described in Section 2.3, the lower limit of γ mm can be estimated from γ pm in Fig. 12(b) and the ratio of d CFS /c crit. .

Here, d CFS is replaced with d 90%
CFS to represent the cleavage fracture unit. c crit. is measured on the fracture surfaces of CVN specimens tested at temperatures 193 K-253 K. Therefore, the ratio of d 90% CFS /c crit. is calculated for each temperature, Fig. 13(a). Since the larger particles are more critical than the smaller ones due to lower local stress required for the microcrack penetration over the particle/matrix interface, the upper limit of the correlation between d 90% CFS /c crit. and temperature will be adopted for the estimation of γ mm . The calculated γ mm is presented in Fig. 13(b) based on γ pm calculated with two different scenarios.

The DBT of the low-carbon ultrahigh-strength steel
The CAFE method implemented with γ mm shown in Fig. 13(b) are employed to simulate the Charpy test with the 5 mm-thick specimens (Fig. 11). The load displacement curves and absorbed energies at different temperatures are obtained. The absorbed energies of the 10 mm-thick specimens at all temperatures are converted following the methodology proposed by Wallin et al. Before doing this, the parameters in the Rousselier model to describe the ductile damage need to be calibrated by reproducing the experimental value of the upper shelf absorbed energy at 313 K. In this phase, the material constants σ 1 and D, and the critical damage variable β F , are calibrated for two different γ mm . The other parameters for ductile damage, L D and f 0 were identified in Section 3. The ductile damage variables are listed in Table 3.
The parameters for cleavage fracture have been measured through the microstructure characterization and fractographic analysis described in Section 3. These parameters are listed in Table 4.
In the CA phase, the number of cells is determined by the element size, ductile cell size, L D , and brittle cell L B . For the material considered, the cell size per linear dimension for ductile (m D ) and brittle (m B ) CA, are set to be 7 and 18, respectively. Each material integration point consists of 343 ductile cells or 5832 brittle cells. The total number of cells implemented into the damage zone is 795760 for ductile and 13530240 for brittle array. Each ductile cell is given a critical damage variable β F chosen according to the distribution of β F listed in Table 3.
Based on the measured grain size distribution presented in Fig. 4 and Table 4, a grain size is allocated to each brittle cell. Once cells within an orthogonal plane of the cell array become dead, the ductile or brittle array is presumed to lose its ability to carry load (Shterenlikht, 2003). Accordingly, the maximum amount of dead cells is N D− max = m 2 D = 49    Note: a three-parameter Weibull distribution, y = r a , is used to depict the grain size (d g ) distribution in the material, where a, r and u are scale, shape and location parameter, respectively.
for a ductile array and N B− max = m 2 B = 324 for a brittle array. When a cell is dead in either array, concentration factors c D and c B are adopted to characterize the strain/stress concentration in the vicinity of a ductile cell and a brittle cell, respectively. Here, c D and c B are set to be 1.4 and 3.5, respectively. In addition, a fraction of cells in each brittle array η = = 5.1 × 10 − 4 is adopted to denote grains carrying particles where a microcrack has already nucleated. It is noted that c D , c B and η are obtained from the calibration of ductile damage parameters for two different γ mm . Fig. 14 presents the variation of predicted absorbed energy with temperature, where the test result shown in Fig. 5(b) is also included for comparison. The scatter in impact toughness is reproduced in transition region by running the CAFE calculations three times at each temperature. The reproduction of the scatter is possible because of the possibility to incorporate the experimentally measured statistical features of microstructure, i.e. β F , d g , and orientations of cells, see Tables 3 and 4, respectively. With the two designated scenarios, the predicted impact toughness in full temperature range can capture the experimental DBT behavior of the low-carbon ultrahigh-strength steel with comparable LSE and USE. The predicted impact toughness in the transition region is within both upper-bond and lower-bond experimental values. Although the predicted absorbed energies are overestimated at 273 K, the prediction based on S.1 is more accurate than that based on S.2, since the predicted absorbed energies of S.1 match the experimental results better at 253 K and 233 K.

Discussion
The CAFE method has clear advantages in solving the long lasting challenges in modelling the DBT of steels. There is still one important benefit that needs to be emphasized. In the brittle cell array, the nucleation of cleavage, the propagation of existing crack in cells, and even the crack arrest can be modelled by the CAFE method . It can numerically reproduce the full process of explicit cleavage propagation, comparing with the local approach, e.g. the Beremin model, where the failure of whole structure is evaluated only as a failure probability .
It is obvious that the temperature-dependent flow stress (Li et al., 2019;Shterenlikht and Howard, 2006) is needed for the prediction of the DBT. However, the CAFE method itself cannot properly capture the DBT of steels without a second temperature-dependent variable. In the present study, an attempt is made to search for a physically-based variable, i.e. the effective surface energy varied with temperature.
According to the multi-barrier model, the effective surface energy associated with cleavage fracture contains two substances, i.e. γ pm and γ mm . With a continuum model based on the shielding effect of dislocation dynamics on crack tip, the temperature-dependent γ pm of the low-carbon ultrahigh-strength steel is theoretically estimated to describe the plastic work dissipated during the cleavage formation (i.e. microcrack nucleation and penetration over the particle/matrix interface) at different temperatures. The temperature dependence of γ pm is different to the conventional treatment, i.e. temperature-independent γ pm . There are two possible reasons for the conventional conception: (i) the measured γ pm is usually obtained from the relationship of local fracture stress and the reciprocal square root of particle size, which is usually measured at very low temperatures of 77 K-143 K (Bowen et al., 1986;San Martin and Rodriguez-Ibabe, 1999). Almost pure cleavage occurs accompanied with only minor plastic deformation at these low temperatures. (ii) The plastic work γ p occurred during the process of cleavage formation can hardly be measured, and it is also unlikely to separate the γ p from the surface energy γ s experimentally. The continuum method adopted in the present work provides a feasible approach to quantify the required γ p for the cleavage formation and capture the energy barrier to be overcome in the second stage. Linaza et al. (1997) and San Martin et al. (San Martin and Rodriguez-Ibabe, 1999) have presented a method to estimate γ mm . However, certain ambiguity exists in their practice, for example, the usage of a constant γ pm and a cleavage island size to describe the cleavage fracture unit. The alternative method to quantify γ mm proposed by Kawata et al. (2018) tends to overestimate γ mm at low temperature. To solve these problems, a method is proposed to estimate the lower limit of γ mm in present study. A fact has been revealed that microcrack trespassing or being arrested at the grain boundary relies on the competition between the particle-size controlled and grain-size controlled cleavage fracture propagation (Li et al., 2019). Further, the mechanism of this competition is represented in this method through establishing the correlation between the critical fracture stress for cleavage propagation across the particle/matrix interface and the one over the grain boundary. An explored temperature-dependent γ pm is introduced. In addition, d CFS is adopted to accurately characterize the cleavage fracture unit. Particularly, a temperature-independent d CFS is obtained experimentally, and a real temperature-dependent c crit. is measured. These two variables to describe the microstructural features of cleavage fracture are the key elements to build up a tailored method to estimate γ mm for a study material.
In order to properly predict the DBT of the studied material, it is important to acquire reliable parameters needed in continuum model to calculate the temperature-dependent γ pm . Following Eq. (11), this can be achieved by reproducing the correlation between K and T c obtained experimentally, particularly with a similar E a . A measured E a of the lowcarbon steel reported in (Tanaka et al., 2009) is approximately adopted as that of the low-carbon ultrahigh-strength steel studied in present work, which is the basis of the parameters identification. Although both scenarios S.1 and S.2 produce similar E a with the experimental one, i.e. the slope of the correlation between K and T c , S.1 can capture the experimental DBT more accurately than S.2, Fig. 14. This implies that the semi-experimental scenario of parameter identification is probably a better alternative to calculate γ pm .
A systematic characterization of an as-quenched lath-martensitic steel was performed to retrieve the microstructural features relevant to the mechanical properties, especially the statistical features. The results presented in Fig. 14 indicate that the methodology explored in the present work is capable of reproducing the transition behavior. Comparing to the preliminary study (Li et al., 2019), an evident improvement on the accuracy of prediction of DBT has been achieved through introducing microstructural parameters, i.e.d 90% CFS /c crit. , θ th , and a semi-experimental scenario to identify the parameters. However, there are still some limitations in present methodology, like the compromise of using experimental E a of a different material from the literature due to the absence of the measurement of E a on the low-carbon ultrahigh-strength steel. Such a measurement will be needed to get a more accurate prediction of DBT of the studied material in the future. Finally, in modelling the Charpy test the adiabatic heating effect and viscoplasticity are omitted in Rousselier model for the material studied, which should also be taken into account in future studies.

Conclusions
Based on the proposed methodologies, one can concluded that following material data are needed to predict ductile-brittle transition in bcc material, as demonstrated here with an ultrahigh-strength steel: (1) material stress strain curves at different temperatures, (2) the activation energy for DBT, E a , (3) d CFS /c crt. measured at the transition region, (4) statistical features of microstructure of the material, i.e. f 0 , void spacing L D , grain size, cleavage facet size d CFS , and the threshold of misorientation θ th of grains where cracks arrest locally.
To demonstrate the application of the proposed methodology, systematic microstructural characterization and also fractographic analysis should be performed to capture the statistical features of a given microstructure, here that of a low-carbon ultrahigh-strength steel, including a temperature-independent distribution of the primary d CFS . Out of the two scenarios employed to identify the parameters for estimation of γ pm , the semi-experimental scenario appears to be a better alternative.

Declaration of competing interest
To the best of our knowledge, the named authors have no conflict of interest, financial or otherwise.
And, (i) no support, financial or otherwise, has been received from any organization that may have an interest in the submitted work; (ii) there are no other relationships or activities that could appear to have influenced the submitted work.