Multiscale modelling of precipitation hardening: a review

Precipitation hardening, a cornerstone of alloy strengthening, finds widespread application in engineering materials. Comprehending the underlying mechanisms and formulating models bear crucial significance for engineering applications. While classical macroscopic theoretical models based on the line tension model have historically guided research efforts, their reliance on simplifications, assumptions, and parameter adjustments limits their predictability and expansibility. Moreover, the challenge of understanding the intricate coupling effects among various hardening mechanisms persists. One fundamental question to achieve the transition of material design paradigms from empirical trial-and-error methods to predictive-and-design approaches is to develop more physics-based multiscale modelling methods. This review aims to elucidate the physical mechanisms governing precipitation hardening and establish a tailored bottom-up multiscale modelling framework to steer the design of new alloys. The physical scenarios of precipitation hardening are firstly summarized, including particle shearing, Orowan bypass, and dislocation cross-slip and climb. Afterwards, an in-depth discussion is given regarding the application of macroscopic models and their correlation with the mechanisms and precipitation characteristics. As for the multiscale modelling methods, we categorize them into three main types: slip resistance based approaches, misfit stress field based approaches, and energy based approaches. By integrating multiscale modelling with the physical scenarios, we systematically addressed the key idea of the multiscale coupling framework, and their scale transfer procedure, applicability, advantages, and limitations. Some examples of coupling different types of multiscale methods and considering precipitates with complicated shapes are also presented. This study not only furnishes insightful comprehension of precipitation hardening, but also guides the development of multiscale modelling methodologies for other types of hardening effects in alloys.


Introduction
Precipitation hardening is one of the most widely used methods to manipulate the properties of metallic alloys, and especially to reach high strength.The interfacial nature (coherent/incoherent), shape, size, and spatial distribution of precipitates all play critical roles in precipitation hardening, collectively affecting the final properties of materials.
Many methods have been proposed to tune the precipitates, for example, through complex solid-state phase transformations induced by thermal and/or thermo-mechanical treatments (Sun et al. 2015;He et al. 2016;Tang et al. 2019;Gwalani et al. 2019;Jang et al. 2021;Jiang et al. 2023).
Pioneers have proposed several classical macroscopic theoretical models aimed at describing the precipitation hardening effects, which are primarily based on the line tension model.Most of these analyses focused on the motion of dislocations gliding in a single slip plane, where dislocations bypass particles by either Orowan bypass or particle cutting (Brown and Ham 1971;Argon 2007).In addition, precipitates can also be overcome through cross-slip and climb mechanisms, which has been used to successfully explain experimental findings, such as the formation of prismatic loops (Hirsch 1957;Humphreys and Hirsch 1970).Undoubtedly, these models have significantly contributed to our understanding of precipitation hardening behaviors over an extended period.Numerous detailed reviews on precipitation hardening can be found in (Kelly 1963;Brown and Ham 1971;Reppich 1982;Ardell 1985;Martin 1998;Courtney 2005;Argon 2007;Cui and Cui 2024).However, classical macroscopic models generally necessitate the adjustment of coefficients to align with experimental findings, relying on multiple simplifications and assumptions.For example, the simplification of the precipitate shape to a spherical form, the presumption of uniform distribution of the precipitate within the matrix, and sometimes neglecting the local interaction between precipitates and dislocations (Ham 1968;Brown and Ham 1971;Bacon et al. 1973;Werner Hüther and Reppich 1978).Additionally, these models are basically limited to analyzing each hardening mechanism independently (Jiang et al. 2017;Wang et al. 2018;Zhang et al. 2019;Ma et al. 2019), which makes it challenging to account for the coupling effect of several mechanisms in concert.Moreover, the development of new high-strength metallic alloys involves complex physical mechanisms that are influenced by numerous factors (Alaneme and Okotete 2019).Therefore, the application of macroscopic models for guiding the design and enhancement of novel high-strength alloys encounters inherent limitations.
In light of these constraints, there is an urgent need to establish more advanced methods for understanding precipitation hardening and further guiding new ultra-strong alloy design.Relying solely on macroscopic models or experiments is no longer sufficient.It is crucial to elucidate the underlying mechanisms of precipitation hardening crossing the macro, micro and nanoscales.This necessitates a thorough comprehension and quantification of the influences of different mechanisms on precipitation hardening.A bottom-up, physically-based theoretical framework for precipitation hardening is expected to significantly enhance design efficiency and cost savings.Furthermore, multiscale models for precipitation hardening could combine the advantages from different scale and link microstructure to properties, signifying the transition from the trial-anderror era to the prediction-and-controlling era.
Many previous studies have been carried out along this line (Shin et al. 2003;Monnet 2006;Takahashi and Ghoniem 2008;Ji et al. 2014;Liu et al. 2017Liu et al. , 2023;;Kim et al. 2017;Santos-Güemes et al. 2020;Chatterjee et al. 2021;Zhang and Sills 2023).The key idea is mainly to develop the bottom-up precipitation hardening model by analyzing the detailed interaction mechanisms between dislocations and the precipitates, such as based on molecular dynamics methods, dislocation dynamics methods.However, to the best of our knowledge, there remains a gap in systematically summarizing efficient approaches within multiscale computational modeling framework for depicting precipitation hardening.In the existing literatures, the effect of precipitates is introduced in multiscale models through multiple ways, such as by introducing the resistance stress, the stress field of the precipitates, or the stacking fault energy.How can these methods benefit most from multiscale modelling?Which kind of underlying mechanisms can be considered for different methods?What are the advantages and limitations of different methods?How to choose the tailored solutions to address specific problem of precipitation hardening?
This review aims at answering these questions and providing a comprehensive framework for employing multiscale modelling approaches to facilitate the design of advanced alloys.In "Physical scenarios of precipitation hardening" section, we commence by delineating the physical diagram of precipitation hardening, as well as essential parameters used for characterizing precipitates including their shapes and eigenstrain.In "Macroscopic models of precipitation hardening" section, the macroscopic models of precipitation hardening are discussed.In "Multiscale model based on slip resistance" -"Multiscale model based on energy" sections, different types of multiscale modelling methods, including the slip resistance based approaches, the misfit stress field based approaches, and energy based approaches, are presented, respectively.Further, "Application" section delves into some examples of applications of these multiscale methods.Finally, "Conclusions" section summarizes this review.
It is worth noting that the multiscale modelling method can provide a universal framework for understanding the effects of voids, bubbles, inclusions, etc., and their interactions with dislocations, although this study focuses on precipitation hardening.

Theory of precipitation strengthening
The objective of this section is to clarify the physical scenarios of precipitation hardening briefly, as schematically shown in Fig. 1.The additional stress required for dislocations to bypass precipitates depends on the underlying dislocation escape mechanism.The mechanism for dislocation planar motion includes particle shearing and Orowan bypass mechanisms, while cross-slip and climb mechanisms correspond to non-planar motion of dislocations.The dominated mechanism depends on various microstructural parameters, such as alloy composition, lattice misfit, precipitate size, morphology, volume fraction, as well as external conditions like applied stress and the operating temperature.Here, we mainly focus on the planar-motion mechanisms, and simply discuss the dislocation cross-slip and climb mechanism.
As shown in Fig. 1, the dislocation is obstructed by the precipitate and then escapes through a combination of mechanisms: (a) particle shearing (Hirsch 1957;Yashiro et al. 2006), where the dislocations cut through the particles; (b) Orowan bypass (Orowan 1948), where the dislocations bow out between particles under the applied stress, leaving dislocation loops around the particles; (c) cross-slip for screw dislocations (Hirsch 1957;Ashby and Johnson 1969;Humphreys and Hirsch 1970;Xiang et al. 2004), for example by surmounting a particle via forming prismatic loops as shown in Fig. 1; (d) climb for edge dislocations (McLean 1985), which involves the diffusion of point defects.
Whether the precipitate and the matrix satisfy the coherency relationship often affects the dislocation-precipitate interaction, which in turn affects the dominant precipitation hardening mechanisms, i.e., particle shearing or Orowan bypass mechanisms, as illustrated in Fig. 1.For the coherent precipitate, both particle shearing and Orowan bypassing are possible, which are independent and compete with each other.Dislocations encountering precipitates tend to select the least obstructed pathway with the lowest critical shear stress for shearing or bypassing these obstacles, which corresponds to the lowest required external energy in a macroscopic scale.The size of the coherent precipitate is a key factor determining the dominant mechanism.The dislocation can easily shear the small-sized precipitates, while the Orowan bypass mechanism dominates when the coherent precipitate exceeds a critical size.For the incoherent precipitate, dislocations usually bypass the particle, leaving behind a dislocation loop around the particle rather than penetrating it, as incoherent interface has no continuity of lattice planes.
Specifically, the long-range and short-range interaction (in Fig. 1) during dislocation approaching the precipitate corresponds to different mechanisms, which influence the hardening effects collectively.The main mechanisms in particle shearing are coherency strengthening �σ CS , modulus strengthening �σ MS and order strengthening �σ OS .Of these, the first two ( �σ CS and �σ MS ) operate in the long range before the possible occur- rence of shearing, while the order strengthening occurs during shearing in the short range.As shown in the process of dislocation approaching precipitate in Fig. 1, the dislocation is obstructed not only by the strain field generated by precipitate owing to lattice misfit and modulus mismatch, but also by the formation of higher energy bond along the slip plane inside the precipitate when shearing, i.e., antiphase boundary (APB) (Milligan and Antolovich 1991;Long et al. 2017).Therefore, throughout the whole physical process of particle shearing (from stage (1) to (3) of particle shearing in Fig. 1), the larger value of �σ CS + �σ MS and �σ OS dominates the degree of hardening in particle shearing (Wang et al. 2018).
Similar to the particle shearing, the Orowan bypass mechanism also involves the longrange and short-range interactions.The different origins of the misfit strain for a coherent precipitate and an incoherent one are demonstrated in Fig. 2, where the effects of modulus mismatch is not considered.For the coherent precipitate, both lattice misfit and modulus mismatch could affect the misfit strain field.As shown in Fig. 2a-b, the interface between the coherent precipitates and the matrix is continuous and their lattice constants are close, so the lattice misfit δ can be calculated by where a p is the lattice constant of precipitates and a is the lattice constant of the matrix.While incoherent precipitates do not have the concept "lattice misfit" as their interface with the matrix is not continuous, shown in Fig. 2c-d.Instead, both volume misfit � = �V /V and modulus mismatch affect the misfit strain field for incoherent precipi- tates, contributing to the long-range interaction, where V is the volume of the uncon- strained hole in the matrix and V − V is the volume of the unstrained precipitate (Porter 2009).Note that the interaction mechanisms shown in Fig. 1 between dislocations and coherent/incoherent precipitates can change under certain conditions, so the escape mechanisms for dislocation gliding on the plane when obstructed by precipitates are not rigidly fixed.On the one hand, incoherent precipitates might be sheared by dislocations under larger deformation (Cepeda-Jiménez et al. 2019).For instance, the incoherent precipitate (θ phase) in Al-Cu alloys, cannot be sheared by dislocations in most cases, but θ phase might be sheared under large strains due to localized stress concentrations at precipitate interfaces caused by dislocation accumulation (Honeycombe 1984).On the other hand, the precipitate may undergo local phase transformation due to interaction with the dislocation, which may result in the precipitate partially sheared and the rest bypassed by dislocation (Bonny et al. 2019).What's more, the interaction between dislocations and precipitates with mixed-type interfaces (some combination of coherent, incoherent, or semi-coherent interfaces), such as θ′ precipitates in the Al-Cu alloy with both coherent and semi-coherent interfaces, is more complicated (Porter 2009).These precipitates with mixed type interfaces prefer to grow as disc-like or plate-like morphologies, where the broad face is typically coherent but the rest is not, leading to different interaction mechanisms (Awan and Khan 2017).Thus, when analyzing the precipitation hardening due to this kind of precipitates, careful consideration and appropriate simplification should be adopted.
Another type of precipitate, i.e., the semi-coherent precipitate, is not shown in Fig. 1 to avoid complicating the diagram and confusing the reader.When the lattice misfit calculated by Eq. ( 1) is relatively large, the difference between the precipitate and the matrix cannot be accommodated by simply stretching.Instead, microstructural defects such as misfit dislocations may form across the semi-coherent interface partially to relax the misfit.The origin of the misfit strain of semi-coherent precipitates is basically identical to that of the coherent precipitates, but semi-coherent precipitates usually cannot be sheared by dislocations.This main difference between coherent and semi-coherent precipitates results in an adjustment in the calculation of the misfit strain (Zbib et al. 2011).Taking spherical precipitate as an example, we consider the dilatational misfit strain caused solely by lattice misfit (and misfit dislocations).For coherent precipitates, the eigenstrain (or transform strain) ε T caused by lattice misfit is given as where I is the identity tensor.The eigenstrain for a semi-coherent spherical particle can be approximated by (Matthews and Blakeslee 1974) where b is the Burgers vector of the misfit dislocation, and is the misfit dislocation spacing. (2) As for the non-planar motion of dislocations, it is demonstrated that cross slip depends critically on both the local stress state (Kang et al. 2014) and the dislocation line length (Xu et al. 2017).Dislocation climb plays a significant role in the dislocation creep behavior of particle-containing materials, where the climb velocity is related to the climb force perpendicular to the original slip plane (Xiang and Srolovitz 2006;Liu et al. 2021).

Basic information about the precipitates
In describing precipitates, fundamental physical parameters such as shear modulus, lattice structure, and lattice constants play a crucial role.These parameters can be utilized to determine whether precipitates are coherent with the matrix, whether lattice misfit or modulus mismatch exists, and so on.The lattice constant and structure of precipitates can not only be measured by experimental techniques, such as transmission electron microscopy (TEM) or X-ray diffraction (XRD), but also be calculated by the ab initio methods, i.e., Density functional theory (DFT).
Furthermore, it is essential to describe the external physical features of precipitates, such as their shape, size, and spatial distribution.The statistics commonly used to describe the spherical precipitates include: the number of spherical particles per unit volume N v , the number of particles per unit area in the plane N s , the volume fraction f, the average particle radius R, the average radius at the intersecting cross-section between the dislocation slip surface and the particles r, the average area of the intersecting cross-section π < r i 2 > (the subscript i denoting the ith cross-section of the particles), and the inter-and intra-polygonal precipitate distances L and L′ respectively.The relationships between the above statistics are schematically illustrated and discussed in Fig. 3.These external features can be determined not only through characterization of microscopic experiments but also through calculations employing multiscale methods.For example, phase-field approach is sometimes used to simulate evolution of precipitations during heat treatment process or phase transformation process (Ji et al. 2014;Kim et al. 2017;Liu et al. 2017;Liu and Nie 2017).In this kind of model, energy parameters controlling the evolution of precipitates, including chemical free Fig. 3 Information used for statistically describing spherical precipitates (Ardell 1985;Chatterjee et al. 2021;Martin 1998).a The average particle radius is R.The average radius intersected by a random slip plane is r = πR/4.The volume fraction is defined as f = 4/3πR 3 N v , where N v is the number of spherical particles per unit volume.Since the number of particles per unit area intersected by a plane is is the radius of the particle on ith planar section, the average < r i 2 > = 2/3R 2 allows us to conveniently measure the volume fraction as the area fraction f = π < r i 2 > N s on any planar section of the material.By means of the identity < r i 2 > = 32 < r i > 2 /(3π 2 ), the volume fraction can also be measured as f = 32r 2 N s /(3π), where we have let r = < r i > .b The particle spacing L on the slip plane, and the inter-particle spacing L′ = L-2r energy, interfacial energy, lattice parameters, and elastic constants, can be obtained through first-principles density functional theory and molecular statics simulations.
Factors influencing the shape of precipitates include interfacial energy E i and elastic energy E e .While spheres or ellipsoids are common simplifications of precipitate shapes, needle-like or plate-like forms are also widely observed.Under the condition of constant precipitate volume, the equilibrium shape corresponds to the minimum values of E i + E e , where E i and E e are proportional to the particle's surface area and volume, respectively (Johnson and Voorhees 1992;Khachaturyan et al. 1988).What's more, as the particle size increases, with a decrease in specific surface area, the elastic energy, proportional to the volume of the particle, becomes the dominant factor that affects the precipitate shape (Wang et al. 2018).
From a quantitative perspective, we can estimate the influence of elastic energy and interfacial energy on the shape of precipitates.For example, in the cubic matrix, the relative importance of the elastic energy and the interfacial energy can be evaluated through the characteristic parameter L , which are defined as (Thompson et al. 1994) where, ǫ represents the strain.ǫ is equal to the lattice misfit δ for the case of dilatational misfit.C 44 is one of the three independent elastic constants of the cubic matrix.l is the characteristic size of particles.s is the average specific interfacial energy.Figure 4 illus- trates the evolution of particle morphology by plotting a R 2 as a function of L , where a R 2 , one of the real parts a R n of Fourier coefficients used in defining a particle shape, is used to classify each family of particle morphology (Thompson et al. 1994).In Fig. 4, L=0, where elastic energy is zero, corresponds to spherical shapes.As L increases, it smoothly (4) L = ǫ 2 C 44 l/s Fig. 4 Morphology evolution of particles with purely dilatational misfit changes with the characteristic parameter L , where the vertical axis is the Fourier coefficient a R 2 to represent different particle shapes.It shows that the bifurcation from the four-fold symmetric cuboid to the two-fold symmetric shapes (plate or needle) occurs at the critical point 'O' with value ( L= 5.6) (Thompson et al. 1994) transitions into cuboidal shapes.The four-fold symmetric shape bifurcates to two-fold symmetric shapes (i.e., needle-like or plate-like) along specific directions at a critical value of L (as marked by point 'O' in Fig. 4).In addition, the presence of such a bifurca- tion implies that the four-fold family of shapes is unstable above the bifurcation point and that two-fold symmetric shapes are the energy-minimizing one (Thompson et al. 1994).As L continues to increase, the particle becomes more plate-like indicating the increasing importance of elastic energy as L increase (Su and Voorhees 1996).
Moreover, lattice misfit plays a decisive role in determining precipitate morphology, especially when particle sizes are comparable.This is because L depends on not only the particle radius R but also the lattice misfit δ , as shown in Eq. ( 4) (Wang et al. 2018), since ǫ is equal to the lattice misfit δ for the case of dilatational misfit strain (Eq. ( 2)).Experimental results confirm a noteworthy trend: an increase in the L value triggers a transition from spherical/ellipsoidal shapes ( ǫ < 0.2%) to cuboidal B2 nanoprecipitates at moderate ǫ ( ǫ ∼ 0.4%), and eventually to weave-like microstructure at large ǫ ( ǫ > 0.6%) (Ma et al. 2018).It's worth noting that the ellipsoidal shapes, as mentioned above, are not depicted in Fig. 4 and require further studies.From an energy equilibrium perspective, the ellipsoidal shape typically occurs at a relatively smaller L .However, the eigenstrain in this case is not dilatational with the same principal values, resulting in the longer axis of ellipsoid oriented perpendicular to the direction of the larger principal eigenstrain (Jou et al. 1997).

Macroscopic models of precipitation hardening
Most of the previous work focuses on developing the precipitation hardening models corresponding to the shearing and Orowan mechanisms mentioned above.Here, we briefly summarize the corresponding commonly used models, as listed in Table 1 (Jansson and Melander 1978;Ardell 1985;Wen et al. 2013;Wang et al. 2018).These equations are essentially based on the theory of line tension and experimental calibration.
To gain a deeper understanding of the transformation of the dominant precipitation hardening mechanism and the primary selection of the most influential mechanism in engineering practice, Fig. 5 illustrates the relationship between precipitation hardening mechanisms and the size of precipitates, by taking spherical coherent precipitates as examples.Here, the volume fraction (f) keeps as a constant.According to Table 1, the variations in yield strength with changing precipitate radius across distinct hardening mechanisms satisfies where, R is the average size of precipitate.�σ OS ∼ R 0 means that order strengthening remains constant despite an increase in R, as depicted by the horizontal line in Fig. 5. Modulus strengthening is generally weaker compared to other mechanisms (Nembach and Neite 1985;Santos-Güemes et al. 2018;Takahashi and Ghoniem 2008), so we approximate it as �σ CS + �σ MS ∼ √ R resulting from �σ CS ∼ √ R and �σ MS ∼ R 0.275 .The dominant mechanism in particle shearing is the one with the larger value between (5) �σ CS + �σ MS and �σ OS in the whole process when dislocations shear throughout the precipitates.As shown in the blue region of Fig. 5a, when the size of precipitates is small ( R < R c1 ), order strengthening dominates, while the lattice strengthening and modulus Table 1 Precipitation hardening models, which are used to calculate the strength increments by both particle shearing mechanism and Orowan bypass mechanism.M = 2.73 for body-centeredcubic structure and M = 3.06 for face-centered-cubic structure (Hosford 2005).α ε = 2.6 (Ardell 1985;Jansson and Melander 1978).(Argon 2007;Jansson and Melander 1978).The power exponent m = 0.85 (Ardell 1985) Mechanisms Equations Parameter explanation (See the table caption for parameter values.)strengthening dominates when the precipitate radius exceeds R c1 .We can understand this phenomenon by considering the area-to-volume ratio as follows.The interaction between dislocations and precipitates in order strengthening is positively correlated with the cross-sectional area in precipitate sheared by dislocation, while coherency strengthening plus modulus strengthening is positively correlated with the volume of precipitates.Therefore, when the precipitate size is relatively small ( R < R c1 ), the ratio of cross-sectional area to volume is larger, making order strengthening effects more pronounced.As the radius of coherent precipitates continues to increase until it exceeds the critical radius R c2 , where �σ Orowan is less than �σ CS + �σ MS , dislocations tend to shear through precipitates in a less obstructed manner, and the dominant mechanism shifts from dislocation shearing to Orowan bypass.Therefore, for coherent precipitates under constant volume fraction conditions, there exists an optimal size at which the increase in yield strength is maximized.This critical threshold results from a balance between Orowan stress (decreasing with precipitate size) and the repulsion caused by modulus and lattice misfit (increasing with precipitate size) (Marquis and Dunand 2002).
In most alloys, the mechanism phase diagram of hardening mechanisms associated with coherent precipitates can follow the relationship illustrated in Fig. 5a.As the size of precipitates increases, the dominant mechanism shifts from �σ OS to �σ CS + �σ MS , and then to �σ Orowan (Jiang et al. 2017;Wang et al. 2018;Zhang et al. 2019;Ma et al. 2019).However, it's also possible that the relationship governing the dominant mechanism as precipitate radius increases may change, due to the unique characteristics of precipitates.For instance, if the antiphase boundary energy γ APB is relatively high in order strengthening, the mechanism phase diagram of hardening mechanisms could shift from (a) to (b) in Fig. 5.In this scenario, within the dislocation shearing mechanism (blue region), �σ OS always surpasses �σ CS + �σ MS , indicating a transition from order strengthening dominated to Orowan bypass hardening dominated, with a corresponding critical radius of R c .Considering the extreme case in Fig. 5b, when γ APB of coherent pre- cipitates tends toward infinity, it can be assumed that these coherent precipitates are not penetrable.For a special precipitate case, the dominant hardening mechanism no longer depends on the size of precipitates, as dislocations exclusively bypass the precipitates by Orowan looping mechanism.Nonetheless, the fundamental principle for determining the dominant mechanism remains unaltered.Curves depicting how precipitation hardening mechanisms change with radius can be obtained based on the characteristic parameters of the actual precipitates listed in Table 1, thus enabling specific analyses.
Most of the above models assume a uniform distribution of precipitates.To illustrate the effect of non-periodic distribution, researchers carried pioneering studies (Friedel 1964;Kocks 1966;Foreman and Makin 1966).The Bacon-Kocks-Scattergood (BKS) model (Bacon et al. 1973) have shown that the critical shear stress required for dislocation passage through non-periodic obstacles is lower than that required for passage through periodic obstacles.
For particle strengthening at high temperatures, dislocations could climb over particles with point defect diffusion, resulting in a decrease in the critical bypass stress.In the 'local climb' model (Brown and Ham 1971;Shewfelt and Brown 1977) and the 'general climb' model (Lagneborg 1973;Hausselt and Nix 1977;Rösler and Arzt 1988), the critical bypass stress is determined by the competition between dislocation climb and the increase of dislocation line length.

Multiscale model based on slip resistance
The macroscopic models for precipitation hardening presented in "Macroscopic models of precipitation hardening" section provide an incomplete description of the precipitatedislocation local interaction.In addition, these models only analyze single hardening mechanisms, making it challenging to account for the impact of multiple mechanisms on hardening.This limits the applicability of the existing models to some extent.
On the other hand, the multiscale modelling methods are developed to overcome the difficulties, which can be classified into: slip resistance based approaches, the misfit stress field based approaches, and energy based approaches."Multiscale model based on slip resistance"-"Multiscale model based on energy" sections will elaborate these multiscale modelling methods in detail.
Many researchers pay key attention to the obstruction for dislocations only inside the area occupied by the precipitates, which corresponds to stage (2) of particle shearing or Orowan looping in Fig. 1.Introducing slip resistance to the movement of dislocation within the precipitate area is a simple and straightforward method for considering the hardening effects.This can be achieved by either applying a critical shear stress τ c when dislocations enter the precipitate, or adopting different drag coefficients in the matrix and the precipitate, as schematically by case 1 or case 2 in Fig. 6 (Bocchini and Dunand 2018;Chen et al. 2022a, b;Monnet 2006;Monnet et al. 2011;Santos-Güemes et al. 2022;Santos-Güemes and LLorca 2023).When the applied resolved shear stress on the dislocations is greater than τ c , the dislocations could penetrate the precipitate, and the direc- tion of the Peach-Koehler force induced by τ c keeps opposite to the dislocation motion.Otherwise, if the dislocation under the applied stress cannot overcome the obstruction, the dislocation will bypass the precipitate, with an Orowan loop left behind around the precipitate, as shown in Fig. 1.
The critical shear stress, required for dislocations to shear the precipitates, physically should correspond to an increase in the yield strength in macroscopic scale, depicting the effect of order strengthening.It tends to ignore the obstruction of dislocation motion before shearing into the precipitate, which corresponds to coherent strengthening and modulus strengthening in macroscopic scale.Of course, the researchers can approximately consider all these strengthening effects, by setting the value of τ c but neglecting the detailed long-range interaction.For incoherent precipitates, the critical shear stress or drag coefficients are typically set infinitely large, keeping the precipitates not penetrable by dislocations, regardless of the magnitude of the applied stress (Monnet 2006;Santos-Güemes et al. 2018;Santos-Güemes and LLorca 2023).When it comes to penetrable precipitates, the critical shear stress or drag coefficients are finite and should be treated with more care (Bocchini and Dunand 2018;Santos-Güemes et al. 2021;Chen et al. 2022a, b).
Here, we summarize slip resistance based methods and categorize them into three cases based on the scale of the simulation methods used (Fig. 6).The first case involves Density Functional Theory (DFT) in combination with Dislocation Dynamics (DD).The second case involves Molecular Dynamics (MD) and DD.The third case involves Crystal Plasticity-Finite Element Method (CP-FEM) in combination with DD or Macroscopic models.Each case will be discussed in detail below.

Case 1: DFT + DD
The shearing of ordered precipitates by the glide of superdislocation in the γ phase, accompanied with the formation of APB, is one of the dominant deformation mechanisms of Ni-base and Co-base superalloys in the low/intermediate temperature regime (Chatterjee et al. 2021;Eggeler et al. 2021).When investigating the hardening behaviors induced by these precipitates, e.g., L1 2 structure precipitate in nickelbased superalloys, order strengthening typically dominates, while effects of modulus mismatch and lattice misfit are weaker (Chatterjee et al. 2021;Chen et al. 2022a, b).Therefore, the APB energy is typically used to express the barrier of penetration into the precipitates for dislocations.When one dislocation segment meets one precipitate, the critical shear stress (CRS) required to penetrate the precipitate can be calculated as (Nembach and Neite 1985;Mohles 2001;Hussein et al. 2017;Bocchini and Dunand 2018) where, b is the Burgers vector, γ APB is the antiphase boundary energy due to disturbed local chemical order of precipitates, which can be calculated from DFT or tested by experiments (Yashiro et al. 2006;Chen et al. 2022a, b).The critical shear stress given by ( 6) c is critical shear stress on slip system α .B is drag coefficient in dislocation mobility law Eq.( 6) can be embedded into the DD simulation framework by defining the stress necessary for one dislocation segment penetration into one precipitate.
Case 2: MD + DD MD could provide the mobility law and the local interaction rule to DD, in which the movement of dislocations is straightly obstructed by the precipitates.In DD, the dislocation mobility can be defined as (Marian et al. 2020): where the dislocation velocity of the elementary dislocation segment is defined as a function of local quantities such as local stress field σ , line orientation ξ , and possibly other parameters like Burges vector b , slip plane normal n and temperature T .For the simple case, when the dislocation motion is phonon drag dominated, the dislocation mobility law is generally described by a linear function as follows where the drag coefficent B is a function of temperature.The researchers distinguish the precipitate with respect to the matrix according to the dislocation mobility law (Eq.( 8)), either by changing the drag coefficient B or considering the resistance stress induced by the precipitate in the calculation of τ .Namely, τ = τ eff − τ c , where τ eff is the resolved shear stress applied on the dislocation induced by the external loading and other defects except for the precipitations (such as other dislocations), τ c is the CRS necessary to shear the precipitate.
On the one hand, tuning the drag coefficient B within the precipitates and the matrix can illustrate the bigger resistance of precipitates to dislocations in comparison to the matrix.The drag coefficient could be calculated by carrying out MD simulations of straight dislocation segments with different characters (Santos-Güemes et al. 2018).On the other hand, the critical shear stress in Eq. ( 8), characterizing the obstruction of precipitates to dislocation, can also be determined by MD, influencing the stress in dislocation mobility law.For instance, it is easy to compute the CRS from the MD simulation using the equation τ c = L ′ (τ max − τ f )/2R , where L ′ stands the internal distance between precipitates, R is the radius of the precipitate, τ max is the maximum applied stress cor- responding to unpinning recorded and τ f is the friction stress for dislocation motion in the matrix (Monnet 2015).The examination of stress-strain history in MD simulation of dislocation-precipitate interaction contributes to defining the critical shear stress.Just as the representative example (Fig. 7) shows, the snapshot (c) (dislocation unpinning) and snapshot (d) in the stress-strain curve correspond to the maximum applied stress τ max and the friction stress τ f , respectively.

Case3: CP-FEM + DD / Macroscopic models
The computational framework of crystal plasticity-finite element method (CP-FEM) is illustrated in Fig. 8. Crystal plasticity model has the capability to simulate the impact of precipitates on deformation at a relatively large scale.Some important parameters within this framework might be acquired from macroscopic models, such as listed in Table 1 "Macroscopic models of precipitation hardening" section, or other simulation methods at lower scales, such as dislocation dynamics (DD).
Crystal plasticity model shows that the deformation gradient F is commonly decom- posed into two parts, where F E represents elastic lattice distortion and F P denotes the crystallographic slip.Then the velocity gradient L can be similarly be decomposed: where, L E and L P are the elastic and plastic velocity gradient, respectively.
For dislocation dominated plasticity, the plastic velocity gradient in the CP-FEM simulation is given by ( 9) (α) ⊗s (α) Fig. 7 Representative example of stress-strain curve for dislocation-precipitate interaction and the corresponding snapshots in MD simulations (Bonny et al. 2019) where the superscripts α denote parameters associated with the αth slip system, includ- ing the slip plane normal n (α) , the slip direction s (α) , and its plastic slip rates γ (α) .The plastic slip rates γ (α) are typically defined by resolved shear stress applied on the slip system, the slip resistance and so on.For example, the plastic slip rate can be simplified from the Arrhenius type (Estrin and Mecking 1984) to its power expression as where n is the rate sensitivity of slip and γ0 is a reference shear rate, and τ (α) and g (α)  c are the resolved shear stress and the slip system strength for dislocation glide on slip system α , respectively.The effects of precipitates on deformation are reflected in the slip sys- tem strength g (α)  c , by adding the CRS necessary for penetrating the precipitates to g (α) c as additional slip resistance (Han et al. 2004;Anjabin et al. 2014;Monnet and Mai 2019;Li et al. 2020Li et al. , 2022)).For example, the slip resistance g (α) c , composed of several contributions can be expressed as follows: where, τ s is the solid solution, τ forest is the forest hardening due to dislocation network, τ HP the Hall-Petch effect related to the grain size, and τ pre is slip resistance induced by precipitates.The value of τ pre could be obtained from macroscopic models, such as listed in Table 1 in "Macroscopic models of precipitation hardening" section, or other simulation methods at lower scales, as elaborated in "Case 1: DFT + DD" and "Case 2: MD + DD" section.( 12) Fig. 8 The Framework of crystal plasticity-finite element method.F,F E and F P are total deformation gradient, elastic deformation gradient and plastic deformation gradient, respectively.L P is the plastic velocity gradient.
The indices α refer to the slip systems, defined by the slip plane normal n (α) and the slip direction s (α) , and γ (α) are the plastic slip rates on each slip system α The additional slip resistance induced by precipitates can be obtained from macroscopic models and then calibrated by experiments (Han et al. 2004;Anjabin et al. 2014;Li et al. 2020Li et al. , 2022)), which are analogous to the equations presented in Table 1.Although the equations that describe the strength of precipitates may differ, they are generally based on line tension theory.In addition, new precipitation hardening rules can be obtained from lower scales.Utilizing DD simulations could refine or extend the macroscopic models (Monnet 2015;Queyreau et al. 2010;Santos-Güemes et al. 2022), which could serve as inputs for CRS to CP-FEM framework (Monnet and Mai 2019).
As depicted in Fig. 8, the plastic deformation gradient F P can be updated after the acquisition of the plastic velocity gradient L P .Subsequently, the new deformation gradient is determined, serving as the starting point for the next cycle of calculations.

Multiscale model based on the misfit stress field
The long-range interaction is mainly induced by the interaction of dislocations and the long-range misfit stress field due to precipitates, as illustrated in Fig. 1.Incorporating the misfit stress field into dislocation dynamics (DD) framework could well simulate the dislocation-precipitate interaction before approaching the precipitates.Cross-slip and climb can be also be better considered in this framework.These methods are illustrated in Fig. 9. "Eshelby's inclusion theory without modulus mismatch" and "Eshelby's inclusion theory with modulus mismatch" section will describe how to calculate the misfit stress field of precipitates based on Eshelby's inclusion theory, neglecting the modulus mismatch effects or not."Numerical methods for calculating misfit stress field" section will describe how to obtain the misfit stress field based on numerical methods, which can provide solutions to the more complicate cases, such as complex shapes of precipitates and with modulus mismatch.
It should be emphasized that both the theoretical and numerical methods, used to solve the stress field caused by modulus mismatch and lattice mismatch between precipitates and matrix, are not limited to understanding precipitate hardening.They can also be used to analyze the stress field of those defects with modulus mismatch relative to the matrix, such as voids, bubbles.

Eshelby's inclusion theory without modulus mismatch
Given the common negligible influence of the hardening effect arising from the modulus mismatch between precipitates and the matrix (Santos-Güemes et al. 2018;Takahashi and Ghoniem 2008), an assumption is generally made that the elastic modulus of precipitates is the same as that of the matrix.For the ellipsoidal inclusion, Eshelby (1957) has given an analytical solution to the stress field inside and outside the inclusion (Eshelby 1957(Eshelby , 1959(Eshelby , 1961)).In Eshelby's inclusion theory, eigenstrain (or called transformation-strain) inside these inclusions refer to the nonelastic strains caused by thermal expansion, lattice misfit, phase transformation and so on.The eigenstrain ε T is assumed to be constant inside the inclusion, and be zero outside.
Taking spherical inclusion as an example, the stress field caused by inclusion can be written in terms of the auxiliary symmetric tensor where I is the second order unit tensor, G , ν represent the shear modulus, Poisson's ratio of matrix and precipitate, respectively.The internal stress field is (Mura 1987) where constant K = 5ν−1 15(1−ν) and H = 4−5ν 15(1−ν) .On the other hand, the stress field of a point x outside the inclusion of radius R is where x ′ = x − C , C is the center of the inclusion and d = x ′ .
The eigenstrain caused by lattice misfit (Eq.( 2)) is typically a dilatational misfit strain, while the eigenstrain caused by phase transform is more complex (Liu et al. 2017).The theoretical solution of internal and external stress field caused by lattice misfit is given as follows.Substituting Eq. (2) into Eqs.( 15) and ( 16), the simplified stress tensor inside the inclusion is The stress inside the particle, therefore, is purely hydrostatic, which can produce climb forces but no glide forces.The stress field outside the inclusion due to lattice misfit is (14 where x, y, z are the Cartesian components of x ′ , as shown in Fig. 10a. For instance, we select δ = 0.6% and draw the stress contours on different cross- sections of a spherical precipitate (z = ± 0.62R, z = ± R and z = 0), as shown in Fig. 10.Here R is the radius of the spherical precipitate.The origin of the coordinates is located at the spherical particle center.It is worth noting that, regardless of the value δ , the shear stress on the plane generated by the precipitates is always zero when the dislocation slips on the equatorial plane (z = 0) of the spherical precipitate along the X axis direction, such as shown in Fig. 10c (Chatterjee et al. 2021).Therefore, the contribution of coherent strengthening to hardening cannot be described by the obstruction of dislocation motion on the equatorial plane.While for the rest of the precipitate (z > 0 or z < 0), dislocations are attracted and repelled at two sides of the precipitate, respectively.At the corresponding positions (x, y, ±z ) on the cross-section at height ±z , the forces on dislocations are equal in magnitude and opposite in direction.For example, the stress field on the cross-section of the same position (z = ± 0.62R, z = ± R) of the positive and negative hemispheres is opposite, and the sign of the stress field depends on the sign of the lattice misfit δ of the precipitate, as illustrated in Fig. 10c.The hindrance to dislocation movement is equivalent on both planes, which causes dislocations easily pinned in the precipitate.This ultimately results in a comparable increase in the macroscopic yield strength.The radius r = π R/4 of the section at the height z = ± 1 − (π/4) 2 R ≈ ±0.62R from the equatorial plane of spherical precipitate just corresponds to the average radius of the cross-sections in the particle intersected by any plane.Therefore, the coherent strengthening is generally described by the results corresponding to the average intersection plane at z = ± 0.62R.Models based on slip resistance only captures interactions of dislocations gliding in planes that directly cross the precipitate.These cases are only appliable for a regular array of precipitates.With a random distribution, it is more complex because the dislocation is influenced by the stress fields of the different precipitates, which may be or not be intersected by the glide plane of the dislocation.Many researchers utilized Eshelby's inclusion theory to straightly integrate the stress field arising from lattice misfit between the precipitates and the matrix into DD framework.For instance, Mahler et al. (Mahler et al. 2021) investigated the collective impact of M 23 C 6 and MX precipitates, voids, and dislocation loops on the hardening behavior of Eurofer97 ferrite-martensitic steel.Liu et al. (Liu et al. 2021) considered the effects of dislocation climb on the interaction between dislocations and precipitates at high temperatures with the athermal mechanisms, where the stress field generated by lattice misfit was introduced into DD using the Eshelby method.Cui et al. (Cui et al. 2022) developed a multiscale framework to study the dislocation evolution during additive manufacturing, where Eshelby method is used in the DD part to study the effect of introducing nanoparticles on mitigating cracks of additively manufactured tungsten.

Eshelby's inclusion theory with modulus mismatch
When the effect of modulus mismatch between precipitate and matrix cannot be ignored, the stress field generated by precipitates will rely on the external stress, making it more complex compared with the case neglecting the modulus mismatch.The precipitate with modulus mismatch is typically called an inhomogeneity.In this case, a uniform applied stress at infinity introduces non-uniform stress field in the neighborhood of the inhomogeneities.Eshelby (Eshelby 1957) for the first time pointed out that the stress disturbance in applied stress due to the presence of inhomogeneity can be simulated when the equivalent eigenstrain of the precipitate is chosen properly.This equivalency is called the equivalent inclusion method (Mura 1987).
Here we briefly introduce the theory of the equivalent inclusion method.Consider an infinite homogeneous material D with the elastic modulus C ijkl containing an ellipsoidal domain with the elastic modulus C * ijkl . is called an ellipsoidal inhomogeneity.The stress and strain perturbations are denoted as σ ij and ε ij , respectively, when D is sub- jected to an external stress σ 0 ij at infinity.According to Hooke's law, the stress-strain relationship for both the interior and exterior of the inclusion is established where σ 0 ij = C ijkl ε 0 kl .The equivalent inclusion method is used to simulate the stress disturbance using the eigenstress resulting from an inclusion that occupies the space .Consider an infinitely extended homogeneous material with the elastic modulus C ijkl everywhere, containing a domain with an equivalent eigenstrain ε * kl to simulate the inhomogeneity problem.Then, Hooke's law yields (19) The necessary and sufficient condition for the equivalency of the stresses and strains from Eqs. ( 19) and ( 20) is, If σ 0 ij is a uniform stress, ε * kl is also uniform in .ε kl = S klmn ε * mn , where S klmn is the component of Eshelby tensor depending on the shape of precipitates (Eshelby 1957).Further, the equivalent eigenstrain of inclusion ε * kl can be determined according to the following equation, so as to solve the stress field inside and outside the inclusion.
Based on the above method, Eshelby has solved the theoretical solution of ellipsoidal inclusion in the isotropic matrix (Mura 1987).When the precipitate has eigenstrain ε p ij itself due to misfit or phase transformation, the inhomogeneity is called inhomogeneous inclusion.Then ε p ij should also be considered in the equivalent inclusion method depicted above (Eqs. (19,20,21 and 22)).The readers are suggested to see (Mura 1987) for more related details.

Numerical methods for calculating misfit stress field
Numerical methods for calculating the stress field are effective in handling intricate precipitate shapes, which is combined with the DD simulation to analyze precipitation hardening.Studying the stress field of a precipitate with arbitrary shape (and/or a modulus mismatch) is intricate.Whilst Eshelby's inclusion theory proves beneficial in resolving the corresponding force field of most precipitations, particularly of spherical inclusions with the same modulus as the matrix (Mura 1987), its effectiveness for solving the inhomogeneous inclusion problem is restricted to ellipsoidal inclusions, where the intrinsic strain within such inclusions is uniform as depicted in "Eshelby's inclusion theory without modulus mismatch" and "Eshelby's inclusion theory with modulus mismatch" sections.
Shin et al. determined appropriate boundary conditions and governing equations for precipitates with different matrix modulus, and solved the stress fields of spherical, cubic and cylindrical precipitates by using the Finite element method (FEM) (Shin et al. 2003(Shin et al. , 2005)).Takahashi and Ghoniem (Takahashi and Ghoniem 2008) identified the remaining challenges when solving the stress field generated by precipitate in 3d: high computational costs due to the requirement of fine mesh around precipitates, and the frequent use of C0 continuous condition in FEM causing a discontinuity in internal stress among adjacent elements.To address the aforementioned challenges, the boundary element method (BEM) (Ghoniem et al. 2000), which requires only boundary discretization, is employed to calculate the stress field and combined with parametric dislocation dynamics.The BEM for discretization of a spherical precipitate is shown in Fig. 11, where different boundary and volume element mesh sizes will result in different (20) spatial resolutions.This method can be utilized to deal with complex shapes of precipitates, by which the calculated stress field meets the continuity condition (Takahashi and Ghoniem 2008).
Another well-used method to solve the stress field is the Fast Fourier Transform (FFT) algorithm.When the slenderness ratio of the precipitate is very large, employing the finite element method for discretization could substantially increase the amount of computation required to achieve precise results.Therefore, Bertin et al. put forward the FFT method to solve the stress field and the boundary value problem of the periodic boundary for both homogeneous and heterogeneous materials (Bertin et al. 2015;Bertin and Capolungo 2018).Here, the mechanical equilibrium equations (Eq.( 23)) will be solved under the periodic boundary conditions within simulation domain V where, C denotes the fourth order elasticity tensor, ε the total strain, ε p the plastic strain, ε T the eigenstrain (or transformation-strain) associated with the precipitate and ∂V stands for the boundaries of domain V with normal n and E is the imposed macroscopic strain.δ(x) is a Dirac delta function that is equal to 1 when x is within the precipitate and 0 otherwise.More details about the FFT solver can be found in (Bertin et al. 2015;Bertin and Capolungo 2018;Santos-Güemes et al. 2018, 2022;Kohnert and Capolungo 2021).

Multiscale model based on energy
Introducing the interaction energy between precipitation and dislocation is another way of describing the resistance effect of the precipitate on the dislocation movement.By calculating the interaction energy, we can outline a continuous force field in the short range.Researchers can employ different multiscale modelling methods to obtain various forms of energy depending on the hardening mechanism of interest, as illustrated in Fig. 12 (Chatterjee et al. 2021;Wu et al. 2022;Lehtinen et al. 2016Lehtinen et al. , 2018)).Energy fitting based method and Eshelby's inclusion theory based method rely on MD for calibration of some coefficients, by which the force can be calculated and embedded in DD simulation.Atomistic simulations (MD or DFT) can also provide the energy to the larger scale to calculate the interaction force.Apart from the above method, stacking fault energy can also be incorporated into the phase field study to investigate the shearing of γ ′ pre- cipitates in Ni-base superalloys (Vorontsov et al. 2010(Vorontsov et al. , 2012)).Nonetheless, dealing with (23) 11 Boundary element meshes for a spherical precipitate with different resolutions (Takahashi and Ghoniem 2008) the intricate interactions between dislocations within phase field framework is difficult.Thus, this review does not discuss the phase field method.Three categories of methods will be elaborated as follows.

Energy fitting method for incoherent precipitates
In order to model a short-range interaction between the precipitate and the dislocation, the most practical approach is to determine the interaction energy form by fitting the DD results to the corresponding results from the microscopic scale, such as MD.By using interaction energy fitting method, simulation efficiency can be achieved on a mesoscale level (DD) while still incorporating the key details of the dislocation-precipitate interaction in a relatively fundamental way.
The form of interaction energy may vary, usually indicating a decay of interaction force with distance (Xiang et al. 2004).For the Gaussian Potential is selected for simplicity (Lehtinen et al. 2016).The parameters required in DD can be obtained from MD, including Gaussian potential, velocity of dislocation, shear modulus, and dislocation core energy (Lehtinen et al. 2016(Lehtinen et al. , 2018)).Specifically, the Gaussian potential U (d) = Ae −d 2 /R 2 is used to consider the interaction between the dislocation segments and the spherical immobile precipitate.The resulting interaction force applied to the discretization nodes connecting each dislocation segment inside a cutoff radius ( R cut−off , see Fig. 13) is expressed by where the parameter A characterizes the strength of the precipitate to hinder the dis- location motion.R is the radius of the spherical precipitate, and d is the distance of the dislocation from the center of the precipitate.
(24) The key parameter A incorporates many interaction details on a microscopic scale, which is obtained by comparing the critical stress τ c for different precipitate strengths A in DD simulations to the MD results, as shown in Fig. 14.The model does not contain the stress field solved based on Eshelby's inclusion theory, but it is suitable for simulating the local pinning effect caused by the incoherent precipitates, either strong or weak obstruction (Lehtinen et al. 2018).In addition, the Gaussian potential parameter A can indirectly represent the impact of the modulus mismatch between the precipitates and the matrix on the hardening, as MD simulation considers the modulus mismatch.Therefore, the effect of hardening of precipitates with different sizes, modulus, and distributions can be calculated in a relatively simple way by this method.
There are still some deficiencies that need to be addressed.While this approach is effective, it is somewhat approximate in reducing the interaction to only one single parameter.The interaction between edge dislocations or screw dislocations may be different, and requires further examination to make it more accurate.

Energy method based on Eshelby's inclusion theory
Interaction energy method is also applied to analyze the hardening effects of coherent precipitates.Similar to the misfit stress field method in "Eshelby's inclusion theory with modulus mismatch" section, the interaction energy between the precipitate and the dislocation can be calculated based on Eshelby's inclusion theory.It is established that Eshelby's inclusion theory simplifies acquiring the misfit stress field created by precipitates.However, the stress field of inhomogeneities (precipitates with modulus mismatch) still presents challenges.In that case, it is feasible to obtain the interaction energy and apply force on the dislocations as an alternative method.Cross-slip and climb can also be considered in this framework.
The procedures of interaction energy method for coherent strengthening and modulus strengthening are as follows.The first step to calculate interaction energy is to get the equivalent eigenstrain (see ε * kl in Eq. ( 22) in "Eshelby's inclusion theory with modulus mismatch" section).The interaction energy between dislocation and precipitate E pre,dis int equals to the total energy E pre + dis when the external load and a precipitate coexist, subtracting the energy E dis when only the dislocation exists and the energy E pre when only the precipitate exists, i.e., E pre,dis int = E pre + dis − E dis − E pre .Further, considering the stress field of dislocations σ D and the equivalent eigenstrain, the interaction energy can be calculated as Figure 15 sketches a precipitate with radius R whose center is placed at distance d from a straight dislocation with deviating angle θ from the slip plane of the dislocation.Consider the interaction force between a precipitate located at (x, y, 0) and an infinite small dislocation segment located at (0, 0, z) with length δz on the dislocation as shown in Fig. 15.The interaction force acted on the infinitesimal dislocation segment is given as the derivative of the energy change where F x is considered to simulate the resistance for dislocation on the slip plane in the framework of (Wu et al. 2022;Liu et al. 2023).F y perpendicular to the slip plane can also be incorporated in the same way in order to consider the cross-slip or climb.After some simplification, the nodal force can be obtained to implement in the DD simulations (Wu et al. 2022;Liu et al. 2023).
Note that Eshelby's inclusion theory alone does not suffice.It is imperative to integrate multiscale modelling since the dislocation stress field at the core region under classical elastic dislocation theory exhibits singularity.The interaction energy, as calculated by Eq. ( 25), becomes singular when a dislocation touches a precipitate (d = R).To overcome this singularity, MD simulation is employed to calibrate the spreading radius r a for iso- tropic materials (Cai et al. 2006).Under the given setup where the length of the initial edge dislocation, the volume of the simulation box, and the dislocation velocity in MD and DD are similar, calibration of r a is achieved by iteratively calculating critical stress by DD until it reproduces the MD simulated critical stress (Wu et al. 2022).

Stacking fault energy method
To deeply understand the mechanisms in the process of precipitates sheared by superpartial dislocations (Liu et al. 2006;Rai et al. 2021;Yang et al. 2018), a new method to compute generalized stacking fault forces is developed (Chatterjee et al. 2021).The motion of super-partial dislocations in precipitates is mainly governed by antiphase boundary (APB), which is one type of the generalized stacking energies (Long et al. 2017;Zhang et al. 2021).In DD framework, the basic idea of this method is to solve the force acted on the dislocations from the APB energy, which is similar to the case of defining critical stress due to antiphase boundary energy in "Case 1: DFT + DD" section (Long et al. 2017;Zhang et al. 2021).
The key to computing the forces exerted by generalized stacking fault on dislocations is to determine the stacking fault energy and save computational costs.Consider a slip plane in the material.Let γ (s(x)) be the generalized stacking fault energy density of that plane at point x, and s(x) be the slip vector at that point.Then the stacking fault energy Fig. 15 Schematic diagram of a precipitate interacting with an infinite straight dislocation lying on slip plane XOZ with dislocation line direction along Z-axis (Wu et al. 2022) of the plane is E SF = γ (s(x)) dA , where dA is an elementary area element of the glide plane.The force per unit length on the dislocation due to generalized stacking fault is therefore F =δ ĖSF /δw , where w is the dislocation velocity.Generalized stacking fault energy density can be obtained by DFT.For example, it could be fitted to the periodic function for nickel-based superalloys (Chatterjee et al. 2021) where s is a local (two-dimensional) slip vector on the glide plane, the k i are wave vec- tors, the M j are matrices representing the symmetries of the lattice, and S i , C i are the corresponding fitting constants.To solve the memory requirement issue due to the mesh resolution or storage and update of the slip vectors s(x) (Hussein et al. 2017), the concept borrowed from complex analysis and known as the winding number of a closed curve about a point (Krantz et al. 1999) is incorporated into the MoDELib (Mechanics Of Defect Evolution Library) code (Po and Ghoniem 2015).Figure 16 helps to understand how to calculate the stacking fault force acting on superpartials inside the precipitates and more details are referred to (Chatterjee et al. 2021).Then once the force is obtained, it can be summed to the Peach-Koehler force routinely computed in DD simulations.

Combining different types of multiscale models
Table 2 discussed the three types of multiscale models introduced in "Multiscale model based on slip resistance" -"Multiscale model based on energy" sections.Each method has limitations and may not be the solution to all problems.As illustrated in Table 2, each method could consider the effects of some mechanisms with ticks on precipitation hardening.For example, the method based on slip resistance typically focuses on slip ( 27 resistance for dislocation on a glide plane, but does not consider the dislocation motion outside the glide plane.This method is therefore not applicable to the analysis of crossslip or climb at high temperature.
Combining different types of multiscale models may be able to compensate for the shortcomings of these methods and thus improve their applicability.Here, we provide an illustrative case of how various approaches can be combined to reveal novel findings on precipitation hardening.
Using the stacking fault energy method in "Stacking fault energy method" section, (Chatterjee et al. 2021) conducted a thorough analysis of the interaction between spherical precipitates and a pair of edge super-partial dislocations in Ni-based superalloy, corresponding to order strengthening.It was found that the dislocation bypass mechanisms depend on the size ( r ) and volume fraction ( f ) of the precipitates, as displayed in Fig. 17.The study demonstrates that a hybrid looping-shearing process is possible and that it acts as the transition from shearing to looping regimes.However, it should be noted that this method solely takes into account the impact of dislocation when it penetrates the precipitates, disregarding the resistance outside the precipitates.
Combining the multiscale modelling methods based on misfit stress field in "Energy method based on Eshelby's inclusion theory" section with the stacking fault energy method in "Stacking fault energy method" section enhances its utility for comprehending long-range interactions arising from lattice misfit, in addition to the APB effects.This method is useful for understanding the effect of misfit on the bypass mechanisms shown in Fig. 17.The precipitate-dislocation interaction becomes more physically realistic, by considering the resistance for dislocation motion both inside and outside the precipitates.Figure 18 presents the impact of a misfit on the bypass stress in three cases: shearing, looping and hybrid shearing-looping regimes.It is suggested that the presence of dislocation pinning at precipitates becomes a characteristic indicator of strengthening by a lattice misfit (Chatterjee et al. 2021).

Multiscale modelling for precipitates with complex shapes
The multiscale modelling approach offers an effective solution for studying precipitates with complex shapes.Various shapes and its influencing factors have been discussed in "Basic information about the precipitates" section.Many macroscopic modelling studies ("Macroscopic models of precipitation hardening" section) tend to  oversimplify the precipitates, treating them as spherical particles with uniform size and distribution, while separately considering the contributions from different hardening factors.To show the effectiveness of multiscale modelling methods, we will present several such methods on precipitates with complex shapes, with or without modulus mismatch, as shown in Fig. 19.
(1) non-spherical precipitates without modulus mismatch Solving the eigenstrain of the precipitate with modulus mismatch presents a challenge in direct application to dislocation dynamics (DD) due to its dependence on the external stress field, as shown in "Eshelby's inclusion theory with modulus mismatch" section.In addition, it is demonstrated that modulus mismatch effect can be ignored compared to other effects in most cases (Nembach and Neite 1985;Takahashi and Ghoniem 2008).As a consequence, researchers typically disregarded the hardening effect due to modulus mismatch when studying the hardening behavior of precipitates.
Zhang and Sills (Zhang and Sills 2023) put forward a novel method to predict the hardening behavior of dislocations bypass plate-like precipitates, which allows for the comparison with the classical BKS model (Bacon et al. 1973) for spherical precipitates, where plate-like precipitates are simplified to high aspect ratio ellipsoids (Fig. 19a).Green's function is employed to solve the misfit stress field of the precipitate (Zhang and Sills 2023).This choice was motivated by the practical challenges in directly applying the analytical solutions derived from Eshelby's inclusion theory within a DD framework.The research suggests that achieving misfit strengthening is more attainable with plate-like precipitates than with spherical precipitates.
(2) non-spherical precipitates with modulus mismatch Shin et al. investigated the interaction between dislocations and spherical, cubic, and cylindrical precipitates with modulus mismatch, using DD in combination with finite element method (FEM), as depicted in Fig. 19b (Shin et al. 2003(Shin et al. , 2005)).In order to improve the limitation of FEM discussed in "Numerical methods for calculating misfit stress field" section, Takahashi and Ghoniem (Takahashi and Ghoniem 2008) proposed the coupling of boundary element method (BEM) and Parametric Dislocation Dynamics (PDD).The coupled method of BEM and PDD is used to solve the stress field caused by modulus mismatch and investigate the modulus strengthening during the interaction between dislocations and precipitates.Figure 19c shows the interaction between the edge dislocation and the slipped copper precipitate.It has been discovered that copper precipitates lose some of their resistance to dislocation motion when they are sheared after being cut by the leading dislocations in a pileup.In addition, when the elastic modulus ratio between the precipitate and the matrix is greater than 3-4, a transition from the particle shear mechanism to the Orowan bypass mechanism occurs, which also depends on the size of the precipitate (Takahashi and Ghoniem 2008).
Substantial systematic research based on DD and fast Fourier transform (FFT) methods has been conducted by researchers on the hardening behavior of plate-like precipitates in Al-Cu alloys, as shown in Fig. 19d, including the unshearable precipitate Al 2 Cu (Santos-Güemes et al. 2018, 2020) and the shearable precipitate Al 3 Cu (Santos-Güemes et al. 2021).In the framework of the coupling of DD and FFT methods (Bertin and Capolungo 2018), it is possible to introduce effects induced by coherent strengthening and modulus strengthening, and it is capable of dealing with precipitates of complex shapes in a relatively simple way.

Conclusions
This paper presents a comprehensive analysis of macroscopic theoretical model and multiscale modelling framework in terms of understanding precipitation hardening.The physical scenarios of dislocation-precipitate interaction process are first discussed, including particle shearing, Orowan bypass, cross-slip, and climb.Both intrinsic and extrinsic features of the precipitates, such as interface nature (coherent, semi-coherent, or incoherent) of precipitates, precipitate size and temperature, determines the dominated interaction mechanisms.
Regarding macroscopic theoretical models, an in-depth analysis of the classical linear tension based theoretical model is carried out.We quantitatively analyze how the dominant hardening mechanism changes as the radius of precipitates varies in a mechanism phase diagram.
As for multiscale modelling for precipitation hardening, we categorize recent research methods into three main types: slip resistance based approaches, misfit stress field based approaches, and energy based approaches.By integrating multiscale modelling with the physical scenarios, we systematically assess the key ideas of multiscale coupling, their applicability, advantages, and limitations.Furthermore, some examples are shown which combine different types of methods to mitigate the shortcomings of a specific type of multiscale model.The modelling of precipitates with complex shapes is also addressed as applications.
This review summarizes the understanding of precipitation hardening mechanisms, enables more precise predictions of hardening behavior, and informs the alloy design strategies by developing advanced multiscale methods.The basic ideas of the multiscale models discussed in this study can also be applicable to comprehend and predict other types of hardening in alloys caused by voids, bubbles, inclusions and so on.

Fig. 1
Fig. 1 Classification diagram of different precipitation hardening mechanisms incorporating planar and non-planar motion of dislocations.The blue solid lines in the figure represent the dislocation lines, and the blue planes are slip planes.The bold blue arrows indicate the gliding directions of dislocations, and the orange arrows indicate the Burgers vectors of dislocations.The grey circles correspond to precipitates

Fig. 2 a
Fig. 2 a and b shows the origin of the misfit strain of a coherent precipitate.a p , a are the lattice constant of precipitate and matrix, respectively.c and d shows the origin of misfit strain for an incoherent precipitate (no lattice matching) (Porter 2009)

Fig. 5
Fig.5Precipitation hardening mechanisms ( �σ OS , �σ CS + �σ MS and �σ Orowan ) as a function of coherent precipitate radius(Cui and Cui 2024).If the antiphase boundary energy γ APB is increased, the phase diagram illustrating the transformation of the hardening mechanism could shift from (a) to (b).The blue region corresponds to the range dominated by dislocation shearing mechanism, while the pink one represents the range dominated by Orowan bypass mechanism.The bold blue line signifies the dominant mechanism.In (a), R c1 is the critical radius of transition from order strengthening (OS) to coherent strengthening (CS) plus modulus strengthening (MS) as the dominant mechanism.R c2 is the critical radius of transition from dislocation shearing mechanism (blue region) to Orowan bypass mechanism (pink region) as the dominant mechanism.In (b), R c represents the critical radius for the shift from OS to the Orowan bypass mechanism as the dominant mechanism

Fig. 6
Fig. 6 Illustration of a computational approach to the multiscale modelling of precipitation hardening based on slip resistance methods.The three main approaches to analyze the hardening effects are case 1: Density Functional Theory (DFT) + Dislocation Dynamics (DD), case 2: Molecular Dynamics (MD) + DD, and case 3: Crystal Plasticity-Finite Element Method (CP-FEM) + DD or Macroscopic models.γ APB in the figure is antiphase boundary energy.τ c is critical shear stress, and τ (α)

Fig. 9
Fig. 9 Illustration of computational approaches to the multiscale modelling of precipitation hardening based on misfit stress field.The two main approaches to solving the stress field embedded in the DD framework are: theoretical solutions based on Eshelby's inclusion theory, and numerical solutions based on FEM/BEM or FFT solver

Fig. 10 a
Fig. 10 a Schematic diagram of the intersection of dislocation slip plane and precipitate.The blue solid line and the grey sphere in figure represent a dislocation line and a spherical precipitate, respectively.The shallow blue plane corresponds to the dislocation slip plane, by which the cross-section in the precipitate intersected is marked as the blue elliptical dotted line.b Spherical precipitation intersects with dislocation slip planes at different heights from the equatorial plane (z = 0, z = ± 0.62R and z = ± R) shown in shallow blue ellipses.c Contour plots of stress component σ 13 /G on the cross-section of the spherical precipitate intersected by dislocation slip planes at different heights ( δ = 0.6% ), corresponding to the slip planes in (b)

Fig. 12
Fig. 12Illustration of a approach to the multiscale modelling of precipitation hardening based on energy.The three main approaches to analyze the hardening effects are: (a) energy fitting based method, (b) Eshelby's inclusion theory based method, and (c) energy calculated from atomistic simulations (DFT or MD) embedded in DD.Sections 6.1, 6.2 and 6.3 in the figure correspond to "Energy fitting method for incoherent precipitates", "Energy method based on Eshelby's inclusion theory", and "Stacking fault energy method"

dFig. 13 a
Fig. 13 a Top view and (b) front view of the schematic diagram of dislocation-precipitate interaction implementation in ParaDis.A force generated by the precipitate is applied to the discretization nodes inside the cutoff radius R cut−off(Lehtinen et al. 2016)

Fig. 16
Fig. 16 Method of calculation of the stacking fault force acting on a/2�110�{111} super-partial dislocations inside a γ ′ precipitate.The super-partial dislocation loop (loop i) has Burgers vectors b i = a/2[110] .Minus Burgers vectors correspond to equivalent slip vectors in γ-surface.Here PC represents the perfect crystal (Chatterjee et al. 2021)

Fig. 17 Fig. 18
Fig. 17 Tree diagram of γ ′ precipitates bypass mechanisms by paired a/2�110�{111} dislocations in nickel-based superalloys.The leading and trailing dislocations are labeled as L and T respectively.Depending on the precipitate radii and volume fraction, L either shears (A) the precipitates or loops (B) around them.If L shears, T always shears.If T is outside the precipitate and remains almost straight during the critical configuration, the bypass mechanism satisfies the weakly-coupled shearing conditions ((A1.1) and (A1.2)).If L and T are inside the same precipitate during the critical configuration, strongly coupled shearing occurs ((A2.1) and (A2.2)).If L loops, T can either form a second loop around the precipitates ((B1.1) and (B1.2)), or trigger the hybrid looping-shearing bypass mechanisms, where the leading loop is collapsed in the critical configuration (B2.1) forming APB, followed by trailing shearing which removes the APB (B2.2)(Chatterjee et al. 2021)

Fig. 19
Fig. 19 The diagram of different multiscale modelling methods for simulating intricate precipitate shapes.a plate-like precipitate without modulus mismatch (Zhang and Sills 2023).b finite element method (Shin et al. 2005), c boundary element method (Takahashi and Ghoniem 2008) and (d) Fast Fourier Transform (Santos-Güemes et al. 2020) could study precipitates with arbitrary shapes with modulus mismatch

Table 2
Three main multiscale modelling methods (resistance based, stress field based and energy based methods) can be used to analyze the mechanisms with ticks in the table