A Parallel Multiscale Method for Efficiently Calculating the Fracture Toughness of Solid Electrolytes

Abstract

A parallel multi-scale calculation method for effectively calculating the fracture toughness of solid electrolytes (GDCs) is disclosed. The method involves: calculating a millimeter-scale model using ABAQUS; transforming the millimeter-scale model into a nanometer-scale model using an equivalent stress intensity factor; and performing multi-scale calculations on the core region of the nanometer-scale model, with boundary conditions provided by the ABAQUS calculation results. This invention enables the study of the uniaxial tensile fracture behavior of GDCs with a central crack at different temperatures. Based on fracture mechanics theory, the macroscopic structure of the GDC is transformed into a microscopic intermediate transition model using the equivalent stress intensity factor method, and a detailed calculation process for analyzing the fracture toughness of the macroscopic structure of the GDC is provided through a multi-scale method. Furthermore, this multi-scale method investigates the fracture toughness of the macroscopic structure of the GDC under a force-electrochemical coupled field.

Description

Technical Field

[0001] This invention relates to a multi-scale calculation method, specifically a method for calculating the fracture toughness of materials by considering the simulation process of atoms to a continuous medium under the action of a force-electrochemical coupling field at high temperatures. Background Technology

[0002] Existing atomic-to-continuous-medium computational methods can compute fracture problems under uniform loading at the nanometer scale[1] and thermal diffusion problems[2]. However, there are no suitable computational models for fracture problems at the millimeter scale and fracture problems under complex boundary conditions.

[0003] [1]Zimmerman, Jonathan A., et al. "Calculation of stress in atomisticsimulation." Modelingandsimulation in materials science and engineering 12.4(2004):S319.

[0004] [2]Wagner,Gregory J.,et al."An atomistic-to-continuum coupling method for heat transfer in solids."Computer Methods inAppliedMechanics andEngineering 197.41-42(2008):3351-3365. Summary of the Invention

[0005] The purpose of this invention is to address the lack of suitable computational models for fracture problems at the millimeter level and under complex boundary conditions. This invention provides a multi-scale computational method from atoms to continuous media to calculate the fracture toughness of electrolytes subjected to force-electrochemical coupling fields at high temperatures.

[0006] To achieve the above objectives, the technical solution adopted by the present invention is as follows:

[0007] A parallel multi-scale calculation method for effectively calculating the fracture toughness of solid electrolytes is proposed, the method being as follows:

[0008] (1) Use ABAQUS to calculate millimeter-level models;

[0009] (2) Use the equivalent stress intensity factor to transform the millimeter-scale model into a nanometer-scale model;

[0010] (3) Multi-scale calculation of the core region of the nanoscale model, with boundary conditions provided by ABAQUS calculation results. Step (3) yields the stress-strain diagram of the core region obtained by multi-scale calculation. Based on this, it is determined how much loading is required to reach the critical state, and then the loading at this point is substituted into the formula to obtain the fracture toughness.

[0011] Furthermore, in step one, such as Figure 1 The model shown is a plane strain model, using a plane eight-node quadrilateral element, uniaxial tension, and calculated using ABAQUS finite element software.

[0012] Furthermore, step (2) specifically includes:

[0013] The formula for calculating the equivalent stress intensity factor is as follows:

[0014]

[0015] Among them, K I Here, is the equivalent stress intensity factor, F is a parameter related to loading and geometry, S is the magnitude of the uniaxial tensile load, and a is half the crack length.

[0016] For the above model, the formula for calculating F is:

[0017]

[0018] α=a / b (3)

[0019] α is the ratio of a to b, which has no practical significance, and b is half the model length in the crack direction;

[0020] During loading, the micrometer model is transformed into a micro-nano model by keeping the equivalent stress intensity factor constant.

[0021] Furthermore, step (3) specifically includes:

[0022] (1) The relationship between atomic quantities and continuous medium quantities is established by using localization functions as follows:

[0023] ρ * (x,t)=∑m α δ[xx α (t)] (4)

[0024] Where, ρ * (x,t) represents the relationship between atomic weights and quantities in a continuous medium, m α Let α be the mass of atom α, δ be the Dirac operator, and x be a node in the continuous medium. α (t) represents the coordinates of atom α at time t;

[0025] (2) By using the weighted residual method, the relationship between atomic mass and nodal quantity in continuous medium is established, and the nodal quantity is expressed in terms of atomic mass as follows:

[0026]

[0027] Where, ρ I Ω is the continuous medium quantity at node I, and Ω is the computational domain. 2 ρ represents the square norm. * These are the atomic quantities at the corresponding nodes, where I is the node I of the continuous medium mesh, and N... I Here, V is the finite element shape function at node I, V is the integration region, and J represents the nodes J and N of the continuous medium mesh, respectively. J It is the finite element shape function at node J, ρ J M is the continuous medium quantity at node J. IJ It is the formula corresponding to the square brackets, that is x α Let N be the coordinates of atom α. Iα The shape function value at the α coordinate of the atom. It is to sum over J and α respectively. For M IJ The inverse matrix,

[0028] Suppose a system of particles starts from a certain moment under ideal (linear or nonlinear) constraints and active forces. Then, for the possible accelerations that satisfy the constraints, establish the constraint function.

[0029]

[0030] Where, m i It is the mass of a point mass, r i It is the displacement of a particle. It is the second derivative of the displacement, F i It is the force acting on the atom;

[0031] Using the Gaussian minimum constraint principle and the Lagrange multiplier method, we have

[0032]

[0033] in, It is the actual force acting on the atoms. The atomic force is calculated from the potential function, λ. I It is a Lagrange multiplier. These are constraint equations, and the generalized form of the constraint equations is:

[0034]

[0035] in, Let α be the atomic weight of atom α. and Let A be the coordinates and velocity of atom α. I It is the continuous medium node quantity at node I;

[0036] Using the Gaussian minimum constraint method, we take the extreme values ​​of equation (7).

[0037]

[0038]

[0039] Where Φ is the potential function;

[0040] Rearranging equation (10), the atomic force can be expressed as two parts: the standard molecular dynamics potential function and the effect of the continuous medium on atoms in a representative volume.

[0041]

[0042] Force coupling constraint means that the total atomic momentum within the volume is consistent with the corresponding nodal momentum:

[0043]

[0044] in, It is the derivative of the momentum at node J;

[0045] A two-way information transfer process between atomic and continuous medium quantities was established. Considering the governing equations of the entire system, the system momentum in equation (5) was decomposed into two parts: the continuous medium region and the atomic region. Then, the time derivatives of both sides of the equation were taken simultaneously.

[0046]

[0047] Ω FE It is a finite element region, N Jα It is the shape function value at the α atom, and ρ is the quantity of the continuous medium. It is the node velocity of a continuous medium;

[0048] The above establishes an algorithm for multi-scale model calculation.

[0049] The advantages of this invention compared to existing technologies are as follows: This invention can study the uniaxial tensile fracture behavior of GDC (Gadolinium doped ceria, where a coefficient before GDC indicates the proportion of gadolinium ions in the cation, e.g., 10GDC indicates that the gadolinium ions account for 10% of the cation) electrolytes with a central crack at different temperatures. Based on fracture mechanics theory, the macroscopic structure of the GDC electrolyte is transformed into a microscopic intermediate transition model using the equivalent stress intensity factor method, and a detailed calculation process for analyzing the fracture toughness of the macroscopic structure of GDC is provided through a multi-scale method. Furthermore, this multi-scale method investigates the fracture toughness of the macroscopic structure of GDC under a mechanochemical coupled field. Attached Figure Description

[0050] Figure 1 Diagram of non-stoichiometric oxygen vacancy calculation models for 10GDC and 20GDC;

[0051] Figure 2 This is a flowchart of the method of the present invention;

[0052] Figure 3 This is a diagram showing the distribution of non-stoichiometric oxygen vacancy concentrations near the crack tip at 10 GDC 800 °C.

[0053] Figure 4 This is a diagram showing the distribution of non-stoichiometric oxygen vacancy concentrations near the crack tip at 10 GDC and 900 °C.

[0054] Figure 5 This is a diagram showing the distribution of non-stoichiometric oxygen vacancy concentrations near the crack tip at 20 GDC and 800 °C.

[0055] Figure 6 This is a diagram showing the distribution of non-stoichiometric oxygen vacancy concentrations near the crack tip at 20 GDC and 900 °C.

[0056] Figure 7 To calculate the stress-strain diagram of 10GDC under a loading of 0MPa and a temperature of 800℃ using a multi-scale method;

[0057] Figure 8 To calculate the stress-strain diagram of 10GDC under a loading of 0.4MPa and a temperature of 800℃ using a multi-scale method;

[0058] Figure 9 To calculate the stress-strain diagram of 10GDC under a loading of 0MPa and a temperature of 900℃ using a multi-scale method;

[0059] Figure 10 To calculate the stress-strain diagram of 10GDC under a loading of 0.4MPa and a temperature of 900℃ using a multi-scale method;

[0060] Figure 11To calculate the stress-strain diagram of 20GDC under a loading of 0MPa and a temperature of 800℃ using a multi-scale method;

[0061] Figure 12 To calculate the stress-strain diagram of 20GDC under a loading of 0.4MPa and a temperature of 800℃ using a multi-scale method;

[0062] Figure 13 To calculate the stress-strain diagram of 20GDC under a loading of 0MPa and a temperature of 900℃ using a multi-scale method;

[0063] Figure 14 To calculate the stress-strain diagram of 20GDC under a load of 0.4MPa and a temperature of 900℃ using a multi-scale method. Detailed Implementation

[0064] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to specific embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and do not limit the scope of protection of the invention.

[0065] Multiscale computational methods, from atom to continuum, combine atomic-level simulations with continuum-level simulations to study the properties and behavior of matter. This approach uses different computational techniques and models to describe behavior at both the atomic and continuum levels, obtaining accurate and comprehensive information at different scales. Atomic-level simulations, typically employing molecular dynamics (MD) or quantum mechanics, can simulate microscopic processes such as interactions between atoms, atomic motion, and chemical reactions. These simulations provide atomic-level detail, but their applicability is often limited by computational resources and time constraints. Continuum-level simulations, using continuum mechanics or the finite element method, treat matter as a continuous process, establishing mathematical models based on macroscopic physical laws. These simulations can describe large-scale material behavior, such as fluid flow and solid deformation, but cannot provide atomic-level detail. Multiscale computational methods combine these two levels of simulation, establishing coupling between the atomic and continuum levels to facilitate information exchange and transfer. This can be achieved by embedding atomic-level simulation results into continuum models or by using coarse-grained models to transfer information from the continuum level to the atomic-level simulations.

[0066] The principle of this invention is as follows: by deriving the distribution of non-stoichiometric oxygen vacancy concentration in the electrolyte, the non-stoichiometric oxygen vacancy concentration at the crack tip is obtained. Using the equivalent stress intensity factor method, the millimeter-level model is transformed into a nanometer-level model while maintaining the stress intensity factor unchanged. A multi-scale calculation method is used for the core region of the nanometer-level model, while other regions are calculated using the finite element software ABAQUS. Displacement is applied to each node of the boundary nodes in the multi-scale calculation region, and the magnitude of the displacement is calculated by the finite element software ABAQUS.

[0067] Example 1:

[0068] Calculate the fracture toughness of a 10GDC electrolyte at high temperature under low oxygen partial pressure on one side.

[0069] like Figure 1 As shown, this is a model diagram for calculating the non-stoichiometric oxygen vacancy concentration distribution in the 10GDC model. The plane strain model has a length and width of 1 mm, a central crack length of 0.04 mm, and a low oxygen partial pressure concentration on the left side of the model of logP(O2) = 18. Figure 2 As shown in the figure, this is a flowchart of the calculation of fracture toughness using a multi-scale method. Figure 2 (A) is a millimeter-scale model. Figure 2 (A) Model and Figure 1 The model is the same. Figure 2 (B) is a nanoscale model, which uses a method that keeps the stress intensity factor constant. Figure 2 (A) transformed into Figure 2 (B), Figure 2 (C) represents the mesh near the crack, and the multi-scale boundary conditions are obtained from the model in the ABAQUS computation diagram (B). Figure 2 (D) is a schematic diagram of multi-scale computation. Using Figure 1 After calculating the non-stoichiometric oxygen vacancy concentration distribution at the crack tip using the model, then... Figure 2 The multi-scale model is used to calculate the region near the crack, and the fracture toughness can be obtained by the final calculation.

[0070] (1) Use the displacement of each calculation node in the ABAQUS model.

[0071] (2) Calculate the non-stoichiometric oxygen vacancy concentration at the nodes by means of node displacement.

[0072] (3) The fracture toughness of the model is obtained by calculating the region near the crack using a multi-scale model.

[0073] right Figure 1The model was loaded with a load range of 0-0.4 MPa and a load increment of 0.05 MPa. The electrochemical boundary conditions were an oxygen partial pressure of logP(O2) = 18 and temperatures of 800℃ and 900℃. The calculation results are as follows: Figure 3 and 4 As shown, Figure 3 The results are calculated from the non-stoichiometric oxygen vacancy concentration distribution at the tip of a 10GDC crack at 800℃. Figure 4 This is the calculated result of the non-stoichiometric oxygen vacancy concentration distribution at the crack tip of a 10GDC crack at 900℃. Loads of 0 MPa and 0.4 MPa, and temperatures of 800℃ and 900℃ were selected. A multi-scale model was used to calculate the stress-strain curves in the crack tip core region at these temperatures. The calculation flowchart is shown below. Figure 2 As shown. The calculation results are as follows. Figures 7-10 As shown, Figures 7-10 The figures represent stress-strain diagrams obtained from multi-scale calculations at loading levels of 0 MPa at 800℃, 0.4 MPa at 800℃, 0 MPa at 900℃, and 0.4 MPa at 900℃, respectively. Each figure also shows the atomic configuration at the initial stage of the simulation, the maximum stress value, and the end of the simulation. The fracture toughness is calculated by substituting the uniaxial tensile loading value at the maximum stress value into the formula. The fracture toughness (K0) of 10GDC at 800℃ and 900℃ is also shown. IC The fracture toughness is shown in the table below:

[0074]

[0075] Example 2:

[0076] Calculate the fracture toughness of a 20GDC electrolyte at high temperature under low oxygen partial pressure on one side.

[0077] like Figure 1 As shown, this is a model diagram for calculating the non-stoichiometric oxygen vacancy concentration distribution in the 20GDC model. The plane strain model has a length and width of 1 mm, a central crack length of 0.04 mm, and a low oxygen partial pressure concentration on the left side of the model of logP(O2) = 18. Figure 2 As shown in the figure, this is a flowchart of the calculation of fracture toughness using a multi-scale method. Figure 2 (A) is a millimeter-scale model. Figure 2 (A) Model and Figure 1 The model is the same. Figure 2 (B) is a nanoscale model, which uses a method that keeps the stress intensity factor constant. Figure 2 (A) transformed into Figure 2 (B), Figure 2 (C) represents the mesh near the crack, and the multi-scale boundary conditions are obtained from the model in the ABAQUS computation diagram (B). Figure 2 (D) is a schematic diagram of multi-scale computation. Using Figure 1 After calculating the non-stoichiometric oxygen vacancy concentration distribution at the crack tip using the model, then... Figure 2 The multi-scale model is used to calculate the region near the crack, and the fracture toughness can be obtained by the final calculation.

[0078] (1) Use the displacement of each calculation node in the ABAQUS model.

[0079] (2) Calculate the non-stoichiometric oxygen vacancy concentration at the nodes by means of node displacement.

[0080] (3) The fracture toughness of the model is obtained by calculating the region near the crack using a multi-scale model.

[0081] right Figure 1 The model was used for recording, with a load range of 0-0.4 MPa and a load increment of 0.05 MPa. The electrochemical boundary conditions were an oxygen partial pressure of logP(O2) = 18 and temperatures of 800℃ and 900℃. The calculation results are as follows: Figure 5 and 6 As shown, Figure 5 The results are calculated from the non-stoichiometric oxygen vacancy concentration distribution at the 20GDC crack tip at 800℃. Figure 6 This is the calculated result of the non-stoichiometric oxygen vacancy concentration distribution at the crack tip of a 20GDC crack at 900℃. Loads of 0MPa and 0.4MPa, and temperatures of 800℃ and 900℃ were selected. A multi-scale model was used to calculate the stress-strain curves in the crack tip core region at these temperatures. The calculation flowchart is shown below. Figure 2 As shown. The calculation results are as follows. Figures 11-14 As shown, Figures 11-14 The figures represent the stress-strain diagrams obtained from multi-scale calculations at loading levels of 0 MPa at 800℃, 0.4 MPa at 800℃, 0 MPa at 900℃, and 0.4 MPa at 900℃, respectively. Each figure also shows the atomic configuration at the initial stage of the simulation, the maximum stress value, and the end of the simulation. The fracture toughness is calculated by substituting the uniaxial tensile loading value at the maximum stress value into the formula. The fracture toughness (K0) of 20GDC at 800℃ and 900℃ is also shown. IC The fracture toughness is shown in the table below:

[0082]

Claims

1. A parallel multi-scale calculation method for effectively calculating the fracture toughness of solid electrolytes, characterized in that: The method is as follows: Step 1: Use ABAQUS to calculate the millimeter-level model; Step Two: Transform the millimeter-scale model into a nanometer-scale model using the equivalent stress intensity factor; Step Two specifically involves: The formula for calculating the equivalent stress intensity factor is as follows: （1） in, Here, is the equivalent stress intensity factor, F is a parameter related to loading and geometry, S is the magnitude of the uniaxial tensile load, and a is half the crack length. For the above model, the formula for calculating F is: （2） （3） This is the ratio of a to b, which has no practical significance. b is half the model length in the crack direction. During loading, the micrometer model is transformed into a micro-nano model by keeping the equivalent stress intensity factor constant. Step 3: Calculate the core region of the nanoscale model at multiple scales, with boundary conditions provided by ABAQUS calculations; Step 3 specifically involves: (1) The relationship between atomic quantities and continuous medium quantities is established through localization functions as follows: （4） in, This relates to the relationship between atomic weights and quantities in a continuous medium. For atoms quality Here, x is a Dirac operator, representing a node in the continuous medium. for The coordinates of the atom at time t; (2) By using the weighted residual method, the relationship between atomic mass and nodal quantity in continuous medium is established, and the nodal quantity is expressed in terms of atomic mass as follows: （5） in, It is the continuous medium quantity at node I. It is the computational region. Denotes the square norm. This represents the atomic quantity at the corresponding node, where I is the node I of the continuous medium mesh. Here, V is the finite element shape function at node I, V is the integration region, and J represents the nodes J and N of the continuous medium mesh, respectively. J It is the finite element shape function at node J. It is the continuous medium quantity at node J. It is the formula corresponding to the square brackets, that is , For atoms coordinates For atoms Shape function values ​​at coordinates It is for J and respectively Summation, for The inverse matrix, ; Suppose a system of particles starts from a certain moment under ideal constraints and the action of active forces. Then, for the possible accelerations that satisfy the constraints, establish the constraint function. （6） Where, m i It is the mass of a point mass, r i It is the displacement of a particle. It is the second derivative of the displacement, F i It is the force acting on the atom; Using the Gaussian minimum constraint principle and the Lagrange multiplier method, we have （7） in, It is the actual force acting on the atoms. The forces acting on atoms are calculated from the potential function. It is a Lagrange multiplier. These are constraint equations, and the generalized form of the constraint equations is: （8） in, For atoms atomic weight and Atoms coordinates and velocity, It is the continuous medium node quantity at node I; Using the Gaussian minimum constraint method, the equations Take the extreme value, （9） （10） in, It is a potential function; Simplify the equations Atomic force can be expressed as two parts: the standard molecular dynamics potential function and the effect of the continuous medium on atoms in a representative volume. （11） Force coupling constraint means that the total atomic momentum within the volume is consistent with the corresponding nodal momentum: （12） in, It is the derivative of the momentum at node J; A two-way information transfer process between atomic quantities and continuous dielectric quantities was established. Considering the governing equations of the entire system, the equations were... The system momentum is decomposed into two parts: the continuous medium region and the atomic region. Then, the derivatives of both sides of the equation with respect to time are taken simultaneously: （13） It is a finite element region, yes The shape function values ​​at the atoms, It is a continuous medium quantity. It is the node velocity of a continuous medium; The above establishes an algorithm for multi-scale model calculation.

2. The parallel multi-scale calculation method for effectively calculating the fracture toughness of solid electrolytes according to claim 1, characterized in that: In step one, a planar eight-node quadrilateral element was used, with uniaxial tension, and the result was calculated using ABAQUS finite element software.

Images