A method for constructing a multi-scale plastic constitutive model of a metal material under multi-axial loading

Abstract

The present application relates to the technical fields of multi-scale plastic constitutive model construction, and especially to a kind of multi-scale plastic constitutive model construction method of metal material under multi-axial loading, and its technical scheme includes: representative volume element model and crystal plastic single crystal constitutive are constructed;Stress triaxiality, Lode angle parameter and non-proportionality are introduced, and its association with critical decomposition shear stress of slip system is established;Macro-micro two-way information transmission mechanism is established, the feedback of stress path to slip system activation and the correction of micro evolution to macro response are realized;Acceleration calculation module is constructed using order reduction modeling strategy;Path characteristic vector internal variable is introduced to describe loading path memory, and numerical implementation is completed combined with implicit integration algorithm;Through two-way coupling and multi-axial state system representation, the prediction accuracy of non-proportional loading is significantly improved, and the calculation efficiency is greatly improved, which is suitable for accurate prediction of metal plastic behavior under complex multi-axial loading.

Description

Technical Field

[0001] This invention relates to the field of multi-scale plastic constitutive model construction technology, and in particular to a method for constructing a multi-scale plastic constitutive model of metallic materials under multi-axis loading. Background Technology

[0002] Metallic materials often endure complex multiaxial loading conditions during engineering service, such as combined tension and torsion, biaxial tension / compression, and even triaxial stress states. Their plastic deformation behavior exhibits significant path dependence and anisotropy. Accurately describing the plastic response of metallic materials under multiaxial loading is of great engineering significance for structural integrity assessment, forming process optimization, and reliability design. To address this need, scholars both domestically and internationally have developed various types of plastic constitutive models from different scales.

[0003] At the macroscopic phenomenological level, researchers have established numerous elastoplastic constitutive models based on the continuum mechanics framework, describing the hardening behavior of materials by introducing isotropic hardening, kinematic hardening, and hybrid hardening rules. For multiaxial non-proportional loading conditions, some macroscopic models characterize the non-proportional additional hardening effect by introducing additional nonlinear terms into the back stress evolution equation and employ implicit integration methods for numerical implementation. However, the constitutive parameters of macroscopic phenomenological models lack clear physical meaning, and their parameter calibration is highly dependent on macroscopic experimental data under specific loading conditions. When the loading path changes, the predictive ability of the model often decreases significantly. Furthermore, such models struggle to reflect the microstructural evolution during the plastic deformation process of materials, such as dislocation density changes, texture evolution, and grain rotation, limiting their ability to provide a deeper description of the material response mechanism under complex multiaxial loading.

[0004] At the microscale, crystal plasticity theory provides a robust theoretical framework for describing the anisotropy and texture evolution of metallic materials. Multiscale simulation techniques based on the crystal plasticity finite element method, which establish single-crystal-scale constitutive models on representative volume element models and then obtain macroscopic responses through homogenization, have become an important approach to connecting the microstructure and macroscopic mechanical properties of materials. For example, Chinese patent CN119252402B discloses a method and system for constructing a multiscale constitutive model of the elastic viscoplasticity of polycrystalline materials. This method establishes a single-crystal-scale constitutive model based on a representative volume element model of the polycrystalline material, obtaining the relationship between stress rate and strain rate of a single grain, and then obtaining the relationship between stress rate and strain rate at the RVE scale. However, this type of method employs a unidirectional information transmission path "bottom-up" (grain scale → RVE scale → macroscopic scale). The stress state and loading path information at the macroscopic scale cannot be effectively fed back to the microscopic scale to influence the activation criteria and evolution of the slip system, resulting in insufficient prediction accuracy of the model under complex multiaxial loading paths.

[0005] To characterize multiaxial loading conditions, existing research has attempted to introduce multiaxial state parameters into the crystal plasticity framework. For example, some scholars have incorporated stress triaxiality into crystal plasticity constitutive models, introducing its influence on yield and rheological criteria in slip systems, thus elucidating the causes of tensile-compressive asymmetry in materials at both macroscopic and mesoscopic scales. However, most existing studies only involve the single parameter of stress triaxiality, failing to systematically incorporate multiaxial state descriptors such as Lode angle parameters and nonproportionality into the crystal plasticity constitutive framework, resulting in an insufficient characterization of multiaxial loading effects. Under nonproportional loading, materials exhibit additional strengthening effects, and their plastic strain is not only related to the current stress state but also closely related to the loading path. Existing models handle the loading path history effect in a relatively simple way, often using the accumulation of scalar internal variables, making it difficult to accurately capture the transient response and subsequent hardening evolution after path changes.

[0006] In terms of computational efficiency, multi-scale simulations based on crystal plasticity finite element methods require iterative calculations within large-scale incremental steps, resulting in high computational costs and severely limiting the application of this method in practical engineering. Although existing studies have attempted to improve computational efficiency by employing data-driven methods or reduced-order modeling strategies, these methods often sacrifice the interpretability of microscopic physical mechanisms and lack effective integration with bidirectional multi-scale information transfer frameworks.

[0007] In summary, existing methods for constructing multi-scale plastic constitutive models still have significant shortcomings in terms of scale coupling mechanisms, characterization of multiaxial loading effects, computational efficiency, and handling of loading path history, making it difficult to meet the practical engineering needs for accurate and efficient prediction of the plastic behavior of metallic materials under multiaxial loading. Therefore, there is an urgent need to develop a method for constructing multi-scale plastic constitutive models that can achieve bidirectional information transfer between macro and micro scales, systematically characterize multiaxial loading effects, and balance computational accuracy and efficiency.

[0008] Therefore, we propose a method for constructing a multi-scale plastic constitutive model of metallic materials under multi-axis loading to address the existing problems. Summary of the Invention

[0009] The purpose of this invention is to provide a method for constructing a multi-scale plastic constitutive model of metallic materials under multiaxial loading, so as to solve the problems mentioned in the background art.

[0010] To achieve the above objectives, the present invention provides the following technical solution: a method for constructing a multi-scale plastic constitutive model of metallic materials under multiaxial loading, comprising:

[0011] Step 1: Obtain the microstructure characteristic parameters of the target metallic material, including grain size distribution, grain orientation, phase composition and initial dislocation density. Based on the microstructure characteristic parameters, construct a representative volume element model of the polycrystalline material and apply periodic boundary conditions to the representative volume element model.

[0012] Step 2: Based on the representative volume element model, establish a single-crystal-scale constitutive model. The single-crystal-scale constitutive model uses crystal plasticity theory as a framework to define the evolution equations of decomposed shear stress, slip resistance and back stress of each slip system, and establishes the rate correlation relationship between the shear strain rate and decomposed shear stress of the slip system based on thermally activated dislocation dynamics.

[0013] Step 3: Establish a multiaxial loading state characterization module. The multiaxial loading state characterization module introduces three multiaxial state descriptors—stress triaxiality, Lode angle parameter, and non-proportionality—into the crystal plastic constitutive framework, and establishes a quantitative correlation function between the critical decomposition shear stress of each slip system and the three multiaxial state descriptors to characterize the influence of multiaxial loading conditions on the activation of slip systems.

[0014] Step 4: Establish a two-way information transfer mechanism between the macroscopic and microscopic scales. The two-way information transfer mechanism includes: transferring the stress tensor, strain tensor, and loading path information of the macroscopic scale down to the single-crystal scale constitutive model as the activation criterion of the slip system and the boundary conditions of the evolution equation; and transferring the texture evolution information, dislocation density evolution information, and slip system activation information of the microscopic scale up to the macroscopic scale to correct the constitutive response of the macroscopic scale.

[0015] Step 5: Construct a path-dependent multi-scale accelerated computing module. The accelerated computing module is based on a single-crystal-scale constitutive model and a two-way information transfer mechanism. It adopts a reduced-order modeling strategy and constructs a multi-scale response database under different loading paths offline. Combined with online interpolation or a lightweight neural network proxy model, it can quickly predict the plastic response of metallic materials under multi-axis loading.

[0016] Step Six: Introduce the loading path memory and evolution model. The loading path memory and evolution model uses the path feature vector as the internal variable and combines the plastic strain amplitude memory effect to describe the transient response and subsequent hardening evolution of the material after the loading path is changed. The multi-scale plastic constitutive model is numerically implemented through the implicit integral algorithm to complete the construction of the constitutive model.

[0017] Furthermore, the rate-dependent relationship between the shear strain rate and the decomposed shear stress of the slip system in step two is expressed in the form of: ,in, Let be the shear strain rate of the α-th slip system. The reference shear strain rate is given by F0, the activation free energy is given by k, the Boltzmann constant is given by T, and the absolute temperature is given by τ.α For the decomposed shear stress of the α-slip system, B α S is the back stress of the α-slip system. α Let τ be the slip resistance of the α-th slip system, τ0 be the reference stress, and p and q be material constants.

[0018] Furthermore, in step three, the stress triaxiality η is defined as the hydrostatic pressure σ. m Equivalent stress to Mises The ratio, The Lode angle parameter θ is defined as a function of the third invariant J3 and the second invariant J2 of the stress deviator. The non-proportionality Φ is defined as the degree of deviation between the current stress path and the proportional loading path. σ is the stress tensor. Let be the stress rate tensor.

[0019] Furthermore, the quantitative correlation function between the critical decomposed shear stress of each slip system and the three multiaxial state descriptors in step three is as follows: ,in, Let be the critical decomposed shear stress of the α-slip system in the current multiaxial state. f1(η) is the critical decomposed shear stress under the reference state, f2(θ) is the stress triaxiality influence function, f3(Φ) is the Lode angle parameter influence function, and f3(Φ) is the non-proportionality influence function. f1(η), f2(θ) and f3(Φ) are calibrated by uniaxial tensile compression test, pure shear test and non-proportional loading test, respectively.

[0020] Furthermore, the bidirectional information transmission mechanism in step four is implemented using an alternating iterative solution strategy: within each incremental step, the current macroscopic strain increment is first used as the boundary condition to drive the mesoscopic solution of the representative volume element model to obtain the mesoscopic stress field and texture evolution information; then the mesoscopic solution results are returned to the macroscopic scale through a homogenization method to update the macroscopic stress tensor and the uniform tangent stiffness matrix; the above process is repeated until the residuals at both the macroscopic and mesoscopic scales simultaneously satisfy the convergence criterion.

[0021] Furthermore, the homogenization method employs the equivalent stress-strain relationship based on strain energy, and the macroscopic equivalent stress... and macro equivalent change The stress σ of each grain in the representative volume element model is respectively determined by the stress σ. (r) and strain ε (r) Obtained by volume fraction weighted average , Where N is the total number of grains in the representative volume element model, and v (r) Let be the volume fraction of the r-th grain.

[0022] Furthermore, in step five, the reduced-order modeling strategy employs the intrinsic orthogonal decomposition method to reduce the dimensionality of the offline multi-scale response database, extract the dominant modes, and construct a reduced-order basis function space. Online interpolation or lightweight neural network surrogate models perform response prediction within the reduced-order basis function space, and the prediction results are physically consistent and corrected through a two-way information transmission mechanism.

[0023] Furthermore, in step six, the path feature vector P is defined as the stress path integral representation of the current loading path. Where ω is the path weight function, used to characterize different stress tensors σ and stress rate tensors. The influence weights on the evolution of the microstructure within the material, where t represents time; the path eigenvector P and the slip resistance S α and back stress B α The evolution equations are coupled, making the material hardening response dependent on the complete loading path history.

[0024] Furthermore, the implicit integration algorithm in step six adopts a radial backtracking mapping scheme, including: an elastic prediction step, which calculates the trial stress based on the current strain increment; a plastic correction step, which updates the stress tensor and internal variables by solving the plastic multiplier increment through Newton-Raphson iteration when the yield function is greater than zero; and a consistent tangent stiffness matrix step, which calculates the consistent tangent stiffness matrix after convergence for use in the next increment step.

[0025] Furthermore, it also includes model verification and calibration steps: after completing the constitutive model construction, multiaxial loading experiments are carried out by designing cruciform specimens or thin-walled circular tube specimens. The multiaxial loading experiments include proportional loading experiments and non-proportional loading experiments to obtain multiaxial stress-strain response data; the model prediction results are compared with the experimental data, and the material parameters in the model are calibrated by the parameter inversion method until the model prediction error meets the preset accuracy requirements.

[0026] Compared with the prior art, the beneficial effects of the present invention are:

[0027] 1. This invention establishes a bidirectional information transmission mechanism: on the one hand, the stress tensor, strain tensor, and loading path information at the macroscopic scale are transmitted downwards to the single-crystal-scale constitutive model as the activation criterion for the slip system and the boundary conditions for the evolution equation; on the other hand, the texture evolution information, dislocation density evolution information, and slip system activation information at the microscopic scale are transmitted upwards to the macroscopic scale to correct the constitutive response at the macroscopic scale. Through an alternating iterative solution strategy, the residuals at both the macroscopic and mesoscopic scales simultaneously satisfy the convergence criterion within each incremental step, overcoming the information loss problem in the traditional unidirectional homogenization method. This constructs a macroscopic-mesoscopic-microscopic bidirectional multi-scale information transmission mechanism, significantly improving the prediction accuracy under complex multi-axis loading conditions.

[0028] 2. This invention fully defines three multiaxial state descriptors: stress triaxiality, Lode angle parameter, and non-proportionality, and establishes a quantitative correlation function between these three state variables and the critical decomposed shear stress of each slip system. This correlation function, by decoupling and calibrating the influence of the three multiaxial state variables, can accurately characterize complex phenomena such as yield surface distortion, non-proportional additional strengthening, and texture evolution under multiaxial loading conditions. For stress states corresponding to different Lode angles at the same equivalent stress level, this invention can accurately distinguish the different slip system activation modes caused by them, with prediction accuracy significantly better than existing models that only consider stress triaxiality. The system introduces a multiaxial state characterization system, achieving accurate characterization of multiaxial loading effects.

[0029] 3. This invention employs intrinsic orthogonal decomposition to reduce the dimensionality of an offline multi-scale response database, extracting dominant modes and constructing a reduced-order basis function space. Response prediction is then performed using online interpolation or a lightweight neural network surrogate model, and the prediction results are physically consistent and corrected through a bidirectional information transfer mechanism. Since the reduced-order basis function space preserves the main physical modes of the crystal plastic constitutive model, the prediction results of the surrogate model still possess interpretability of the microscopic physical mechanism. Compared to traditional crystal plastic finite element simulation, this invention reduces computation time under a single multi-axis loading path while maintaining the same prediction accuracy, making the multi-scale plastic constitutive model feasible for engineering analysis of large and complex structures. A path-dependent multi-scale accelerated computation module has been developed, significantly improving computational efficiency while maintaining the interpretability of the microscopic physical mechanism.

[0030] 4. This invention defines a path feature vector, which completely preserves the historical information of the loading path in tensor form and directly couples it with the evolution equations of slip resistance and back stress. This path feature vector differs from the traditional scalar accumulation method for accumulated plastic strain; its tensor form carries the directional information of the stress path. When the loading path undergoes abrupt changes (such as a tension-torsion-tension path transition), the path feature vector can remember the path characteristics of the previous stage and influence the current hardening modulus in real time, thereby accurately reproducing the Bauschinger effect and cross-hardening effect. By introducing variables within the path feature vector, accurate prediction of transient response and hardening evolution after a change in the loading path is achieved.

[0031] 5. The rate-dependent relationship between shear strain rate and decomposed shear stress in this invention is established based on the thermally activated dislocation dynamics mechanism, and includes parameters with clear physical meaning such as activation free energy, reference stress, and material constants. Compared with the traditional power-law rate-dependent model, this thermally activated form can accurately describe the transrate plastic flow behavior from quasi-static to high strain rate, overcoming the problems of physical distortion under low stress and high strain rate conditions and excessively rapid saturation under high stress conditions in the power-law model. It establishes a rate-dependent constitutive framework based on thermally activated dislocation dynamics, thus broadening the strain rate applicability range of the model.

[0032] 6. The radial backtracking mapping scheme of this invention, through a three-step strategy of elastic prediction, plastic correction, and uniform tangent stiffness matrix, ensures the numerical stability of the constitutive integral under large incremental steps. For complex stress paths commonly encountered under multiaxial loading conditions, this algorithm can effectively avoid stress drift and convergence oscillation problems common in explicit integration methods. The explicit calculation of the uniform tangent stiffness matrix further guarantees the quadratic convergence speed when solving large-scale boundary value problems, enabling the constitutive model to be robustly embedded in finite element software for structural-level analysis.

[0033] 7. The model verification and calibration steps of this invention involve conducting proportional and non-proportional loading experiments using cruciform or thin-walled circular tube specimens to obtain multiaxial stress-strain response data, and then systematically calibrating the model parameters using parameter inversion methods. This verification and calibration process ensures that the model parameters have deterministic physical meaning and a repeatable experimental calibration path, avoiding the problems of parameter overfitting and insufficient generalization ability caused by pure phenomenological fitting, and guaranteeing the engineering applicability of the model under different material systems and different multiaxial loading conditions. Attached Figure Description

[0034] Figure 1 This is a schematic diagram illustrating the workflow of a method for constructing a multi-scale plastic constitutive model of metallic materials under multiaxial loading according to the present invention. Detailed Implementation

[0035] The technical solution of the present invention will be further described below with reference to the accompanying drawings and specific embodiments.

[0036] like Figure 1 As shown, a method for constructing a multi-scale plastic constitutive model of metallic materials under multi-axis loading includes the following steps: material preparation and microstructure characterization, construction of a representative volume element model, establishment of a single-crystal scale constitutive model, establishment of a multi-axis loading state characterization module, realization of a two-way information transmission mechanism, construction of a path-dependent multi-scale accelerated calculation module, realization of a loading path memory and evolution model and implicit integral algorithm, and model verification and parameter calibration.

[0037] Example 1: Material Preparation and Microstructure Characterization

[0038] In this embodiment, 316L austenitic stainless steel was selected as the target metal material. First, the microstructural characteristic parameters of the material were obtained.

[0039] (1) Sample preparation. Metallographic samples are cut from the material to be tested, and after mechanical grinding and polishing, electrolytic polishing or vibratory polishing is performed to remove the surface deformation layer and obtain a surface quality suitable for electron backscatter diffraction (EBSD) analysis.

[0040] (2) EBSD characterization. The crystal orientation and microstructure of the sample were characterized using a scanning electron microscope equipped with an EBSD detection system. The scanning area should cover a sufficient number of grains (ideally no less than 500 grains). The scanning step size was determined based on the average grain size, typically set to 1 / 10 to 1 / 5 of the average grain size. The acquired data included: Euler angles (or orientation matrix) of each grain, grain size distribution, grain boundary type and distribution, phase composition, and phase fraction.

[0041] (3) Extraction of microstructure characteristic parameters. Post-processing analysis of the EBSD data yielded the following microstructure characteristic parameters: average grain size... The parameters include the size distribution function, texture (orientation distribution function, ODF), initial dislocation density ρ0 (obtainable via X-ray diffraction or transmission electron microscopy), and phase volume fraction (if the material is multiphase). These parameters serve as input data for constructing a representative volume element model.

[0042] Example 2: Construction of a representative volume element model

[0043] After obtaining the microstructure characteristic parameters, the following steps are used to construct a representative volume element (RVE) model of the polycrystalline material.

[0044] (1) Voronoi mosaicking to generate a polycrystalline geometric model. Based on the average grain size and grain size distribution obtained in step one, a three-dimensional polycrystalline geometric model containing N grains is generated using the Voronoi mosaicking algorithm. The number of grains N should be selected to meet the statistical representativeness requirements of RVE. Generally, for equiaxed materials, N of 50 to 200 grains is sufficient to meet the representativeness requirements; for materials with significant texture or anisotropy, N should be appropriately increased to 200 to 500. The three-dimensional Voronoi model can be generated using the Neper open-source software, the Digimat commercial software, or a self-written Python / MATLAB program.

[0045] (2) Finite element mesh generation. The generated Voronoi polycrystalline geometric model is meshed using the finite element method. To ensure computational accuracy, each grain should contain at least 4-6 tetrahedral or hexahedral elements, with appropriate mesh refinement near grain boundaries to capture strain gradient effects. The mesh quality must meet the basic requirements of finite element analysis (e.g., Jacobian determinant greater than 0.3, aspect ratio less than 5, etc.). For an RVE model containing 512 grains, the total number of elements is typically between 100,000 and 500,000.

[0046] (3) Assigning grain orientations. The grain orientation data obtained from EBSD in step one is batch assigned to the corresponding grains in the RVE model using a Python script. If the number of grains in the EBSD data is inconsistent with the number of grains in the Voronoi model, random sampling or clustering methods can be used for matching. For materials with significant texture, it should be ensured that the assigned orientation distribution is consistent with the orientation distribution function (ODF) measured by EBSD.

[0047] (4) Apply periodic boundary conditions. Apply periodic boundary conditions to the opposite faces of the RVE model to ensure the continuity of displacement and stress on the RVE boundaries during deformation. The application of periodic boundary conditions is achieved by constraining the degrees of freedom of corresponding nodes on the opposite faces of the RVE model: ,in and These represent the displacements of corresponding nodes on the opposite surface. Let Δx be the macroscopic strain tensor. j The distance between the opposing surfaces is denoted as . These constraints can be automatically generated by a Python script and written into the input file of the finite element software.

[0048] Example 3: Establishment of a Constitutive Model at the Single Crystal Scale

[0049] Based on the RVE model, a single-crystal-scale constitutive model based on crystal plasticity theory is established.

[0050] (1) Kinematic description. The total deformation gradient F is decomposed into the elastic part F using the deformation gradient multiplication decomposition. e and plastic part F p F = F e ·F p Plastic deformation gradient rate Shear strain rate of each slip system The relationship between them is: , where n s The total number of slip systems (12 {111} in FCC crystal). <110> Slip systems (BCC crystals have 24 or 48 slip systems, while HCP crystals have 24 or 48 slip systems depending on the specific structure), s α and m α Let be the slip direction and the normal vector of the slip surface of the α-th slip system, respectively. Let it be the Schmid tensor.

[0051] (2) Decomposition of shear stress. Decomposition of shear stress τ on the α-slip system. α The Kirchhoff stress tensor τ is obtained by projecting it onto the Schmid tensor: .

[0052] (3) Rate-dependent slip shear strain rate evolution equation. The rate-dependent form based on thermally activated dislocation dynamics is adopted: ,in, Let τ be the shear strain rate of the α-th slip system. α For the decomposed shear stress of the α-slip system, B α S is the back stress of the α-slip system. α The slip resistance of the α-th slip system is given by the reference shear strain rate. Usually 10 is taken. -3 ~10 -4 s -1 F0 is the thermal activation free energy (approximately 0.4~0.6 eV for FCC copper and approximately 0.5~0.8 eV for 316L stainless steel), and k is the Boltzmann constant (1.38 × 10⁻⁶ eV). -23 J / K), T is the absolute temperature, τ0 is the reference stress, and p and q are material constants describing the shape of the activation energy barrier (usually taken as 0≤p≤1, 1≤q≤2).

[0053] (4) Evolution equation of slip resistance. Slip resistance S of the α-th slip system α The evolution takes into account the interactions between the various slip systems: ,in, Let β be the rate of change of the slip resistance of the α-th slip system, and β be the index of the other slip system. αβ Let be the slip hardening modulus matrix, where the diagonal terms (α=β) are the self-hardening moduli, and the off-diagonal terms (α≠β) are the latent hardening moduli. The self-hardening modulus typically follows the Voce-type hardening law: h αα = h0(1 - S α / S sat ) a Where h0 is the initial hardening modulus, S sat Let be the saturated slip resistance, and 'a' be the hardening exponent. The latent hardening modulus is taken as h. αβ = qh αα (α≠β), where q is the latent hardening coefficient (usually taken as 1.0~1.4). R s This is the dynamic response coefficient.

[0054] (5) Equation of back stress evolution. The back stress B of the α-slip system. α An Armstrong-Frederick type evolution equation is used to describe the kinematic hardening effect: ,in, Let c be the rate of change of the back stress in the α-slip system. b d is the back stress hardening coefficient. b This represents the dynamic recovery coefficient of the back stress.

[0055] Example 4: Establishment of a Multi-Axis Loading State Characterization Module

[0056] This embodiment details the construction process of the multi-axis loading state characterization module.

[0057] (1) Calculation of stress triaxiality η. In the current increment step, the macroscopic stress tensor obtained by RVE homogenization is... Calculate hydrostatic pressure Equivalent stress to Mises ,in Let I be the stress deviatoric tensor and I be the unit tensor. Then the stress triaxiality is: .

[0058] (2) Calculation of the Lode angle parameter θ. Calculate the second invariant of the stress deviatoric tensor. Given the third invariant J3 = det(s), the Lode angle parameters are: The range of θ is 0 ≤ θ ≤ π / 3, where θ = 0 corresponds to a pure tensile shear mode, θ = π / 6 corresponds to a pure shear mode, and θ = π / 3 corresponds to a pure compressive shear mode.

[0059] (3) Calculation of nonproportional coefficient Φ. Within the current increment step, the nonproportional coefficient Φ is calculated using the stress tensor σ and the stress rate tensor. Calculate the non-proportional degree: Where t represents time. In the actual numerical implementation, integration is achieved by accumulating over each increment step: , where ΔΦ acc = ||σ n :Δσ| - |σ n |·|Δσ||,ΔΦ norm = |σ n |·|Δσ|, where Δσ is the stress increment tensor within the current increment step, and σ n This represents the stress tensor at the start of the nth increment step. Φ=0 indicates fully proportional loading, Φ>0 indicates non-proportional loading, and the larger the value, the higher the degree of non-proportional loading.

[0060] (4) Multiaxial correction of critical decomposed shear stress. The critical decomposed shear stress of each slip system is corrected according to the following formula: ,in, Let be the critical decomposed shear stress of the α-slip system in the current multiaxial state. Let f1(η) be the critical decomposed shear stress of the α-slip system under the reference state (η=0, θ=π / 6, Φ=0). The specific forms of each influence function are determined through experimental calibration: f1(η) is calibrated by obtaining the yield stress at different η values ​​through uniaxial tensile tests (using round bar specimens with different notch radii) under different stress triaxialities, and fitting the result to obtain f1(η)=1+k1η+k2η. 2The values ​​of f2(θ) and f3(Φ) are calibration parameters. The yield stress at different Lode angles is obtained through pure tension, pure torsion, and pure compression experiments, and the result is f2(θ) = 1 + k3cos(3θ), where k3 is the calibration parameter. The yield stress of f3(Φ) is obtained by comparing the yield stress from non-proportional loading experiments (such as a tension-torsion non-proportional path) and proportional loading experiments, and the result is f3(Φ) = 1 + k4Φ, where k4 is the calibration parameter.

[0061] After the above calibration is completed, in the calculation of each incremental step, the critical decomposed shear stress of each slip system is updated in real time according to the current multiaxial state variables, thereby achieving an accurate characterization of the multiaxial loading effect.

[0062] Example 5: Implementation of a Two-Way Information Transmission Mechanism

[0063] This embodiment details the alternating iterative solution strategy for the bidirectional information transfer mechanism between macro and micro scales.

[0064] (1) Information transfer from macro to micro (downward transfer). At the beginning of each increment step, the macroscopic strain increment of the current increment step is provided by the macroscopic finite element solver. and the current macroscopic stress state The above information is passed to the RVE model as boundary conditions: macroscopic strain increment. The periodic boundary conditions are transformed into displacement increments at nodes on the RVE boundary; macroscopic stress state. Used to calculate the current multiaxial state variables (η, θ, Φ), and then correct the critical decomposed shear stress of each slip system through the correlation function. .

[0065] (2) Solving at the mesoscale. The macroscopic strain increment is used as the boundary condition to drive the RVE model to solve the crystal plasticity finite element problem at the mesoscale. Within each mesoscale increment step, constitutive integrals are performed on each grain to update the stress increment Δσ of each grain. (r) Strain increment Δε (r) Texture evolution (grain orientation update) and dislocation density evolution.

[0066] (3) Information transfer from micro to macro (upward transfer). After the microscopic solution converges, the stress, strain, texture, and dislocation density of each grain in the RVE model are transferred to the macroscopic scale through homogenization. The macroscopic equivalent stress and equivalent strain are obtained by weighted average of volume fractions: , Simultaneously, calculate the macroscopically consistent tangent stiffness matrix. This is for use in macroscopic finite element solvers. Among them: For macroscopic stress tensor, representing the average stress response of the entire RVE under the current loading at the end of the (n+1)th increment step; The macroscopic strain tensor represents the volumetric average strain of the entire RVE at the end of the (n+1)th increment step; r is the grain index within the RVE, N is the total number of grains in the RVE model, and v... (r) Let be the volume fraction of the r-th grain; For microscopic (grain-level) stress tensor, representing the actual stress of the r-th grain at the end of the (n+1)-th increment step, calculated by the crystal plastic constitutive integral; For microscopic (grain-level) strain tensor, representing the strain tensor of the r-th grain at the end of the (n+1)-th increment step.

[0067] (4) Alternating Iteration Convergence Criterion. The above macroscopic → microscopic → macroscopic transmission process constitutes a complete information loop. Since there may be inconsistencies between the macroscopic strain increment and the microscopic homogenization result, an alternating iteration strategy is required. Within each macroscopic increment step, steps (1) to (3) are repeated until the following two residuals simultaneously satisfy the convergence criterion: , Where the superscript (k) indicates the k-th alternation iteration, e σ and e ε The preset convergence tolerance (usually 10) -4 ~10 -6 ).

[0068] Example 6: Construction of a path-dependent multi-scale accelerated computing module

[0069] This embodiment details the implementation process of the order reduction modeling strategy.

[0070] (1) Construction of offline database. For the target metallic material, sampling is performed in a typical multi-axis loading path space to construct a multi-scale response database. The sampling of the loading path space covers the following dimensions: stress triaxiality η∈[−1,1], Lode angle parameter θ∈[0,π / 3], nonproportionality Φ∈[0,1], and equivalent plastic strain range. ∈[0,0.5]. Representative sampling points are selected in the above parameter space (e.g., using the Latin hypercube sampling method, the number of sampling points is typically 500~2000). For each sampling point, the corresponding macroscopic stress-strain response data are obtained through a complete RVE crystal plasticity finite element simulation, and a database is constructed. , where X i Y is the input parameter vector (containing load path parameters and current state variables). i To output the response vector (containing stress and uniform tangent stiffness matrix), N s This represents the total number of sampling points in the database.

[0071] (2) Intrinsic orthogonal decomposition dimensionality reduction. Snapshot the output response matrix in the database. Perform intrinsic orthogonal decomposition (POD). First, calculate the correlation matrix C=Y. T Y, then solve the eigenvalue problem Cϕ i =λ i ϕ i ϕ i Let λ be the i-th eigenvector. i Let be the i-th eigenvalue. Sort the eigenvalues ​​in descending order and select the top r dominant modes (that satisfy...). Construct a reduced basis function space Φ r =[ϕ1,ϕ2,...,ϕ r Any response vector can be approximated as Y ≈ Φ r a, where a is a reduced-order coordinate vector (r << N) s ).

[0072] (3) Construction of the online proxy model. In the reduced basis function space, construct the mapping relationship from the input parameter X to the reduced coordinate a. One of the following two methods can be used:

[0073] Radial basis function interpolation: using sampling points X in the database i For interpolation nodes, the corresponding reduced-order coordinates To achieve the interpolation objective, a radial basis function interpolation model is constructed. w i Let ϕ be the interpolation weight of the i-th radial basis function, and let ϕ() be the radial basis function whose value depends only on the distance between the input points. Common forms include Gaussian function, multiple quadratic function, and inverse multiple quadratic function.

[0074] Lightweight Neural Network: Construct a feedforward neural network with 2-3 hidden layers, using X as input and a as output, and train it using a database. The number of neurons in the hidden layers is typically 32-128, and the activation function is either ReLU or tanh.

[0075] (4) Physical Consistency Verification and Correction. The results predicted by the surrogate model are physically consistent through a two-way information transmission mechanism. Specifically, the macroscopic response predicted by the surrogate model is used as the initial value, and a complete RVE mesoscopic solution (rather than a complete iteration) is performed to verify the following physical constraints: symmetry of the stress tensor, satisfaction of the yield criterion, and the second law of thermodynamics (non-negativity of plastic power). If the prediction result deviates from the physical constraints by more than a preset threshold, a fine correction calculation based on RVE is triggered, and the correction result is fed back to the database for online updates of the surrogate model.

[0076] Example 7: Implementation of Loading Path Memory and Evolution Model and Implicit Integration Algorithm

[0077] This embodiment details the updating of the path feature vector, its coupling with the evolution of slip resistance and back stress, and the implicit integration algorithm of the constitutive model.

[0078] (I) Implementation of Loading Path Memory and Evolution Model

[0079] (1) Update of path eigenvector P. In each incremental step, based on the current stress tensor σ... n and stress rate tensor Calculate path weight function Then update the path feature vector:

[0080]

[0081] The specific form of the path weight function is designed as follows:

[0082]

[0083] Among them, P n+1 P is the path feature vector at the end of the current increment step. n σ is the path feature vector at the end of the previous increment step. n Let n be the stress tensor at the beginning of the nth increment step. Let σ be the stress rate tensor within the current increment step, Δt be the time increment of the current increment step, ω0 be the weight magnitude parameter, and σ be the stress rate tensor within the current increment step. ref For the reference stress state (usually the stress state at the most recent unloading or path change), σ w ω is the weight width parameter. This form ensures that ω is assigned a larger weight when the stress state changes significantly (such as path turning), thereby strengthening the path memory effect.

[0084] (2) Coupling of path eigenvectors and slip drag evolution. The evolution equation of slip drag incorporates a correction to the path eigenvectors:

[0085]

[0086] Where g(P) is the path history correction factor, specifically in the form of:

[0087] g(P) = 1 + g0·||P - P prop ||

[0088] Among them, P prop This represents the path feature vector corresponding to pure proportional loading, and g0 is the calibration parameter. When the loading path deviates from proportional loading (||P - P...), prop (|| > 0), the hardening modulus increases, thus accurately describing the non-proportional additional strengthening effect.

[0089] (3) Coupling of path eigenvectors and back stress evolution. The evolution equation of back stress also incorporates the influence of path eigenvectors:

[0090]

[0091] in, This is an additional term for back stress caused by path change, used to describe the transient response during abrupt path changes. p0 This is an additional factor for the back stress path. This represents the rate of change of the path feature vector.

[0092] (ii) Implicit Integration Algorithm

[0093] The constitutive model is fully implicitly integrated using a radial backtracking mapping scheme.

[0094] (1) Elastic prediction step. Assuming no plastic deformation occurs within the current increment step, calculate the trial stress:

[0095] τ trial = τ n + C e :Δε

[0096] Where, τ trial To test the Kirchhoff stress tensor, τ n Let C be the Kirchhoff stress tensor at the end of the nth increment step. e Let Δε be the elastic stiffness tensor and Δε be the strain increment. Simultaneously, calculate the trial decomposed shear stress. Define the yield function:

[0097]

[0098] Among them, f α Let be the yield function value of the α-th slip system. The back stress of the previous convergence step, The sliding resistance of the previous convergence step, This represents the critical decomposed shear stress under the current multiaxial state.

[0099] If f of all slip systems α If the value is ≤ 0, then this increment step is a purely elastic step, the tested stress is the actual stress, and the calculation of this increment step ends.

[0100] (2) Plasticity correction steps. If f exists α For slip systems with a slip value greater than 0, plastic correction is applied. A system of nonlinear equations is established that includes all activated slip systems:

[0101] , α∈A

[0102]

[0103]

[0104] Among them, R α For the plastic flow residual of the α-activated slip system, Δγ α Let be the shear strain increment of the α-th slip system. The shear strain rate function defined in Example 3 is the current decomposed shear stress τ. α Slip resistance S α and back stress B α The function, where A is the set of currently active glide systems. For the residual updated by slip resistance, ΔS α (Δγ) represents the increment of slip resistance. For the residual of the back stress update, ΔB α (Δγ) represents the back stress increment. The above equations are solved using a Newton-Raphson iterative method. In each iteration, the residual vector R and the Jacobian matrix are calculated. (where x = [Δγ]) α , S α B α ] T Solve the linear equation system J·δx=−R to update the solution vector until the norm of the residuals is less than the preset tolerance (usually taken as 10). -8 ~10 -10 ).

[0105] (3) Stress update and state variable update. After the plastic correction converges, the stress tensor is updated:

[0106]

[0107] Update slip resistance Back stress and path feature vector P n+1 Simultaneously update grain orientation (through rotation tensor) ).

[0108] (4) Calculation of the uniform tangent stiffness matrix. To ensure the quadratic convergence rate of the macroscopic finite element solution, the uniform tangent stiffness matrix is ​​calculated after convergence:

[0109]

[0110] in Let be the Schmid tensor of the α-th slip system, and H be the tangential modulus matrix of plastic flow. The uniform tangential stiffness matrix is ​​passed to the macroscopic finite element solver in tensor form.

[0111] Example 8: Model Validation and Parameter Calibration

[0112] This embodiment illustrates the process of verifying and calibrating the constitutive model through multi-axis loading experiments.

[0113] (1) Multiaxial loading experimental design. Select the appropriate specimen type according to the form of the target material:

[0114] For bars or tubes: use thin-walled circular tube specimens and achieve a multiaxial stress state through a combination of axial force, internal pressure, and torque. The ratio of the inner diameter to the outer diameter of the thin-walled circular tube is recommended to be greater than 0.9 to ensure the thin-wall assumption holds.

[0115] For sheet metal: cross-shaped specimens are used, and loads of different proportions are applied in two mutually perpendicular directions using a biaxial tensile testing machine.

[0116] (2) Loading path design. The following multi-axis loading paths were designed for model validation and calibration:

[0117] Proportional loading paths include uniaxial tension, uniaxial compression, pure torsion, and biaxial proportional tension (ε). 11 :ε 22 = 1:1) and biaxial non-uniform stretching (e.g., ε) 11 :ε 22 = 2:1).

[0118] Non-proportional loading paths: including tension-torsion non-proportional paths (tension before torsion, torsion before tension), and biaxial non-proportional paths (e.g., tension before torsion, torsion before tension). 11 Directional stretching and ε 22 Directional stretching (path with a 90° turn) and multi-segment complex paths.

[0119] Cyclic loading paths: including symmetric and asymmetric cyclic loading, used to verify the accuracy of the back stress evolution model.

[0120] (3) Data Acquisition. Experiments were conducted under each of the above loading paths, and load-displacement curves or stress-strain curves were recorded. For thin-walled circular tube specimens, axial force, torque, internal pressure, and corresponding axial strain, shear strain, and circumferential strain were recorded. For cruciform specimens, forces and displacements in both loading directions were recorded.

[0121] (4) Model parameter calibration. Based on experimental data, the unknown parameters in the model are calibrated using the parameter inversion method. Parameters to be calibrated include: single-crystal elastic constant (C 11 C 12 C 44 Initial critical decomposition shear stress τ c0 Thermal activation parameters (F0, τ0, p, q), hardening parameters (h0, S) sata, q), back stress parameters (c) b d b ), multi-axis influence function parameters (k1, k2, k3, k4) and path memory parameters (ω0, σ). w g0, c p0 ).

[0122] The objective function for parameter inversion is the weighted squared error between the experimental data and the model predictions:

[0123]

[0124] Where Θ is the parameter vector to be calibrated, and N exp N represents the number of experiments. data The number of data points for each experiment, w ij These are the weighting coefficients. The measured stress value is for the i-th experiment and the j-th data point. Given a parameter vector Θ, the constitutive model predicts the stress value for the i-th experiment and the j-th data point. An optimization algorithm (such as a genetic algorithm, particle swarm optimization, or Bayesian calibration method) is used to minimize the objective function.

[0125] (5) Model Validation. The calibrated model is used to predict the material response under uncalibrated loading paths, and the prediction results are compared with experimental data. The model prediction error is defined as:

[0126]

[0127] Where N val The total number of experiments used for validation. Let be the stress tensor predicted by the model in the i-th validation experiment. Let be the experimentally measured stress tensor in the i-th verification experiment. When the prediction error is less than the preset accuracy requirement (usually 5%~10%), the model is considered to have passed verification and can be applied to engineering prediction.

[0128] Example 9: Software Implementation and Engineering Application of Constitutive Models

[0129] (1) Software Implementation. The above constitutive model is written as a user material subroutine that can be called by mainstream finite element software (such as ABAQUS, ANSYS, etc.). For ABAQUS, the UMAT or VUMAT subroutine is written in FORTRAN; for ANSYS, the USERMAT subroutine is written in C++ or FORTRAN. The main functions of the subroutine include:

[0130] Receive the strain increment Δε and state variables from the finite element master program;

[0131] Execute the implicit integration algorithm described in Example 7 to update the stress and state variables;

[0132] Calculate and return the uniform tangent stiffness matrix C alg ;

[0133] Update and store state variables (including slip resistance, back stress, cumulative shear strain, grain orientation, path eigenvectors, etc. for each slip system).

[0134] (2) Overview of the calculation process. The entire construction and calculation process of the multi-scale plastic constitutive model is as follows:

[0135] Offline phase: Complete EBSD characterization of materials → Construct RVE model → Establish single-crystal constitutive model → Construct multi-axis state characterization module → Establish bidirectional information transmission mechanism → Construct offline database through RVE microscopic simulation → Perform POD dimensionality reduction to construct proxy model.

[0136] Online phase: The macro finite element solver provides strain increments → queries the surrogate model to obtain initial predictions → performs physical consistency verification through a two-way information transmission mechanism → triggers RVE fine correction when necessary → updates macro stress and consistent tangent stiffness matrices → returns to the macro solver.

[0137] (3) Engineering Application Scenarios. The multi-scale plastic constitutive model constructed in this invention can be applied to the following engineering scenarios:

[0138] Aero-engine turbine disks and blades: strength assessment and life prediction under complex multiaxial thermo-mechanical coupling service conditions;

[0139] Nuclear reactor pressure vessels and piping: Structural integrity assessment under complex stress conditions;

[0140] Automotive body and chassis structural components: Prediction of deformation behavior under collision and fatigue loads;

[0141] Metal forming processes: prediction of forming limits and process optimization under complex multiaxial stress conditions.

[0142] The above specific embodiments are merely several preferred embodiments of the present invention. Based on the technical solutions of the present invention and the relevant teachings of the above embodiments, those skilled in the art can make various alternative improvements and combinations to the above specific embodiments.

Claims

1. A method for constructing a multi-scale plastic constitutive model of metallic materials under multi-axis loading, characterized in that, include: Step 1: Obtain the microstructure characteristic parameters of the target metallic material, including grain size distribution, grain orientation, phase composition and initial dislocation density. Based on the microstructure characteristic parameters, construct a representative volume element model of the polycrystalline material and apply periodic boundary conditions to the representative volume element model. Step 2: Based on the representative volume element model, establish a single-crystal-scale constitutive model. The single-crystal-scale constitutive model uses crystal plasticity theory as a framework to define the evolution equations of decomposed shear stress, slip resistance and back stress of each slip system, and establishes the rate correlation relationship between the shear strain rate and decomposed shear stress of the slip system based on thermally activated dislocation dynamics. Step 3: Establish a multiaxial loading state characterization module. The multiaxial loading state characterization module introduces three multiaxial state descriptors—stress triaxiality, Lode angle parameter, and non-proportionality—into the crystal plastic constitutive framework, and establishes a quantitative correlation function between the critical decomposition shear stress of each slip system and the three multiaxial state descriptors to characterize the influence of multiaxial loading conditions on the activation of slip systems. Step 4: Establish a two-way information transfer mechanism between the macroscopic and microscopic scales. The two-way information transfer mechanism includes: transferring the stress tensor, strain tensor, and loading path information of the macroscopic scale down to the single-crystal scale constitutive model as the activation criterion of the slip system and the boundary conditions of the evolution equation; and transferring the texture evolution information, dislocation density evolution information, and slip system activation information of the microscopic scale up to the macroscopic scale to correct the constitutive response of the macroscopic scale. Step 5: Construct a path-dependent multi-scale accelerated computing module. The accelerated computing module is based on a single-crystal-scale constitutive model and a two-way information transfer mechanism. It adopts a reduced-order modeling strategy and constructs a multi-scale response database under different loading paths offline. Combined with online interpolation or a lightweight neural network proxy model, it can quickly predict the plastic response of metallic materials under multi-axis loading. Step Six: Introduce the loading path memory and evolution model. The loading path memory and evolution model uses the path feature vector as the internal variable and combines the plastic strain amplitude memory effect to describe the transient response and subsequent hardening evolution of the material after the loading path is changed. The multi-scale plastic constitutive model is numerically implemented through the implicit integral algorithm to complete the construction of the constitutive model.

2. The method for constructing a multi-scale plastic constitutive model of metallic materials under multi-axis loading according to claim 1, characterized in that: In step two, the rate-dependent relationship between the shear strain rate and the decomposed shear stress of the slip system is expressed in the form of: ,in, Let be the shear strain rate of the α-th slip system. The reference shear strain rate is given by F0, the activation free energy is given by k, the Boltzmann constant is given by T, and the absolute temperature is given by τ. α For the decomposed shear stress of the α-slip system, B α S is the back stress of the α-slip system. α Let τ be the slip resistance of the α-th slip system, τ0 be the reference stress, and p and q be material constants.

3. The method for constructing a multi-scale plastic constitutive model of metallic materials under multi-axis loading according to claim 1, characterized in that: In step three, the stress triaxiality η is defined as the hydrostatic pressure σ. m Equivalent stress to Mises The ratio, The Lode angle parameter θ is defined as a function of the third invariant J3 and the second invariant J2 of the stress deviator. ; Non-proportionality Φ is defined as the degree of deviation between the current stress path and the proportional loading path. σ is the stress tensor. Let be the stress rate tensor.

4. The method for constructing a multi-scale plastic constitutive model of metallic materials under multi-axis loading according to claim 1, characterized in that: The quantitative correlation function between the critical decomposed shear stress of each slip system and the three multiaxial state descriptors in step three is as follows: ,in, Let be the critical decomposed shear stress of the α-slip system in the current multiaxial state. f1(η) is the critical decomposed shear stress under the reference state, f2(θ) is the stress triaxiality influence function, f3(Φ) is the Lode angle parameter influence function, and f3(Φ) is the non-proportionality influence function. f1(η), f2(θ) and f3(Φ) are calibrated by uniaxial tensile compression test, pure shear test and non-proportional loading test, respectively.

5. The method for constructing a multi-scale plastic constitutive model of metallic materials under multi-axis loading according to claim 1, characterized in that, In step four, the bidirectional information transfer mechanism is implemented using an alternating iterative solution strategy: within each increment step, the current macroscopic strain increment is first used as the boundary condition to drive the mesoscopic solution of the representative volume element model to obtain the mesoscopic stress field and texture evolution information; then the mesoscopic solution results are returned to the macroscopic scale through a homogenization method to update the macroscopic stress tensor and the uniform tangent stiffness matrix; the above process is repeated until the residuals at both the macroscopic and mesoscopic scales simultaneously satisfy the convergence criterion.

6. The method for constructing a multi-scale plastic constitutive model of metallic materials under multi-axis loading according to claim 5, characterized in that: The homogenization method employs the equivalent stress-strain relationship based on strain energy, and the macroscopic equivalent stress. and macro equivalent change The stress σ of each grain in the representative volume element model is respectively determined by the stress σ. (r) and strain ε (r) Obtained by volume fraction weighted average , Where N is the total number of grains in the representative volume element model, and v (r) Let be the volume fraction of the r-th grain.

7. The method for constructing a multi-scale plastic constitutive model of metallic materials under multi-axis loading according to claim 1, characterized in that: In step five, the reduced-order modeling strategy uses the intrinsic orthogonal decomposition method to reduce the dimensionality of the offline multi-scale response database, extract the dominant modes, and construct a reduced-order basis function space. Online interpolation or a lightweight neural network surrogate model performs response prediction in the reduced-order basis function space, and the prediction results are physically consistent and corrected through a two-way information transmission mechanism.

8. The method for constructing a multi-scale plastic constitutive model of metallic materials under multi-axis loading according to claim 1, characterized in that: In step six, the path feature vector P is defined as the stress path integral representation of the current loading path. Where ω is the path weight function, used to characterize different stress tensors σ and stress rate tensors. The influence weights on the evolution of the microstructure within the material, where t represents time; the path eigenvector P and the slip resistance S α and back stress B α The evolution equations are coupled, making the material hardening response dependent on the complete loading path history.

9. The method for constructing a multi-scale plastic constitutive model of metallic materials under multi-axis loading according to claim 1, characterized in that, Step six, the implicit integration algorithm, employs a radial backtracking mapping scheme, which includes: an elastic prediction step, calculating the trial stress based on the current strain increment; a plastic correction step, when the yield function is greater than zero, solving the plastic multiplier increment through Newton-Raphson iteration to update the stress tensor and internal variables; and a consistent tangent stiffness matrix step, calculating the consistent tangent stiffness matrix after convergence for use in the next increment step.

10. The method for constructing a multi-scale plastic constitutive model of metallic materials under multi-axis loading according to claim 1, characterized in that, It also includes model verification and calibration steps: After completing the constitutive model construction, multiaxial loading experiments are carried out by designing cruciform specimens or thin-walled circular tube specimens. Multiaxial loading experiments include proportional loading experiments and non-proportional loading experiments to obtain multiaxial stress-strain response data; the model prediction results are compared with the experimental data, and the material parameters in the model are calibrated by the parameter inversion method until the model prediction error meets the preset accuracy requirements.

Images