Crystal plasticity finite element simulation method based on dislocation climb mechanism

Abstract

The application discloses a kind of crystal plastic finite element simulation methods based on dislocation climb mechanism, establishes the crystal plastic finite element model (CPFEM) based on dislocation, studies the relationship between microstructure and creep behavior.The model considers the dislocation climb caused by vacancy thermal motion, and considers that dislocation climb will affect the attempt frequency of dislocation bypass mechanism and the evolution of dislocation density.After comparing experimental and simulation results, the results show that the crystal plastic finite element model can better describe the creep response of the material at high temperature.

Description

Technical Field

[0001] This invention relates to the field of high-temperature creep mechanical behavior simulation and prediction technology, and in particular to a crystal plasticity finite element simulation method based on dislocation climb mechanism. Background Technology

[0002] Energy utilization has long been a key challenge for human development. Researching and utilizing energy, and continuously improving the utilization rate of new energy sources, has always been an important research topic for countries worldwide. With the rise of a new industrial revolution, addressing the environmental impacts of climate change has increasingly become a focus of global attention, and green, low-carbon, and sustainable development has become a primary goal of energy research. As research into energy development technologies continues, researchers have discovered a significant increase in energy conversion efficiency at high temperatures. Simultaneously, with continuous economic and technological development and progress, the service temperature of modern energy infrastructure structures is increasing year by year. In recent years, northern my country has begun to experience smog and extreme weather events, leading to growing emphasis on environmental protection and green development. Against this backdrop, my country has proposed building a green, safe, low-carbon, and sustainable modern energy system, which requires vigorous development of new energy technologies while improving existing energy utilization technologies. In the future, with the development of new energy research, installed capacity of new energy sources will gradually increase, and research on new energy units will occupy an increasingly important position.

[0003] Nuclear energy, as a safe, economical, and efficient clean energy source, plays a unique role and is of great significance in implementing the new development philosophy, building a new development pattern, and leading the world to achieve a green and low-carbon recovery in the post-pandemic era. With the continuous development of nuclear energy, nuclear power, along with hydropower and coal power, now constitutes one of the three pillars of the world's energy supply. At the same time, nuclear energy is also one of the important foundations for meeting energy supply needs and ensuring national security. As of 2019, nuclear power accounted for 10.4% of global electricity generation. Globally, 449 commercial nuclear power reactors were in operation in 30 countries, with a total installed capacity of 396 GW. There were 55 nuclear power units under construction, with an installed capacity of 57 GW. China's nuclear power plant development began in the early 1980s. By the end of 2019, 47 nuclear power units were in operation, with a total capacity of 48 GW, second only to the United States and France. There were 13 units under construction, with a total installed capacity of 13 GW, ranking first in the world in terms of installed capacity.

[0004] Historically, nuclear power development can be divided into four periods. The first generation of nuclear power plants consisted of early prototype reactors, primarily developed in the 1950s and 60s, with light water reactors being the main representative. During this decade, approximately 40 nuclear power units were put into operation globally, such as the UK's Magnox graphite gas-cooled reactor, the US's Shiping Harbor pressurized water reactor, and the Dresden boiling water reactor. The second generation of nuclear power plants, mainly developed from the late 1960s to the 1990s, consisted of large-scale commercial nuclear power plants built upon the foundation of the first generation. This period saw rapid development in the nuclear power industry, coinciding with the global oil crisis and the increasing economic viability of nuclear power. Over 240 units were put into operation in 20 years, including the Canadian reactor and the Soviet pressurized water reactor. By the early 2010s, most of the world's operating commercial nuclear power plants belonged to the second generation. The third generation of nuclear power plants mainly refers to advanced light water reactor nuclear power plants, those that began operation from the late 1990s to the early 21st century. This generation of nuclear power plants offers better economics and higher safety performance, mainly including advanced boiling water reactors and European pressurized water reactors. Fourth-generation nuclear energy systems are nuclear power plants under development, with the goal of building advanced nuclear energy systems that are ready for practical use by 2030. The main advantages of fourth-generation nuclear energy systems are higher economics, better safety performance, sustainable development, and effective prevention of nuclear proliferation.

[0005] Generation IV nuclear energy systems represent the technological forefront of advanced nuclear energy systems and are the future direction of nuclear power development. Generation IV nuclear energy systems mainly include six reactor types: lead-bismuth cooled fast reactor systems, molten salt reactor systems, sodium-cooled fast reactor systems, supercritical water reactor systems, ultra-high temperature gas-cooled reactor systems, and gas-cooled fast reactor systems. Lead-bismuth cooled reactors, as one of the candidate reactor types for Generation IV nuclear energy systems, have promising development prospects, and research on lead-bismuth cooled reactors is extensive worldwide. The latest information released by the Generation IV International Forum indicates that lead-bismuth cooled fast reactors are expected to be the first Generation IV nuclear energy system to achieve commercial application. Among the reactor components of Generation IV nuclear energy systems, the operating conditions of fuel element cladding are the most demanding. Taking lead-bismuth cooled fast reactors as an example, because the fuel element cladding of lead-bismuth cooled fast reactors operates in a high-temperature environment, significant creep occurs, becoming a crucial evaluation factor threatening the safety of lead-bismuth cooled fast reactors. Therefore, conducting research on the creep phenomenon of cladding materials at high temperatures and analyzing the creep evolution law is of great reference value and guiding significance for solving the problem of material selection for the cladding structure of lead-bismuth cooled fast reactors.

[0006] The core of creep models is to obtain the stress-strain relationship of materials in the study and analysis of high-temperature structural materials. Since Norton (1929) and Bailey (1935) established the creep constitutive equation, scholars have developed hundreds of creep models to describe the curves of creep strain changing with time, stress, temperature, and other factors. The simplest creep equation can be expressed as follows:

[0007] ε=Kσ n (1-1)

[0008] ε=Bt m (1-2)

[0009] Where σ and ε are creep stress and creep strain, respectively; t is creep time; m and n are exponential constants; and K and B are material parameters.

[0010] Equations (1-1) and (1-2) mainly reflect the power-law relationship between creep strain and stress and time. The creep-time curve divides the creep process into three stages, and the creep strain rate is correspondingly divided into three stages: the first stage, where the creep strain rate decreases with increasing time, is called the creep deceleration stage; the second stage, where the creep strain rate remains essentially constant, is called the steady-state creep stage; and the third stage, where the creep strain rate continuously accelerates with the accumulation of creep damage, is called the accelerated creep stage. The creep strain rate proposed by Norton satisfies the following equation:

[0011]

[0012] Where Q is the thermal activation energy; R is the Boltzmann constant; T is the temperature; and a is the material constant.

[0013] Subsequently, Johnson et al. revised the above relationship in 1963:

[0014]

[0015] Where b, n1, and n2 are constants characterizing the material properties; Q A and Q B These represent different activation energies.

[0016] In 1995, Bartsch further improved Norton's equations by introducing an exponential correction term for stress:

[0017]

[0018] Where a1, a2, b1, b2, m1, and m2 are material constants used to correct for the effects of stress and time on creep evolution.

[0019] The microscopic process of creep fracture manifests as the nucleation and growth of pores and microcracks. Studying the fracture behavior of creep from a microscopic perspective aligns with its physical essence. Because metallic crystals inherently contain numerous defects, the microscopic inhomogeneity of the matrix material distribution leads to localized stress concentrations in the structure. Furthermore, the formation and growth of pores vary across different parts of the material, making it difficult to accurately represent the overall fracture process using the evolution of pores in localized areas. Therefore, accurately describing the overall microscopic evolution of the material is a key challenge in creep research. By considering the microscopic evolution of materials during fracture from the perspective of fracture mechanics, establishing a mechanical model of temperature, stress, and crack propagation can be well applied to engineering practice. The ASME standard introduces a parameter related to stress and temperature to establish a creep strain rate evolution equation:

[0020]

[0021] Where Ω represents the material parameter; The initial creep strain rate is denoted as .

[0022] Considering the evolution of microscopic inhomogeneities and discontinuities (such as micropores and microcracks) in materials, microscopic models can be established to describe them. However, the direct calculation of microscopic models introduces enormous complexity, which makes them lack the simplicity and efficiency of phenomenological models. To account for the complexity brought about by the microscopic evolution mechanism, damage mechanics has been introduced into creep modeling.

[0023] With in-depth research on creep damage, damage theory has gradually developed. Kachanov first proposed the concept of continuity and gave the evolution relationship of creep damage, realizing the prediction of creep life. Subsequently, Rabotnov established the Kachanov-Rabotnov damage equation by introducing the concept of effective stress, which promoted the development of creep damage constitutive equations into a new stage. Its main equations are as follows:

[0024]

[0025]

[0026] Where ω and These represent damage evolution parameters and damage evolution rates, respectively; A, C, q, p, and These are material parameters.

[0027] During the service life of high-temperature components, the formation, growth, and connection of intergranular pores are the main creep damage mechanisms. Furthermore, the location of creep damage exhibits localized characteristics; for example, pores form in the weld at welded joints, leading to intergranular fracture. Although creep damage is localized, considering this locality within phenomenological creep damage theory remains a difficult problem. In other words, the locality of creep damage is not reflected in the aforementioned creep damage constitutive equations. Liu Yan, employing the basic principles of mesomechanics, proposed an improved KR model. This model treats the material as a composite of damaged and undamaged components at the microscale, while still considering it as a continuous medium at the macroscale. A simple volume-averaging method is used to construct the damage constitutive equations for macroscopic quantities.

[0028] For a uniaxial creep tensile specimen under high temperature conditions, the total strain consists of three parts: elastic strain, plastic strain, and creep strain.

[0029]

[0030] Where ε e , ε p , ε c These are the elastic strain tensor, plastic strain tensor, and creep strain tensor, respectively. The evolution of elastic strain follows Hooke's Law:

[0031] σ=C:ε e (2)

[0032] Where C is the elastic tensor, which can be represented by a 6×6 matrix as follows:

[0033]

[0034] In the matrix, λ and μ are the first and second Lamé constants, respectively, which can be represented by the elastic modulus E and Poisson's ratio v:

[0035]

[0036]

[0037] For plastic strain, considering that the applied stress is small under creep loading, plastic strain can be ignored; creep strain, however, requires a specific evolution equation to describe it.

[0038] The improved KR model is as follows:

[0039]

[0040]

[0041] D cr = 1 - (1 - g) (1 / (φ+1)) (8)

[0042] Where σ e Mises equivalent stress; S is the creep strain tensor; ij Deviatoric stress tensor; D damage variable; σ1 is the first principal stress; D cr The critical damage value is when D = D cr When , it indicates that the material has reached its creep failure life. B, n, A, m are material constants related to the minimum creep strain rate and fracture; φ is the damage constant; ρ, g are constants that consider damage non-uniformity, where ρ is the volume occupied by all damage elements in the material, and g comprehensively reflects the types and proportions of damage elements in the material; α is a material constant characterizing the multiaxial failure criterion.

[0043] While numerous macroscopic creep models (including the Norton model and the improved KR model) have successfully predicted creep response trends, they do not consider the correlation between mechanical response and underlying deformation mechanisms. Currently, an increasing number of high-temperature resistant materials are being used in specialized high-temperature environments, such as cladding materials for lead-bismuth fast reactors. In such high-temperature environments, the diffusion of vacancies and interstitial atoms within the metallic material, and the resulting dislocation climb, play crucial roles in the material's plastic deformation. Furthermore, dislocation climb is considered a key mechanism of creep in metallic materials, significantly impacting equipment performance and service life under high-temperature conditions. Therefore, it is necessary to establish constitutive models that more physically meaningfully consider dislocation mechanisms.

[0044] In recent years, crystal plasticity models have been used to describe the creep behavior of various materials. Unlike classical plasticity theory, crystal plasticity theory introduces the concept of crystal slip systems and posits that the plastic deformation of crystalline materials mainly originates from shear strain on activated slip systems, thus naturally capturing the anisotropic mechanical behavior of crystalline materials. These models are based on dislocation slip theory, describing the microscopic plastic deformation mechanism of materials through the concept of dislocations. For example, Staroselsky et al. established a creep-plastic composite constitutive relation to describe the low-cycle and thermomechanical fatigue behavior of single-crystal Ni. While this creep model can well describe the creep fatigue response of materials, its constitutive relation is phenomenological and cannot reveal the dislocation climb phenomenon at high temperatures. Yuan et al. further established a crystal plasticity model of face-centered cubic single-crystal slip / climb coupling, in which the authors considered the contribution of dislocation climb. Through this model, they revealed and elucidated the evolution process of dislocations during creep. However, this model only describes the creep response of single-crystal materials and did not compare it with experimental results. Summary of the Invention

[0045] To address the problems existing in the prior art, the purpose of this invention is to provide a crystal plasticity finite element simulation method based on dislocation climb mechanism. This invention is mainly used for the simulation and prediction of the high-temperature creep mechanical behavior of key structural materials in lead-bismuth fast reactors.

[0046] To achieve the above objectives, the technical solution adopted by this invention is: a finite element simulation method for crystal plasticity based on dislocation climb mechanism, comprising:

[0047] In finite deformation, the overall deformation gradient tensor F is decomposed into the elastic part F by multiplication. e and the plastic part F p :

[0048] F = F e ·F p

[0049] In the above equation, the plastic deformation gradient tensor F p The mapping of material points from a reference configuration to an intermediate configuration corresponds to the plastic deformation caused by dislocation motion in the crystal plasticity model, and the elastic deformation gradient tensor F. e It describes the mapping of material points from intermediate configurations to the current configuration, corresponding to lattice elastic stretching and rigid body rotation in the crystal plasticity model;

[0050] intermediate configuration plastic velocity gradient L P Represented as:

[0051]

[0052] When dislocation slip is the only source of plastic deformation, L P It can be expressed as the shear strain rate on all slip systems of the crystal. Superposition:

[0053]

[0054] Where the unit vector m α and n α Describe the slip direction and slip plane normal of the edge dislocation in slip system α in the intermediate configuration; for the dislocation-dominant model, the shear strain rate on its α slip system. Represented as:

[0055]

[0056] Where v e This is the average velocity of dislocation motion in the crystal. The above equation can be rewritten in the following form:

[0057]

[0058] in It is the effective decomposed shear stress τ on the slip system α. sol It is the strength of the solid solution, b s v is the length of the Burgers vector of the slip, v0 is the velocity of the dislocation slip, v climb It is the speed of dislocation climb, L c The physical meaning of Q is the obstacle characteristic size of the climb bypass mechanism. This parameter controls the strength of the climb-enhanced mobility of edge dislocations. s It is the energy required to activate dislocation slip, k B is the Boltzmann constant, T is the temperature, and p and q reflect the fitting parameters of the barrier energy barrier;

[0059] Microstructure evolution is caused by edge dislocation density ρ e With dipole dislocation density ρ d Characterization. When positive and negative dislocations move to a certain range d α When dislocations are inside, they annihilate, causing a decrease in dislocation density; at the same time, when in a state of... When the range is wide, dislocation dipoles will form and will not exhibit polarity outwards;

[0060] The final evolution of the edge dislocation density is given by the following equation:

[0061]

[0062] in The mean free path of dislocation slip; the evolution of dislocation dipole density is given by the following equation:

[0063]

[0064] choose The stress tensor is solved as the starting point for calculation. After the stress tensor is solved and updated, the variables at the microstructure level are iteratively updated.

[0065] As a further improvement of the present invention, effective shear stress is... The calculation is as follows:

[0066]

[0067] Where τ α The shear stress in the slip system α, and the resistance that dislocation slip needs to overcome. The definition is as follows:

[0068]

[0069] Where G is the shear modulus, ξ αα′ It is the mutual influence coefficient between different slip systems α and α′.

[0070] As a further improvement of the present invention, the average velocity v of dislocation motion in the crystal is... e The calculation method is as follows:

[0071]

[0072]

[0073] The velocity of dislocation climb, v climb The calculation method is as follows:

[0074]

[0075] Where D0 is the self-diffusion coefficient of the face-centered cubic structure, Ω is the volume activated by dislocation climb, and Q... c It is the activation energy for dislocation climb. When two dislocations can form a dislocation dipole, It is the maximum distance between the two sliding planes; d α It is the minimum distance between the annihilation of two edge dislocations; and d α The calculations are as follows:

[0076]

[0077] d α =C anni b s

[0078] Where C ann i is a fitting parameter. Strain hardening is achieved through the mean free path of dislocations. To describe:

[0079]

[0080]

[0081] Where d is the average grain size, ξ αα′ It is the mutual influence coefficient between different slip systems α and α′.

[0082] The beneficial effects of this invention are:

[0083] 1. This invention is mainly used for the simulation and prediction of high-temperature creep mechanical behavior of key structural materials in lead-bismuth fast reactors. Starting from the mesoscale, this invention uses the crystal plastic finite element method to study the essence of material plastic deformation, namely the influence of dislocation mechanism on the macroscopic creep response of materials, so that the model is closer to the real metal plastic deformation in a physical sense, providing theoretical support for the high-temperature creep application of engineering materials;

[0084] 2. This invention can effectively simulate and predict the high-temperature creep response of key structural materials in complex environments of lead-bismuth fast reactors, obtaining the steady-state creep rate under relevant conditions to predict material lifetime. Furthermore, compared to traditional phenomenological models that lack physical meaning, this model introduces the concept of dislocation density, considering the dislocation climb phenomenon caused by vacancy movement due to thermal disturbances in high-temperature environments. It considers the dislocation climb mechanism to be non-negligible, which is clearly observed in subsequent simulation results. Moreover, compared to recent creep crystal plastic constitutive models, this model can effectively simulate the high-temperature creep response of polycrystalline materials at different temperatures. Attached Figure Description

[0085] Figure 1 This is a deformation gradient evolution diagram in an embodiment of the present invention;

[0086] Figure 2 This is a diagram illustrating dislocation annihilation and dislocation dipole formation in an embodiment of the present invention;

[0087] Figure 3 This is a diagram illustrating the dislocation transition and dislocation climb mechanism in an embodiment of the present invention;

[0088] Figure 4 This is a clockwise closed-loop process for stress calculation in this embodiment of the invention;

[0089] Figure 5 This is a two-dimensional 50-grain Vorinoi model diagram from an embodiment of the present invention;

[0090] Figure 6 for Figure 5 A schematic diagram of the boundary conditions in a 2D 50-grain Vorinoi model;

[0091] Figure 7 This is a Mises stress cloud diagram of each grain in the embodiments of the present invention;

[0092] Figure 8 This is a dislocation density cloud map of each grain in the embodiments of the present invention;

[0093] Figure 9 The figures show the uniaxial tensile stress-strain curves of 316 austenitic stainless steel at different temperatures in the embodiments of the present invention.

[0094] Figure 10 This is a graph showing the yield strength and tensile strength of 316 austenitic stainless steel as a function of temperature in an embodiment of the present invention.

[0095] Figure 11 These are three types of sawtooth waves in the embodiments of the present invention;

[0096] Figure 12 This is a uniaxial tensile stress-strain curve diagram of T91 in an embodiment of the present invention;

[0097] Figure 13 This is a creep-time curve of 316 austenitic stainless steel 873K in an embodiment of the present invention;

[0098] Figure 14 This is a creep-time curve of 316 austenitic stainless steel 973K in an embodiment of the present invention;

[0099] Figure 15 This is a high-temperature uniaxial creep tensile simulation diagram of stainless steel at 773K in an embodiment of the present invention;

[0100] Figure 16 This is a high-temperature uniaxial creep tensile simulation of stainless steel at 873K in an embodiment of the present invention. Detailed Implementation

[0101] The embodiments of the present invention will now be described in detail with reference to the accompanying drawings.

[0102] Example

[0103] This embodiment establishes a dislocation-based crystal plasticity finite element model (CPFEM) to investigate the relationship between microstructure and creep behavior. The model considers dislocation climb caused by vacancy thermal motion, assuming that dislocation climb affects the attempt frequency of dislocation bypass mechanisms and the evolution of dislocation density. Experimental and simulation results are then compared, showing that the CPFEM model can effectively describe the creep response of materials at high temperatures.

[0104] Although dislocation slip in crystals is not continuous, it can be considered continuous macroscopically due to the large number of dislocation slips present. For the kinematics and constitutive equations of high-temperature creep, the classical crystal plasticity theory framework developed by Rice and Hill et al. is adopted. In finite deformation, the overall deformation gradient tensor F can be multiplicatively decomposed into the elastic part F0. e and the plastic part F p :

[0105] F = F e ·F p

[0106] In the above equation, the plastic deformation gradient tensor F p The mapping of material points from a reference configuration to an intermediate configuration corresponds to the plastic deformation caused by dislocation motion in the crystal plasticity model, and the elastic deformation gradient tensor F. e The mapping of material points from intermediate configurations to the current configuration corresponds to lattice elastic stretching and rigid body rotation in the crystal plasticity model, and its evolution is shown in Figure 1.

[0107] intermediate configuration plastic velocity gradient L P It can be represented as:

[0108]

[0109] When dislocation slip is the only source of plastic deformation, L P It can be expressed as the shear strain rate on all slip systems of the crystal. Superposition:

[0110]

[0111] Where the unit vector m α and n α Describe the slip direction and slip plane normal of the edge dislocation in slip system α in the intermediate configuration. For this dislocation-dominant model, the shear strain rate on the α-slip system is... This can be expressed using the Orowan formula:

[0112]

[0113] Where v e This is the average velocity of dislocation motion in the crystal, and the equation can be rewritten in the following form:

[0114]

[0115] in It is the effective decomposed shear stress τ on the slip system α. sol It is the strength of the solid solution, b s v is the length of the Burgers vector of the slip, v0 is the velocity of the dislocation slip, v climb It is the speed of dislocation climb, L c The physical meaning of Q is the obstacle characteristic size of the climb bypass mechanism. This parameter controls the strength of the climb-enhanced mobility of edge dislocations. s It is the energy required to activate dislocation slip, k B q is the Boltzmann constant, T is the temperature, and p and q are fitting parameters of the barrier energy barrier.

[0116] Effective shear stress The calculation is as follows:

[0117]

[0118] Where τ α The shear stress in the slip system α, and the resistance that dislocation slip needs to overcome. The definition is as follows:

[0119]

[0120] Where G is the shear modulus, ξ αα′ It is the mutual influence coefficient between different slip systems α and α′.

[0121] Because 316 austenitic stainless steel has a high stacking fault energy at high temperatures, dislocation slip is the primary mode of plastic deformation. Therefore, the contributions of crystal twinning and phase transformation to plasticity are not considered here. In this model, microstructure evolution is mainly determined by the edge dislocation density ρ. e With dipole dislocation density ρ d Characterization. When positive and negative dislocations move to a certain range d α When this happens, dislocations annihilate, causing a decrease in dislocation density. Simultaneously, when in a state of... When the range is wide, dislocation dipoles will form and will not exhibit polarity outwards, such as Figure 2 As shown.

[0122] like Figure 3 As shown, the evolution of the final edge dislocation density is given by the following formula:

[0123]

[0124] in Mean free path of dislocation slip The specific expression will be given later. The evolution of the dislocation dipole density is given by the following formula:

[0125]

[0126] As previously known, dislocation slip and dislocation climb are the main controlling modes of creep in 316 austenitic stainless steel at high temperatures. Assuming the velocity of dislocations in the crystal approaches infinity, the average velocity of dislocation slip mainly depends on the time the material waits before overcoming an obstacle. Dislocations overcome obstacles primarily in two ways: first, due to thermal activation, edge dislocations can directly leap over obstacles; second, through dislocation climb, when encountering an obstacle, edge dislocations move perpendicularly and eventually overcome it. Therefore, the contribution of material dislocation climb to dislocation mobility is also considered in the Orawaan model.

[0127]

[0128]

[0129] The velocity of dislocation climb is given by the formula:

[0130]

[0131] Where D0 is the self-diffusion coefficient of the face-centered cubic structure, Ω is the volume activated by dislocation climb, and Q... c It is the activation energy for dislocation climb. When two dislocations can form a dislocation dipole, It is the maximum distance between the two sliding planes; dα It is the minimum distance between two edge dislocations annihilating. and d α The calculations are as follows:

[0132]

[0133] d α =C anni b s

[0134] Where C anni It is a fitting parameter. Strain hardening is achieved through the mean free path of dislocations. To describe:

[0135]

[0136]

[0137] Where d is the average grain size, ξ αα′ It is the mutual influence coefficient between different slip systems α and α′.

[0138] It is important to note that in this model, dislocation barriers are uniformly simplified to forest dislocation barriers. This is because it is very difficult to consider all barrier types within a single model. Furthermore, forest dislocation barriers are only a conceptual barrier; they can also be equivalently viewed as dislocation walls, dislocation knots, or other barrier types.

[0139] The core of solving the crystal plasticity finite element method is calculating the stress response at the integration point under a given strain in each increment step. In the above description, the equations for plastic flow are based on the rate form of physical quantities, thus requiring time integration. Furthermore, the flow criteria also depend on the current state of the microstructure; the evolution equations for microstructure variables such as dislocation density, twinning, and martensite volume fraction are also rate-form, similarly requiring time integration. Assuming the starting point of an increment step is t0, after a small time increment Δt, the ending point is t = t0 + Δt. Given the deformation gradient, stress tensor, and microstructure variables at time t0, the stress tensor and microstructure variables need to be updated through integration based on the deformation gradient at the end point t. However, the crystal plasticity constitutive model developed in this paper is based on a complex physical mechanism and exhibits strong nonlinearity, making analytical calculation of the above integration process impossible. Therefore, a numerical iterative method is required, with the Newton-Raphson iterative method being a commonly used approach.

[0140] In crystal plasticity models, iterative calculations of corresponding forces typically use methods such as... Figure 6-16 The process shown can theoretically be completed by using any physical quantity in a clockwise closed loop as the starting point, but generally, it is preferred to start from... Begin. This is because from... Any physical quantity other than that will be involved The calculation in the crystal plasticity model These are often power functions or exponential functions. Through these types of functions, even a small change in S will cause... The drastic fluctuations result in poor iterative convergence. Conversely, from Starting from point S and ending at point S provides better iterative stability in stress calculations. Therefore, under large deformation conditions, from... It would be better to start by calculating convergence, but this would also lead to a larger computational cost because the dimension of the Jacobian matrix in the Newton-Raphson iteration method depends on the number of independent variables of the physical quantities being iterated. In the physical quantities, F... p There are 8 independent variables, F e There are 9 independent variables, and S has 6 independent variables. In comparison, There are 12 FCC crystals, resulting in a very large iteration matrix. Considering both of the above factors, ensuring computational convergence is more important than computational time; therefore, [the appropriate method is chosen]. As the starting point for calculation, such as Figure 4 As shown.

[0141] After the stress tensor is solved and updated, it is also necessary to iteratively update the variables at the microstructure level. This model adopts the two-level iterative method based on intermediate configuration proposed by Kaldindi et al.

[0142] Model building:

[0143] For polycrystalline simulation of creep in 316 austenitic stainless steel, a two-dimensional Vorinoi model with 50 grains was used, such as... Figure 5 As shown in Figure 6, the boundary conditions are set as follows: Assuming a quadrilateral ABCD, the normal displacements of sides AB and BC are constrained respectively, and all points on side CD except point D are coupled to point D, so that the load is applied to point D, which facilitates subsequent data extraction and post-processing.

[0144] Model Calculation: The calculation of this model utilizes the DAMASK crystal plasticity calculation platform developed by the Max Planck Institute for Steel Research in Germany to couple the constitutive model with the commercial finite element software ABAQUS. The calculation process requires three main files: XX.inp, material.config, and abaqus_v6.env. These represent the model file, material parameter file, and the file used by DAMASK to call ABAQUS for calculation, respectively. The model is the Vorinoi polycrystalline model established in Part 2, which can be created using the open-source software Neper. The Material.config file is primarily used to define the material parameters.

[0145] like Figure 7 and Figure 8 As shown, this model can effectively simulate and predict the high-temperature creep response of key structural materials in complex environments of lead-bismuth fast reactors, obtaining the steady-state creep rate under relevant conditions to predict material lifetime. Furthermore, compared to traditional phenomenological models that lack physical meaning, this model introduces the concept of dislocation density, considering the dislocation climb phenomenon caused by vacancy movement due to thermal disturbances in high-temperature environments. It considers the dislocation climb mechanism to be non-negligible, which is clearly observed in subsequent simulation results. Moreover, compared to recent creep crystal plastic constitutive models, this model can effectively simulate the high-temperature creep response of polycrystalline materials at different temperatures.

[0146] In high-temperature creep, the dislocation mechanism is the main plastic deformation mechanism. The role of dislocation climb caused by thermal activation becomes more and more obvious at high temperatures. Therefore, it is necessary to consider the changes in dislocation density and activation frequency caused by the dislocation climb mechanism.

[0147] In this model, the microstructure evolution is mainly driven by the edge dislocation density ρ e Characterized by the dipole dislocation density ρd. When positive and negative dislocations move to a certain range d... α When this happens, dislocations annihilate, causing a decrease in dislocation density. Simultaneously, when in a state of... Within a certain range, dislocation dipoles are formed and do not exhibit polarity outwards. Furthermore, the model considers the contribution of dislocation climb mechanism caused by vacancy diffusion to creep. This contribution mainly consists of two parts: first, the effect of dislocation climb on dislocation density, which alters the multiplication rate of dislocation dipoles. The specific formula is as follows:

[0148] ① The effect of dislocation climb on dislocation density:

[0149]

[0150] The second contribution is the contribution of dislocation climb to the activation frequency of dislocation crossing obstacles. Dislocation slip is generally classified into two types: continuous slip and discontinuous slip (jerky glide). Continuous slip is typically used under high stress conditions, where dislocations can directly pass through obstacles; discontinuous slip usually occurs under low stress conditions, where dislocations stop upon encountering an obstacle and then bypass it due to thermal fluctuations or diffusion. This model adopts the second type, discontinuous slip, and assumes that dislocations reach the next obstacle at extremely high speeds when bypassing it. Therefore, the dislocation velocity mainly depends on the waiting time before the obstacle. The activation frequency of crossing obstacles consists of two main parts: one is the direct crossing of obstacles caused by thermal disturbances, and the other is the crossing of obstacles through dislocation climb, as shown in the following formula:

[0151] ② The contribution of dislocation climb to the activation frequency when dislocations overcome obstacles:

[0152]

[0153]

[0154] It is important to note that all dislocation hazards are uniformly and averaged as forest dislocation hazards, which have an obstructive effect on dislocations. That is, dislocation hazards caused by other types of dislocations, such as dislocation knots and dislocation tangles, are simplified and averaged as forest dislocation hazards. This is because it would be extremely difficult to take all types of hazards into account.

[0155] The following high-temperature creep experiment will further illustrate this embodiment:

[0156] Experimental materials:

[0157] 316 is a stainless steel grade. AISI 316 is the corresponding American designation, and SUS 316 is the corresponding Japanese designation. In my country, the unified numerical code is S31603, and the standard grade is 022Cr17Ni 12Mo2 (new standard), while the old grade was 00Cr17Ni 14Mo2. This indicates that it mainly contains Cr, Ni, and Mo, with the numbers representing the approximate percentages. The material properties and chemical composition of 316 steel at room temperature are shown below:

[0158] Table 1. Basic mechanical properties of 316 austenitic stainless steel at room temperature (σ) 0.2 For yield strength, σ u (where A is tensile strength, Z is elongation, and A is reduction of area)

[0159]

[0160] Table 2 Chemical composition (wt%) of 316 stainless steel material

[0161]

[0162] According to the national standard GB / T 228-2002, uniaxial tensile tests were conducted on 316 stainless steel at 296K, 723K, 773K, 823K, 873K, 923K, and 973K to obtain the uniaxial tensile stress-strain curves and mechanical property indicators of 316 stainless steel. According to the national standard GB / T 2039-2012, high-temperature creep tests were conducted on 316 stainless steel at 773K-973K. Different stresses were selected at each temperature for high-temperature creep tests to obtain the creep-time curves and creep deformation characteristics of 316 stainless steel under different high-temperature stress levels.

[0163] Experimental principle:

[0164] High-temperature uniaxial tensile and creep tests were conducted on an RPL 100 electronic creep fatigue testing machine. During the experiment, the load sensor and displacement sensor on the testing machine converted the sensed displacement signals into electrical signals and sent them to the controller. After amplification and analog-to-digital conversion, the signals were sent to the computer, and the processed data were synchronously displayed on the screen to form a load-displacement curve. The experimental data can be stored and printed.

[0165] Experimental steps:

[0166] The high-temperature uniaxial tensile and creep tests on steel are conducted using the following procedures:

[0167] 1. Sample preparation: Number the specimens, measure the diameter of the three cross sections at both ends and the middle of the gauge length with vernier calipers, and take the average value once in two mutually perpendicular directions at each cross section. Calculate the original cross-sectional area S0 of the specimen using the smallest of the three average values, and set and record the experimental conditions and parameters.

[0168] 2. Heating: Install the specimen and extensometer. Heate the specimen in sections (usually three sections) in the experimental furnace at a heating rate of 20-40K / min to the specified temperature. After reaching the specified temperature, hold for 5 minutes. Collect the temperature data of the specimen while heating. During the heating process, apply a small constant force to both ends of the specimen to ensure that the pressure effect generated during the thermal expansion of the metal can be eliminated. After the holding period is completed, the load value changes slowly. At this time, zero the extensometer and hold for a certain period of time to prepare for loading.

[0169] 3. Loading: ① Single tension: After heat preservation, set the loading method according to the strain rate, setting the loading rate to 1 mm / min. Simultaneously collect load, displacement, and temperature data. ② Creep: After heat preservation, set the loading method according to stress loading, setting the working conditions through computer software. First, set the loading segment (not exceeding 10 minutes), then maintain the load for an isothermal creep test, collecting load, displacement, and temperature data simultaneously. The experiment ends when the specimen has fully creeped to the point of fracture.

[0170] 4. Disassembly and preparation: After the experiment, when the specimens have cooled down to room temperature naturally, the experimental specimens are disassembled and renumbered, the experimental characteristics are recorded, and the experimental data are processed.

[0171] Experimental results:

[0172] ① High-temperature uniaxial tensile test of stainless steel

[0173] First, uniaxial tensile tests were conducted at different temperatures (293K, 723K, 773K, 823K, 873K, 923K, 973K) using a CRIMS RPL 100 testing machine, with a strain rate of 1.3×10-3s-1. Figure 9 and Figure 10 The single tensile stress-strain curves and yield strength and tensile strength variation curves of 316 austenitic stainless steel at different temperatures are shown. The yield strength and tensile strength of 316 stainless steel at different temperatures are shown in Table 3.

[0174] Table 3. Yield strength and tensile strength of 316 austenitic stainless steel at different temperatures.

[0175]

[0176] Tests revealed that this type of 316 austenitic stainless steel has a room temperature yield strength (R0.2) of 278 MPa and an ultimate strength of approximately 608 MPa, while also possessing nearly 50% toughness. However, its yield strength and ultimate strength decrease significantly with increasing temperature. When the temperature rises to 723K-923K, a distinct sawtooth fluctuation appears on the stress-strain curve. This is mainly due to the significant dynamic strain aging (DSA) phenomenon of 316 austenitic stainless steel within a certain temperature and strain rate range. DSA is a phenomenon of strengthening caused by aging during plastic deformation. Essentially, within a certain temperature range, solute atoms in the alloy tend to cluster towards dislocations, forming gas clusters that pin dislocations, thus hindering dislocation movement. As the metal temperature gradually increases from low to high within the DSA temperature range, the clustering effect of solute atoms intensifies, and the process of forming pinned dislocation gas clusters becomes increasingly intense. However, when the temperature rises to a certain critical value, the activation energy reaches a certain level, and absorption sources that absorb solute atoms or vacancies appear inside the metal, thereby weakening the pinning effect of the gas cluster on dislocations. As the temperature continues to rise, this effect becomes increasingly stronger, eventually causing the sawtooth wave to disappear.

[0177] These sawtooth patterns can be broadly classified into three types: A, B, and C. Type A sawtooth patterns typically occur under low temperature or high strain rate conditions, exhibiting sawtooth fluctuations above the average stress value. As the temperature increases, Type A will transform into Type B sawtooth fluctuations, which fluctuate around the average value. Type C is the complete opposite of Type A fluctuations, with its fluctuations mainly below the average value.

[0178] Depend on Figure 11 and Figure 12 The results also revealed that the material exhibits strong toughness at room temperature, but its toughness significantly decreases when the material is within the DSA temperature range of 723K to 923K. Furthermore, when the experimental temperature exceeds 923K, the solute atom strengthening effect caused by the PLC (Portevin-LeChatlier) phenomenon disappears, and the material softens rapidly, resulting in improved toughness.

[0179] The main reason for the repeated changes in material toughness with temperature is the "blue brittleness" effect induced by DSA (Dielectric Strength Assay). Within the DSA temperature range, necking tends to concentrate within a very small range of the effective length of the specimen, resulting in reduced strain after necking and ultimately a decrease in total elongation and elongation rate. The "blue brittleness" effect is highly sensitive to strain rate; it can occur at high temperatures or at room temperature when the strain rate changes. Studies have also found that the fracture pattern changes when the "blue brittleness" effect occurs, exhibiting cleavage fracture morphology characteristics of a brittle fracture. However, it is important to note that the "blue brittleness" effect does not necessarily indicate brittle fracture; the fracture surface may still be ductile.

[0180] ②High-temperature creep test of 316 stainless steel

[0181] like Figure 13 and Figure 14 The results show that the creep of 316 austenitic stainless steel mainly consists of three stages. In the first stage, plastic deformation occurs due to dislocation slip within the crystal under external force. However, as deformation progresses, dislocation pile-up and pinning within the crystal hinder plastic deformation, causing the creep rate to decrease over time. When the internal hardening (dislocation pile-up, pinning, etc.) and softening (dislocation climb, etc.) mechanisms reach equilibrium, the material enters a steady-state creep stage. However, at high temperatures, the material softens due to recovery, disrupting the equilibrium and increasing the creep rate, ultimately leading to failure – the third stage of creep. Simultaneously, the experiment revealed that 316 austenitic stainless steel exhibits extreme temperature and stress sensitivity; that is, temperature and stress significantly influence the material's creep response, and its creep rate is positively correlated with both temperature and stress.

[0182] At the same temperature, as stress increases, the first stage of creep gradually disappears, the third stage of creep arrives prematurely, and the material rapidly fails. When stress is constant, the creep rate of the material increases with increasing temperature. The creep-time curve of 316 stainless steel exhibits sawtooth-shaped fluctuations, indicating that dynamic strain aging not only affects the material's tensile response but also influences its high-temperature creep response under specific temperature and rate conditions.

[0183] Simulation results are as follows Figure 15 and Figure 16 As shown

[0184]

[0185]

[0186] The embodiments described above are merely illustrative of specific implementations of the present invention, and while the descriptions are detailed, they should not be construed as limiting the scope of the present invention. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of the present invention, and these modifications and improvements all fall within the scope of protection of the present invention.

Claims

1. A finite element method for simulating crystal plasticity based on dislocation climb mechanism, characterized in that, include: In finite deformation, the overall deformation gradient tensor Decomposed into elastic components by multiplication and plastic parts : ； In the above formula, the plastic deformation gradient tensor The mapping of material points from a reference configuration to an intermediate configuration corresponds to the plastic deformation caused by dislocation motion in a crystal plasticity model, and is represented by the elastic deformation gradient tensor. It describes the mapping of material points from intermediate configurations to the current configuration, corresponding to lattice elastic stretching and rigid body rotation in the crystal plasticity model; intermediate configuration plastic velocity gradient Represented as: ； When dislocation slip is the only source of plastic deformation, It can be expressed as the shear strain rate on all slip systems of the crystal. Superposition: ； Where the unit vector and Describe the slip system in the intermediate configuration respectively The slip direction and slip surface normal of the edge dislocation; for the dislocation-dominated model, its Shear strain rate on slip system Represented as: ； in This is the average velocity of dislocation motion in the crystal. The above equation is rewritten in the following form: ； in It is a slip system Effective decomposition of shear stress, It is the strength of the solid solution. It is the length of the sliding Burgh's vector. It is the velocity of dislocation slip. It is the speed of dislocation climb. The physical meaning of this parameter is the obstacle characteristic size of the climb bypass mechanism. This parameter controls the strength of the climb-enhanced mobility of edge dislocations. It is the energy required to activate dislocation slip. It is the Boltzmann constant. It's temperature. and Fitting parameters reflecting the barrier energy barrier; Microstructure evolution from edge dislocation density With dipole dislocation density Characterization; when positive and negative dislocations move to a certain range When dislocations are inside, they annihilate, causing a decrease in dislocation density; at the same time, when in a state of... When the range is wide, dislocation dipoles will form and will not exhibit polarity outwards; The final evolution of the edge dislocation density is given by the following equation: ； in The mean free path of dislocation slip; the evolution of dislocation dipole density is given by the following equation: ； choose The stress tensor is solved as the starting point for calculation. After the stress tensor is solved and updated, the variables at the microstructure level are iteratively updated.

2. The finite element method for simulating crystal plasticity based on dislocation climb mechanism according to claim 1, characterized in that, Effective shear stress The calculation is as follows: ； in Sliding system The shear stress on the surface, the obstacle resistance that dislocation slip needs to overcome. The definition is as follows: ； in It is the shear modulus. Different slip systems and The mutual influence coefficient.

3. The finite element simulation method for crystal plasticity based on dislocation climb mechanism according to claim 2, characterized in that, The average velocity of dislocation motion in a crystal The calculation method is as follows: ； ； speed of dislocation climb The calculation method is as follows: ； in It is the self-diffusion coefficient of face-centered cubic geometry. It is the volume activated by dislocation climb. It is the activation energy for dislocation climb; When two dislocations can form a dislocation dipole It is the maximum distance between the two sliding planes; It is the minimum distance between the annihilation of two edge dislocations; and The calculations are as follows: ； ； in It is a fitting parameter; strain hardening passes through the mean free path of dislocations. To describe: ； ； in It is the average grain size. Different slip systems and The mutual influence coefficient.

Images