A Computationally Efficient Multiscale, Multi-Phase Modeling Approach Based on CPFEM to Assess the Effect of Second Phase Particles on Mechanical Properties

Abstract

Crystal plasticity finite element (CPFEM) modeling of metals that can be age hardened consisting of second phase particles is extensively performed based on representative volume element (RVE) models. The RVE model is generated for ferritic low-carbon steel using the data obtained from microstructural observation through optical microscopy (OM) and electron backscatter diffraction (EBSD). The generated RVE is required to statistically represent the original material in terms of grain topology and texture in microscale, as well as the configuration of second phase particles in submicron scale. The multiscale, multi-phase nature of the generated RVE leads to a computationally expensive modeling procedure. To overcome this issue, an alternative multiscale modeling approach based on a homogenization scheme is proposed, in which the effect of second phase particles on deformation behavior is accounted for with no need for the explicit presence of particles in RVE. Lastly, a thorough parametric analysis is performed to investigate the sensitivity of the mechanical properties to the second phase particles in terms of size, volume fraction, geometrical distribution, and deformable or non-deformable properties of precipitates in the investigated material. The results show that the proposed multiscale modeling approach successfully accounts for the effect of second phase particles on deformation behavior, while the computational cost is reduced by more than 99%. In addition, the simulations show that the configuration of second phase particles at a microscale plays an important role in defining the mechanical behavior of the material.

1. Introduction

Metallic materials that can be age hardened possess superior mechanical properties when compared to alloys that cannot be age hardened [1]. The superior performance is related to the microstructures that consist of secondary phases of precipitates distributed over the primary phase matrix. One of these materials is carbon steels, which have been used extensively in various industries, and thus attracted widespread attention in the scientific community. Carbon steels derive their strength through the formation of iron-carbon compounds [2]. The macroscopic properties of carbon steels are highly dependent on the properties of second phase precipitates as well as their interaction with the matrix phase at lower scales [3,4]. For instance, the size, spacing, volume fraction, geometrical distribution, and mechanical properties of precipitates at the microscale [5], as well as particle interfacial properties, homogeneity of chemical composition, and the lattice misfit strain at the submicron-scale [6,7] could play an important role in defining the properties of the material at the macroscale. Therefore, a careful design of type, size, spacing, and distribution of second phase particles through a precise heat treatment scheme may lead to a wide range of material properties for various applications.

Modern product design procedures mainly rely on the use of simulation tools to optimize and accelerate the design cycle. To that end, advanced computational simulation tools and material models are required to accurately represent the performance of the final products [8]. The desired precision cannot be achieved unless the behavior of second phase particles is carefully considered, even though it is very challenging to characterize and model them [8,9]. Several material models have been proposed to consider the effect of second phase particles on the initial yield strength, anisotropy behavior, and kinematic hardening of the metallic materials [10,11,12,13]. For instance, Tanaka and Mori [10] developed a constitutive model, in which the effect of rigid inclusions on the work hardening behaviors of crystals has been considered. In a study by Barlat and Liu [11], a constitutive material model was formulated to consider the effect of second phase particles on plastic anisotropy and kinematic hardening in a binary Al-Cu containing hard θ՛ particles. In another study, Han et al. [12] developed a material model formulated in the crystal plasticity setting to account for the effect of precipitates with pure elastic behavior on the elastic-plastic deformation behavior of metallic crystals. Furthermore, Bonfoh et al. [13] proposed a micromechanical modeling approach to model the elastic-plastic behavior of a polycrystalline containing non-shearable particles.

In the recent literature, the effect of second phase particles on mechanical behavior has been formulated through CPFEM modeling. For instance, Wang and Jiang [14] adopted CPFEM to model the fatigue failure of 304 stainless steel under cyclic loading conditions, in which the effect of second phase particles was evaluated. Furthermore, Keshavarz and Ghosh [4] developed a multiscale constitutive model to simulate the mechanical behavior of a Ni-based superalloy through CPFEM simulations considering the effects of material characteristics at different scales, including microstructures and second phase particles. Ali et al. [5] assessed the effect of γ՛ precipitate on the mechanical properties of a Ni-based super alloy through performing a CPFEM simulation of hardness experiment, as well as crystal plasticity coupled with the phase-field method for creep experiment. In another study carried out by Li et al. [15], the effect of T₁/T₂ precipitates in an Al-Li alloy on the mechanical behavior as well as the microstructure and texture evolution has been investigated through a CPFEM simulation of the uniaxial compression test.

In this study, a CPFEM simulation is performed on an RVE model of ferritic low carbon steel consisting of a matrix phase and second phase particles to investigate the effect of precipitates on the deformation behavior of materials that can be age hardened. To that end, a microstructure-based RVE model is generated using the experimental data, including the properties and distribution of precipitates in submicron scale, as well as the grain topology and orientation distribution of grains in a microscale. The multiscale, multi-phase nature of the generated RVE makes the model computationally very expensive. To overcome this issue, a novel computationally efficient multiscale modeling approach based on a homogenization scheme is developed for the calibration of a CP model, through which the explicit presence of second phase particles is no longer necessary in the RVE model. The primary goal of the multiscale modeling approach is to develop a computationally efficient technique that improves the CPFEM simulations in multiple ways: (a) A significant reduction in computational cost, (b) avoiding the non-unique CP model parameters, (c) examination of the macroscopic properties of the material through incorporating the experimentally validated data of precipitates, and (d) accounting for precipitation hardening in a 3D RVE model with several hundred grains. Furthermore, a thorough parametric study was performed to investigate the sensitivity of the mechanical properties to the second phase particles in terms of size, volume fraction, geometrical distribution, and mechanical properties of precipitates on the macroscopic properties.

2. Experimental Procedure

In this study, a low carbon steel tube with an outer diameter of 38 mm and thickness of 1.24 mm was investigated. The as-received tube was solution heat treated at 900 °C for 75 min and then air cooled to obtain a fully ferritic initial microstructure and uniform chemical composition of Fe—0.13 wt.% C—0.3 wt.% Mn—0.02 wt.%P—0.15 wt.% Si—0.3 wt.% S. To assess the mechanical properties of tubes, a uniaxial tensile test was performed based on the E8-E8M-09 ASTM standard [16]. To that end, the dog-bone specimens were cut in a longitudinal direction with a gauge length and width of 32 mm and 6 mm, respectively. In addition, a pair of special grips was fabricated to accommodate for the curvature of the tubular samples during tensile testing. The uniaxial tensile tests were performed with an INSTRON machine with a 3 mm/min grip speed at room temperature. Further details on sample preparation and testing of tubular materials can be found in the authors’ previous publications [17,18]. Throughout the testing, a 3D Trilion Quality Systems GOM ARAMIS Digital Image Correlation (DIC) machine was used to record the displacement field, while a sensor installed on the tensile machine was used to record the force data. Then, the synchronized force and displacement data were utilized to generate engineering strain-stress curves [19].

To observe the microstructure of the investigated tubes, the specimens were cut along both the hoop and longitudinal directions, so that the surface of both cross sections from almost the same location could be accessed. Then, the samples were mounted, polished to up to 0.05 μm roughness, and etched by Nital detergent containing 100 mL ethanol and 5 mL nitric acid to observe the microstructure through Optical Microscopy (OM). The collected OM images were analyzed through the Dream.3D image processing toolbox [20] to obtain the size distribution of precipitates. Moreover, the crystallographic texture and grain topology of the tubes were measured by means of electron backscatter diffraction (EBSD) imaging. To that end, the EBSD measurements were conducted by a ThermoScientic^TMApreoLoVac field emission EBSD microscope equipped with an EDAX Hikari XP camera. The raw EBSD data were analyzed by means of MATLAB toolbox MTEX [21] to obtain diverse microstructural data, including distribution of grain aspect ratio, grain size distribution, and the crystallographic orientation. To obtain the grain structure, a misorientation angle threshold of 15° was adopted to reveal the grain boundaries. Further details on sample preparation, testing procedure, and analyzing techniques can be found in the authors’ previous publications [22,23].

3. Crystal Plasticity Model

In this study, a rate-independent crystal plasticity (CP) model is employed, in which the yield function for a single crystal is defined as follows [24]:

$$f\left(\mathit{\sigma }\right)=\frac{1}{m}ln\left\{{{\displaystyle \sum }}_{\alpha =1}^{N}exp\left[m\left(\frac{\left|\mathit{\sigma }:{\mathit{P}}^{\alpha }\right|}{{\tau }_y^{\alpha }}-1\right)\right]\right\}$$

(1)

where σ is the Cauchy stress tensor, ${\tau }_y^{\alpha }$ is the critical resolved shear stress (CRSS) on the slip system α, $m$ is a parameter to determine the closeness to the condition, $\left|\mathit{\sigma }:P^{\alpha }\right|={\tau }_y^{\alpha }$, N is the number of slip systems, and $P^{\alpha }$ is the Schmid tensor for the α-th slip system. Here, tensor variables $\mathit{\sigma }$ and ${\mathit{P}}^{\alpha }$ are represented by bold symbols. In this model, the slip rate (${\dot{\gamma }}^{\beta }$) is defined as [25],

$${\dot{\gamma }}^{\alpha }=\lambda \frac{\partial f\left(\mathit{\sigma }\right)}{\partial {\tau }_y^{\alpha }}=\lambda \frac{\frac{sgn\left(\mathit{\sigma } : {\mathit{P}}^{\alpha }\right)}{{\tau }_y^{\alpha }}exp\left[m\left(\frac{\left|\mathit{\sigma }:{\mathit{P}}^{\alpha }\right|}{{\tau }_y^{\alpha }}-1\right)\right]}{{{\displaystyle \sum }}_{\beta =1}^{N}exp\left[m\left(\frac{\left|\mathit{\sigma }:{\mathit{P}}^{\beta }\right|}{{\tau }_y^{\beta }}-1\right)\right]}$$

(2)

where $\alpha$ and β are the index of summation for slip systems.

In the present CP model, a dislocation density-based hardening rule is adopted to account for the impact of microstructure barriers on the evolution of CRSS. The CRSS, ${\tau }_y^{\alpha }$, is decomposed into three terms as follows [8,25]:

$${\tau }_y^{\alpha }={\tau }_0+{\tau }_{HP}+{\tau }_{forest}^{\alpha }$$

(3)

Here, ${\tau }_0$ is the summation of the initial friction stress and the strengthening contribution of precipitation hardening, which is assumed to be a constant value in the present study. ${\tau }_{HP}$ accounts for the Hall-Petch effect, in which the impact of grain structure on the CRSS is considered. ${\tau }_{HP}$ can be calculated by the following equation [25],

$${\tau }_{HP}=\mu H^{\alpha } \sqrt{\frac{b}{d_{MFP}}}$$

(4)

where μ is the shear modulus, $H^{\alpha }$ is the Hall-Petch parameter, b is the Burgers vector, and $d_{MFP}$ is the size of the dislocation mean free path (MFP), which is a function of grain or subgrain microstructures. Lastly, ${\tau }_{forest}$ considers the effect of dislocation density (ρ) on the CRSS through the following relationship [25],

$${\tau }_{forest}^{\alpha }=\mu b\sqrt{{{\displaystyle \sum }}_{\beta =1}^N{h}^{\alpha \beta } {\rho }^{\beta }}$$

(5)

The evolution of the dislocation density, ${\rho }^{\alpha }$, on the slip system α is calculated based on Kocks’ evolution law as follows [26]:

$$d{\rho }^{\alpha }=\frac{1}{b}\left(\frac{\sqrt{{{\displaystyle \sum }}_{\beta =1}^N {\rho }^{\beta }}}{k_a^{\alpha }}-k_b^{\alpha }{\rho }^{\alpha }\right)\left|d{\gamma }^{\alpha }\right|$$

(6)

In this equation, the generation, evolution, and annihilation of dislocation density are controlled by the $k_a$ and $k_b$ parameters. These parameters mainly control the hardening behavior of the material in the plastic deformation region, which thus should be calibrated against experimental data.

It is worthwhile to note that local effects of geometrically necessary dislocations (GND) and discrete dislocation plasticity are neglected from the present model to improve the computational efficiency. However, the effects of these local mechanisms could play an important role on the plastic deformation characteristics within the submicron scale [25,27]. Therefore, neglecting the effect of GND and discrete dislocation plasticity was compensated for through the statistically stored dislocations (SSD) in the material calibration procedure.

4. Calibration of Crystal Plasticity Model

4.1. Virtual Tensile Testing

The CP model explained in the previous section contains various model parameters that should be determined based on the material being investigated.

Table 1. The material constants used in CPFEM model [22,25].

Parameter	Value
Type	BCC
$c_{11}$	231 GPa
$c_{12}$	135 GPa
$c_{44}$	116 GPa
μ	69.3 GPa
b	2.48 × 10−4 μm

provides a portion of the material parameters that are constants, in which the elastic constants of $c_{11}$, $c_{12}$, and $c_{44}$, as well as shear modulus, μ, and Burgers vector, b, are provided for ferrite with a body-centered cubic (BCC) structure. The rest of the parameters, including ${\tau }_0$, $H^{\alpha }$, $k_a$, and $k_b$, are the calibration parameters that are defined through calibration of the model against the experimental data. Among these parameters, ${\tau }_0$, $H^{\alpha }$, and the initial dislocation density, ${\rho }_0^{\alpha }$, contribute to the initial yield of the material, while the remaining parameters ($k_a$ and $k_b$) control the hardening behavior.

As mentioned earlier, a portion of the CP model parameters are defined by calibration of the CPFEM results against the experimental data. It is known that the proper calibration process requires an adequate set of experimental data to achieve accurate and unique CP model parameters. However, the requirement for the type of experimental data depends on the crystal structure of the material and the established modeling setup. The macroscopic uniaxial tensile data is generally used for the calibration of CP model parameters for a material with BCC and FCC crystal structures [25,28]. The author’s previous study showed that at least two stress-strain curves at different orientations are required to properly calibrate the CP model parameters in BCC single crystals [29]. In polycrystalline BCC materials, which capture the heterogeneity/anisotropy of the BCC material, a single stress-strain curve may be sufficient for the calibration. However, the use of a second experimental data point for calibration is highly recommended to avoid non-unique CP parameters. In the present study, due to the lack of availability of much experimental data, a single uniaxial tensile test datum obtained along the extrusion direction (ED) of the steel tube was used for the initial calibration of the CP model parameters using the 3D RVE.

To calibrate the CP model parameters, the virtual tensile test on the representative volume element (RVE) along the ED of the as-received tube was performed. Dream.3D software (Dream 3D-6.5.171) was utilized to generate the 3D RVE model. To that end, a dual phase RVE model consisting of the matrix phase and second phase particles was generated. The grain size, aspect ratio, and crystal orientation for the matrix phase as well as the size distribution of precipitates for second phase particle were obtained from microstructure observations and input to the Dream.3D software for generation of the corresponding RVE. More details on the development of the microstructure-based RVE model is provided in the next section. The generated RVE model was then imported to the Abaqus finite element (FE) commercial package to perform the virtual tensile testing. To perform a virtual tensile test on the RVE along the extrusion direction (ED), the symmetry boundary conditions were applied to three adjacent surfaces, while a tensile displacement boundary condition was defined on the opposite surface along ED, as illustrated in Figure 1. Later, the summation of the reaction force along ED on the surface with a symmetry boundary condition was divided by the initial surface area to obtain the engineering stress, while the average tensile displacement on the opposite surface along ED was divided by the initial RVE size in the same direction to obtain the engineering strain. The synchronized stress-stress data throughout the virtual tensile test would provide the flow behavior of the RVE model. Also, the eight-node linear brick element with reduced integration (C3D8R) was used to mesh the RVE, and a fixed mass scaling factor of 10⁹ is applied to reduce the computational time. To run the simulation, the CP model was incorporated into the commercial Abaqus/Explicit FE code using the user-defined material subroutine VUMAT.

4.2. Microstructure-Based RVE Model

In this section, a 3D RVE model that comprehensively represents the microstructural characteristics of the investigated material is generated to calibrate the dislocation density-based CP model. To that end, Dream.3D software is employed to generate a 3D RVE model consisting of two phases: the matrix phase ferrite, and the second phase particles. For the matrix phase, the distributions of grain size, aspect ratio, and crystal orientation obtained from EBSD measurements are incorporated. The EBSD maps collected from axial and hoop cross-sections of the as-received tube is illustrated in Figure 2. The microstructure on both cross-sections depicts a uniform distribution of fine grains. In addition, these microstructures are comprised of equiaxed grains with some grains slightly elongated toward ED and TD directions, which could be related to the prior production procedure. For this phase, the elastic-plastic deformation condition is considered, in which a total of 12 slip systems of {110} <111> and 12 slip systems of {112} <111> in BCC crystals are assumed to be activated. It is shown that the {110} <111> and {112} <111> slip systems in the BCC crystal equally contribute to the deformation, since the critical resolved shear stresses of the two active slip systems are very similar [30].

As for the second phase particles, the configuration and properties of the precipitates should be correctly defined to effectively model the mechanical behavior at the macroscale. To that end, the second phase particles in the RVE model are generated based on the experimental data, such as average diameter and volume fraction obtained from OM images. As an example, the OM image collected from the hoop cross-section of the as-received tube, along with the corresponding Dream.3D image processing analysis result, is shown in Figure 3. In Figure 3b, the identified precipitants are shown in red, which are unevenly distributed throughout the matrix. As seen, the distribution of second phase particles is random in some regions, while it is aggregated on the grain boundaries in the remaining regions. To simplify the modeling approach, the second phase particles are assumed to be equiaxed and randomly distributed throughout the matrix phase in the RVE. In addition, an isotropic elasticity model with an elastic modulus of 327.58 GPa and Poisson’s ratio of 0.334 was assumed for precipitates [28].

The generated 3D RVE should statistically represent the original material according to three characteristics:

• Grain topology of the matrix phase,
• Orientation distribution of the matrix phase, for which the RVE should contain hundreds of grains to accurately represent the texture,
• Configuration of the second phase particle, for which the RVE’s resolution should be sufficiently fine to account for the small size precipitates.

Figure 4 depicts a 3D RVE model of the as-received steel tube comprised of 670 grains for the matrix phase and 0.5 million second phase particles. For the sake of computational efficiency, this model was found to be the smallest RVE that could accurately represent all the microstructural characteristics of the original material. The size of the generated 3D RVE is determined to be 64 μm × 64 μm × 64 μm, with a total of 32.8 million equal-sized (0.2 × 0.2 × 0.2 μm³) cubic elements. To verify the accuracy of representation, it is crucial to compare the RVE characteristics to the corresponding experimental data. In this regard, the grain size distribution in the matrix phase is compared to the experimental EBSD data in Figure 5, which is found to be accurate. As for the second phase particles, the average equivalent particle diameter and volume fraction obtained from the RVE model are 0.294 μm and 0.028, which is in close agreement with the corresponding experimental values of 0.283 μm and 0.029, respectively. Furthermore, the pole figures of the <001>, <011> and <111> crystal directions for the matrix phase computed from the generated RVE is in good agreement with the measured distributions of the averaged grain orientations obtained from the EBSD data, as depicted in Figure 6. Based on the provided comparisons, it can be confirmed that the generated 3D RVE precisely represents the as-received steel tube in terms of the grain topology and orientation distribution of the matrix phase, as well as the configuration of the second phase particle.

To calibrate the CP model, the calibration parameters are adjusted based on a trial-and-error technique, until the virtual tensile test data performed on the generated 3D RVE matches the experimental uniaxial tensile data. Figure 7 shows the comparison between the engineering flow curves obtained from the CPFEM simulation of virtual tensile test after calibration and the experimental uniaxial tensile test data in ED of the as-received tube. As can be seen, the CPFEM prediction is in reasonable agreement with the experimental data using the material constants listed in Table 1 and the calibration parameters listed in

Table 2. The obtained calibration parameters for dislocation density-based CP model.

Parameter	Value
${\mathit{\tau }}_{\mathbf{0}}$	5 MPa
${\mathit{H}}^{\mathit{\alpha }}$	0.17
${\mathit{\rho }}_{\mathbf{0}}^{\mathit{\alpha }}$	0.1/μm2
${\mathit{k}}_{\mathit{a}}$	186
${\mathit{k}}_{\mathit{b}}$	0.009 μm

. Since the deformation behavior beyond the uniform elongation limit basically depends on the sample geometry, the CPFEM simulation is solely performed in the uniform elongation region, and the post-uniform deformation behavior is neglected in the present study.

As mentioned earlier, the RVE model depicted in Figure 4 should represent the original material in terms of the configuration of second phase particles in the submicron scale as well as the grain topology and texture in the microscale. In the submicron scale, the RVE element size of 0.2 × 0.2 × 0.2 μm³ is sufficiently small to precisely account for the fine size of precipitates with an experimental average particle size of 0.283 μm. On the other hand, in microscale, the RVE model should be large enough to contain a sufficient number of grains to represent the original material in terms of grain topology and orientation distribution. To that end, the size of the generated RVE was determined to be 64 μm × 64 μm × 64 μm. According to the extensive literature, the typical number of elements needed in an RVE model to represent a material is less than one million. However, due to the very small size of the precipitates, the RVE model in Figure 4 has 32.8 million elements to represent the carbon steel tube material. The reason for limiting the number of elements in RVE model is that the computational cost rapidly increases with the increase in the number of elements in FE models. It is worthy to note that the CPU time required for the CPFEM simulation with the RVE model depicted in Figure 4 was 1441-h (60 days) using a 24-core 3.00 GHz processor, while the calibration process required multiple rounds of simulation to optimize the CP model parameters.

4.3. Computationally Efficient Algorithm Based on Simplified RVE Model

In the previous section, a 3D RVE was generated to calibrate the CP model against the experimental uniaxial tensile test data by performing a virtual tensile test. Although the characteristics of the created model was in good agreement with the experimental data in terms of the distribution of precipitates, as well as the grain topology and grain orientation distribution of the matrix phase, it was computationally expensive and inefficient. To improve the efficiency in the developed RVE model without sacrificing the accuracy, a multiscale modeling procedure based on a homogenization scheme is proposed to account for the effect of precipitation strengthening.

The proposed multiscale modeling procedure is based on the superposition of the strengthening effects in the metallic materials. In this approach, the effect of precipitation strengthening on the macroscopic mechanical properties is merged with the material’s intrinsic strength. In other words, an increase in the amount of the material’s intrinsic strength could compensate for the contribution of the precipitation strengthening, through which the explicit presence of particles in the RVE model is no longer necessary. The decoupling of the second phase particles from the RVE model would allow improving the efficiency of the computations by reducing the number of finite elements needed in the RVE model. To precisely account for the effect of precipitates on the mechanical properties through the proposed homogenization scheme, the calibration process is broken down into three sequential steps as described below:

Step 1: homogenization—RVE with precipitates—A microscale RVE model with a high resolution consisting of the matrix phase and second phase particles is generated, as depicted in Figure 8a. This RVE model is determined to be 20 μm × 20 μm × 20 μm with 1 million equal sized cubic elements. The average equivalent diameter and volume fraction of second phase particles are 0.288 μm and 0.029, which is in close agreement with the corresponding experimental values of 0.283 μm and 0.029, respectively. Similar to the RVE model in Figure 4, an isotropic elasticity model with an elastic modulus of 327.58 GPa and Poisson’s ratio of 0.334 was assumed for the second phase particles [28]. In addition, the developed RVE model only contains 30 grains with an average grain diameter of 7.24 μm, which is in good agreement with the experimental value of 7.73 μm. It is worth noting that since the generated RVE model does not contain enough grains, it cannot represent the grain size distribution and orientation distribution of the original material, and thus the Hall-Petch effect is excluded from this calculation. In addition, random orientations are assigned to the grains in the RVE model. The generated RVE model is used for the calibration of the CP model parameters against the experimental uniaxial tensile data in ED of the as-received tube. As can be seen in Figure 9, the computationally obtained flow data (blue curve) is in good agreement with the experimental data. The calibration parameters obtained from this calibration are listed under the Step 1 column in

Table 3. The obtained calibration parameters for dislocation density-based CP model through the homogenization scheme for the effect of precipitation strengthening.

Parameter	Step 1	Step 2	Step 3
${\mathit{\tau }}_{\mathbf{0}}$	120 MPa	120 MPa	5 MPa
${\mathit{H}}^{\mathit{\alpha }}$	0	0	0.19
${\mathit{\rho }}_{\mathbf{0}}^{\mathit{\alpha }}$	0.1/μm2	0.1/μm2	0.1/μm2
${\mathit{k}}_{\mathit{a}}$	338	180	180
${\mathit{k}}_{\mathit{b}}$	0.005 μm	0.009 μm	0.009 μm

.

Step 2: homogenization—RVE without precipitates—The RVE model developed in Step 2, shown in Figure 8b, has similar characteristics and an identical size and resolution as that in Step 1, except that the second phase particles are removed. This model contains 30 grains with an average grain diameter of 7.85 μm, which is in close agreement with the experimental value of 7.73 μm. Similarly, because of the insufficient number of grains, random orientations are assigned to the grains in the RVE and the Hall-Petch effect is excluded from the calculations. To determine the contribution of the precipitation strengthening on the material’s strength, the CPFEM simulation of the virtual tensile test on the RVE model without the precipitates (Figure 8b) is performed using the CP model parameter obtained from Step 1. The obtained flow data from this simulation (green curve) is shown in Figure 9. The difference between the two flow curves, corresponding to the RVE with precipitates and without precipitates, is identified as the contribution of the precipitation strengthening on the material’s strength. It is worth noting that in the CPFEM simulation, these precipitates act as hard barriers to the free motion of dislocations and cause the material to strain harden. Therefore, by removing these precipitates from the CPFEM simulation, the primary effect will be a reduction on the strain hardening portion of the stress-strain curve. As can be seen in Figure 9, the contribution of the precipitation hardening is significant. To compensate for the contribution of the precipitation strengthening, the CP model parameters obtained from Step 1 are recalibrated against experimental tensile test data by performing simulations with the RVE model developed in Step 2 (Figure 8b). Among the CP model parameters, $k_a$ and $k_b$ are the parameters that contribute to the shape of the hardening curve by controlling the speed of dislocation evolution and annihilation. To obtain a perfect agreement between the predicted hardening curve and experimentally measured stress-strain curve, both $k_a$ and $k_b$ are adjusted. The calibrated flow data (red curve) with adjusted calibration parameters (see the Step 2 column of Table 3) are shown in Figure 9. As can be seen, the obtained flow curve is in reasonable agreement with the experimental flow data.

Step 3: homogenization—incorporate the Hall-Petch effect—As mentioned earlier, the Hall-Petch effect was excluded from the simulations during homogenization of the precipitation strengthening due to the insufficient number of grains in the generated microscale RVE models developed in Step 1 and Step 2. In Step 3, the Hall-Petch parameter,$H^{\alpha }$, and the initial friction stress, ${\tau }_0$, are determined by recalibrating the initial yield point. To obtain these two parameters, the virtual tensile test is performed on a 3D RVE without precipitates, but with 2500 grains to accurately represent the grain structure and orientation distribution of the original material. Figure 10a illustrates the generated RVE for Step 3, in which the size of the model is determined to be 105 μm × 105 μm × 105 μm with an element size of 1.75 μm × 1.75 μm × 1.75 μm. The grain size distribution of the generated RVE is compared to the experimental data in Figure 10b, which are found to be very accurate. In addition, the pole figure computed from the generated RVE is presented in Figure 10c, which is also in close agreement with the corresponding experimental data shown in Figure 6b. Performing the virtual tensile test on the generated RVE model using the calibrated CP model leads to the flow curve illustrated in Figure 11. The comparison between the engineering flow curves obtained from the CPFEM simulation and the experimental uniaxial tensile test data shows a close agreement. The Step 3 column in Table 3 shows the calibrated Hall-Petch coefficient along with the rest of the hardening parameters. As can be seen, the three hardening parameters (${\rho }^{\alpha }$, $k_a$, and $k_b$) related to plastic deformation beyond the initial yield were not changed. Based on the observed results, the calibration parameters tabulated in Table 3 corresponding to Step 3 can be adopted for the simulations based on the developed CP model for the as-received material.

The proposed multiscale modeling technique could enhance the CPFEM modeling performance in two ways. Firstly, it could significantly reduce the computation cost at different stages of the modeling procedure, including the generation of an RVE model, the verification of the established model, and conducting simulations to predict the mechanical behavior of the material under complex deformation modes. Based on the comparisons, it is observed that the RVE models depicted in Figure 4 and Figure 10a represent the same material with identical characteristics. The only difference is that the second phase particles are removed from the RVE model in Figure 10a. However, the homogenization scheme is employed to account for the effect of precipitation strengthening on the material’s mechanical properties. In addition, the comparison of results presented in Figure 7 and Figure 11 indicate that carrying out virtual tensile tests with these two RVE models results in similarly accurate predictions of the flow stress of the steel material. However, the CPU time required for completion of the virtual tensile test with the RVE model in Figure 10a takes 12 h with a 24-core 3.00 GHz processor, which is significantly lower than the 1441-h CPU time required for the CPFEM simulation with the RVE model in Figure 4. This result conveys the fact that the proposed multiscale modeling procedure based on a homogenization scheme could successfully improve the computational efficiency of the CPFEM modeling without sacrificing accuracy.

Secondly, the proposed modeling technique is beneficial toward distinguishing the individual effect of precipitation strengthening on each of the CP calibration parameters by breaking down the calibration process into multiple steps. It is well known that the precipitation hardening may impact the elastic, initial yield, and plastic hardening behavior of the material. It was mentioned earlier that the developed CP model has multiple calibration parameters, among which three parameters of ${\tau }_0$, $H^{\alpha }$, and ${\rho }_0^{\alpha }$ define the initial yielding, and two parameters of $k_a$ and $k_b$ define the hardening behavior at the plastic deformation region by evolving the dislocation density of each slip system, ${\rho }^{\alpha }$. Through the proposed multiscale modeling approach, the calibration process is broken into three sequential steps, such that a portion of the calibration parameters is obtained at each step. Using this multiscale technique, the impact of precipitation strengthening on each of the calibration parameters is identified separately, and the mixed effect of second phase particles on the CP model parameters is avoided. This is not possible when using the traditional calibration technique based on the RVE model in Figure 4.

5. Sensitivity of the Macroscopic Mechanical Properties to the Precipitates

As seen earlier, precipitation hardening affects the mechanical properties of metallic materials. The magnitude of the effect on the macroscopic mechanical properties is highly sensitive to the properties of second phase particles at lower scales. It is well known that the microscale properties of precipitates such as size, volume fraction, geometrical distribution, and elasto-plastic properties [5], as well as the submicron scale properties such as particle interface properties, homogeneity of chemical composition, and lattice misfit strain could significantly affect the mechanical properties at larger scales [6,7]. In the present study, the submicron scale parameters are neglected since they are smaller than the resolution of the generated RVE models. However, the influence of microscale precipitate properties on the macroscopic mechanical properties of the material is investigated in this section. To that end, the CPFEM simulation of the uniaxial tensile test is performed on the RVE models with a diverse range of the characteristics of second phase particles to investigate the sensitivity of the macroscopic mechanical properties to the precipitate. The generated RVE models are single grain to eliminate any possible effect caused by grain structure and texture on the mechanical properties. Therefore, the individual effect of the second phase particle on the mechanical properties, as the only varying parameter, is investigated. An example of such RVE models is shown in Figure 12, in which the RVE model is comprised of one grain as the matrix as well as second phase particles with a size and volume fraction of 0.26 μm and 0.025, respectively.

The mechanical behavior of the second phase particles plays a crucial role in determining the contribution of precipitate strengthening on the macroscopic properties. The behavior of the precipitate is basically governed by the alloy system of material. However, it is reported that the size, coherency, and chemical composition of precipitates are the factors that could influence the mechanical behavior of precipitates [31]. The behavior of second phase precipitates are generally assumed to be non-deformable with elastic mechanical properties for simplicity [32], but the deformable behavior with elasto-plastic mechanical properties may better describe the behavior of second phase precipitates under particular conditions [31]. In this regard, the sensitivity of the macroscopic properties to the mechanical behavior of second phase precipitates is investigated through performing a CPFEM simulation of the uniaxial tensile test under various mechanical properties of precipitates.

The effect of the mechanical behavior of second phase precipitates on the macroscopic flow properties is illustrated in Figure 13, in which the effects of deformable precipitates with elasto-plastic behavior and non-deformable precipitates with elastic behavior are shown. Figure 13a depicts the macroscopic flow properties obtained through considering various elasto-plastic properties for deformable second phase precipitates. The elasto-plastic behavior of precipitates is manipulated by increasing the yield strength of the base material at different precipitate to base strength ratios ($P/B={\sigma }_y{}_{Precipitate}/{\sigma }_y{}_{Base}$), while the hardening behavior beyond the initial yield point is considered similar to the base material, as illustrated in Figure 14. Based on the generated elasto-plastic flow data, the CPFEM simulation of the virtual tensile test is performed on the dual phase RVE model with a second phase particle diameter of 0.26 μm and volume fraction of 0.025, as shown in Figure 12. As shown in Figure 13a, the obtained results are compared to the macroscopic properties of the base material, as well as the material with pure elastic precipitates, in which the elastic modulus and Poisson’s ratio of 171.2 GPa and 0.34, respectively, are considered for precipitates. At the precipitate to base strength ratio of 1 (P/B = 1), in which identical mechanical properties are assigned to both base material and second phase precipitates, a slight increase in the macroscopic strength is observed. further increasing the strength of precipitates (higher P/B) leads to an additional increase in the macroscopic strength of the material. The precipitation strengthening ultimately saturates at P/B = 10, in which the precipitates behave like non-deformable precipitates and demonstrate a similar macroscopic behavior to the pure elastic behavior. As observed, the macroscopic properties are quite sensitive to the elasto-plastic behavior of deformable precipitates. The primary reason is that the dislocations cut through and plastically deform the precipitates during deformation of materials with deformable precipitates. Thus, in the materials with deformable precipitates, the amount of precipitation strengthening is basically controlled by the intrinsic mechanical properties of the precipitates.

Figure 13b depicts the macroscopic flow properties obtained through considering various elastic properties for non-deformable second phase precipitates. In carbon steels, various Fe-C compounds with a diverse range of properties and crystal structures can be formed. In a study performed by Liu et al. [32], the properties of the Fe-C compounds generally found in carbon steels have been investigated. The summary of the elastic properties and crystal structure of the proposed Fe-C compounds are tabulated in

Table 4. The structure and elastic properties of the Fe-C compounds generally found in carbon steels [2].

Compound	Crystal Structure	Elastic Modulus (GPa)	Poisson’s Ratio
Fe2C	Orthorhombic	171.2	0.34
Fe3C	Hexagonal	327.58	0.334
Fe5C2	Triclinic	193.51	0.351
Fe5C2	Monoclinic	323.45	0.348
Fe7C3	Hexagonal	227.44	0.336

[32]. These elastic properties are assigned to the second phase particle in the RVE model that is shown in Figure 12 to perform the virtual tensile test. The result of this investigation (Figure 13b) illustrates that the effect of the elastic properties of non-deformable precipitates on the macroscopic mechanical properties is insignificant. This is because the dislocation could not cut through the non-deformable precipitates during the deformation, and could only bypass the precipitates. Therefore, in the materials with non-deformable precipitates, the contribution of precipitates on the macroscopic strength of material is independent of the precipitate properties [31], and is mainly controlled by the topology of precipitates, including size, volume fraction, and configuration of the precipitates across the matrix.

To study the sensitivity of macroscopic properties to the topology of second phase precipitates, the CPFEM simulation of the uniaxial tensile test is performed on RVE models with varying sizes and volume fractions of precipitates. The summary of the RVE models generated for this investigation is provided in

Table 5. The summary of the RVE models generated to study the effect of the size and volume fraction of precipitates on macroscopic properties.

Case #	Eq. Diameter of Particles (µm)	Volume Fraction of Particles
Case 1	N/A	0
Case 2	0.26	0.025
Case 3	0.44	0.025
Case 4	0.61	0.025
Case 5	0.26	0.05
Case 6	0.26	0.1

. Both non-deformable and deformable behavior of precipitates are considered when performing the CPFEM simulations. In the case of non-deformable precipitates, the elastic properties of Fe₂C provided in Table 4 are assigned to the second phase precipitates. On the other hand, elasto-plastic behavior with a precipitate to base strength ratio of 2 (P/B = 2) is considered for the deformable precipitates.

The flow curves obtained from simulations with non-deformable precipitates on the RVE models introduced in Table 5 are illustrated in Figure 15. First and foremost, the presence of second phase particles in the RVE model leads to the precipitation strengthening, regardless of the size and volume fraction of the particles. According to the data presented in Figure 15a, the increase in the size of the precipitates at a constant volume fraction, i.e., a lower number of particles in the RVE model, would weaken the effect of precipitation hardening. The precipitation strengthening is basically caused by the ability of precipitates to resist against the movement of dislocations on the slip systems. In this case, although the few large particles could locally surge the resistance against dislocation movements, resulting in the overall precipitation strengthening, a more uniform distribution of finer particles may lead to a stronger precipitation strengthening effect. Furthermore, Figure 15b depicts the effect of increasing the volume fraction of precipitates at a constant particle size on the macroscopic mechanical properties. As shown, the increase in the volume fraction in the RVE model not only elevates the strength of the material, but also increases the strain hardening effect. Like the effect of particle size, the increase in the volume fraction of precipitates leads to a huge number of particles spread throughout the microstructure. This number of particles significantly enhances the resistance against the movement of dislocations, and creates dislocation forests behind the precipitates. The dislocation forests may lead to a more intense strain hardening behavior, in addition to an increased material’s strength.

The flow curves obtained from simulations with deformable precipitates on the RVE models introduced in Table 5 are illustrated in Figure 16. Like the behavior previously observed for the material with non-deformable precipitates in Figure 15, the presence of second phase particles in the RVE model leads to precipitation strengthening, regardless of the size and volume fraction of the particles. Figure 16a demonstrates the effect of particle size at a constant volume fraction on macroscopic properties, in which it is seen that the magnitude of precipitation strengthening is independent of the particle size. Figure 16b depicts the effect of an increase in the volume fraction of precipitates at a constant particle size on the macroscopic mechanical properties. Like the behavior observed for non-deformable precipitates in Figure 15, the increase in the volume fraction in the RVE model not only elevates the strength of the material, but also increases the strain hardening effect. However, the amount of precipitation strengthening in the case of deformable precipitate behavior is less significant. To quantify this comparison, the sensitivity of the macroscopic strength to the size and volume fraction of precipitates, considering both deformable and non-deformable precipitate behaviors, is plotted in Figure 17. In this figure, the macroscopic strength at 10% strain is plotted against a wide range of precipitate size-to-volume fraction ratios. As seen, two different behaviors for the sensitivity of the macroscopic strength to the size-to-volume fraction ratio are observed for both deformable and non-deformable precipitate behaviors. Below the size-to-volume fraction ratio of 10, the increase in size and reduction in volume fraction of the precipitates leads to a reduction in the magnitude of precipitation strengthening. In this region, the sensitivity of the macroscopic strength to the topology of precipitates in non-deformable precipitates is more significant in comparison to deformable precipitates. As mentioned earlier, since the dislocations could cut through the deformable precipitates in deformable precipitates, these precipitates cannot act as a rigid barrier against the movement of dislocations. Therefore, the amount of precipitation strengthening is controlled by the intrinsic properties of the precipitates rather than the size-to-volume fraction of the precipitates. However, in non-deformable precipitates, the dislocation could only bypass the precipitates. Therefore, the contribution of precipitates on the macroscopic strength of material is independent of the precipitate properties, and is mainly controlled by the topology of precipitates, including the size-to-volume fraction. Beyond the size-to-volume fraction ratio of 10, the precipitation strengthening becomes independent of the topology of precipitates. The main reason is that the presence of precipitates becomes ineffective in creating resistance against the movement of dislocations on the slip systems at extremely large size or low volume fractions of precipitates. In this condition, the magnitude of the precipitation strengthening reaches a plateau, and a minimal impact of precipitates on the macroscopic strength is observed.

As reported in the literature, the position of the second phase particles in the microstructure may affect the macroscopic mechanical properties [1,32]. To assess the magnitude of this effect, various 3D RVE models with a different distribution of precipitants is constructed and examined through virtual tensile testing. To that end, poly-grain RVE models are constructed in a way that a fraction of the second phase particles is placed on the grain boundaries, while the remainder is spread inside the grains. The 3D RVE models with non-deformable second phase particles constructed for this investigation are shown in Figure 18. As can be seen in Figure 18a, all the second phase particles are randomly distributed inside grains, while none are placed on the grain boundaries. On the other hand, Figure 18b displays the RVE model, in which 50% of precipitates are placed on the grain boundaries, and the remainder is evenly distributed inside the grains. The equivalent diameter of particles, volume fraction of particles, and equivalent grain diameter in both RVEs are identical and equal to 0.36 μm, 0.08, and 7.75 μm, respectively, while random grain orientations are assigned in both models. Figure 19 depicts the flow curves obtained from the virtual tensile tests performed on the RVE models presented in Figure 18. As can be seen, placing more second phase particles on the grain boundaries would slightly increase the strength of material, which may be due to the binding of grain boundaries by the precipitates during the deformation [32]. However, the magnitude of strengthening is less than 2% for the increasing of the fraction of second phase particles placed on the grain boundaries from 0% to 50%, which is determined to be insignificant.

6. Conclusions

In the present study, a CPFEM modeling approach was adopted to investigate the effect of the second phase particles on macroscopic mechanical behavior. To that end, a microstructure-based RVE model was generated using the experimental data, including properties and distribution of precipitates in a submicron scale, as well as the grain topology and orientation distribution of grains in a microscale. Because of the multiscale, multi-phase nature of this investigation, a computationally efficient multiscale modeling approach based on a homogenization scheme was proposed to account for the effect of second phase particles on deformation behavior. Lastly, a parametric analysis was performed to investigate the sensitivity of the macroscopic mechanical properties to the configuration and properties of the second phase particles. The main findings of the present study are as follows,

• By comparing the simulation results to corresponding experimental data, it was observed that the proposed multiscale modeling approach successfully accounted for the effect of second phase particles on the deformation behavior. In addition, through the proposed algorithm, the computational cost was reduced by more than 99% for an identical simulation, confirming the vastly improved efficiency of the model.
• It was observed that the topology as well as the mechanical properties of second phase precipitates may heavily affect the macroscopic mechanical properties. It was observed that the amount of precipitation strengthening is controlled by the intrinsic properties of the precipitates in deformable precipitates, while also controlled by topology, including size and volume fraction of precipitates in non-deformable precipitates. In addition, the increase in size and reduction in volume fraction of the precipitates leads to a reduction in the magnitude of precipitation strengthening in both deformable and non-deformable precipitates.
• It was found that the geometrical distribution of second phase particles throughout the microstructure might have limited impact on the macroscopic mechanical behavior of the material.