Calculation method for adding irregular sub-domain without flux boundary condition in computational domain
By adding an irregular subdomain method with no flux boundary conditions within the computational domain, model building is simplified, computational efficiency and accuracy are improved, the simulation problem of complex substrates or impurity phases in the prior art is solved, and better prediction of actual physical phenomena is achieved.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2025-11-28
- Publication Date
- 2026-04-07
AI Technical Summary
When adding substrate or impurity phases in the computational domain, existing technologies are limited by complex and time-consuming structural mesh generation, resulting in low computational efficiency and inability to meet accuracy requirements. They cannot effectively simulate substrate or impurity phases with complex geometries in actual engineering fields.
A computational method is adopted that adds an irregular subdomain without flux boundary conditions within the computational domain. By drawing the shape image of the substrate/impurity phase, processing it into data, establishing computational equations, discretizing the difference, writing the phase field pre-program, and importing the initial conditions, the model building is simplified, complex meshing is avoided, and computational efficiency and accuracy are improved.
It enables better prediction of actual physical phenomena, improves computational efficiency and accuracy, and can better simulate electrochemical energy storage processes on complex substrates or impurity phases.
Smart Images
Figure CN121809130A_ABST
Abstract
Description
TECHNICAL FIELD
[0001] The application belongs to the technical field of electrochemical energy storage, and particularly relates to a calculation method of adding an irregular subdomain with a flux-free boundary condition in a calculation domain. BACKGROUND
[0002] Electrochemical energy storage technology is strategic and pioneering, and is a focus of the academic and industrial communities in maintaining national energy security, and is also a key supporting technology for energy revolution. A comprehensive understanding of the morphology evolution in electrodeposition is crucial for designing the next generation of batteries with metal lithium anodes. The application of theoretical calculation to the research of lithium metal negative electrodes helps to fundamentally understand and solve the above-mentioned key problems. Based on the progress of theoretical algorithms, modeling simulation and computer technology, materials, batteries, devices and battery packs in the field of lithium-ion batteries are gradually combined through large-scale data sharing, which will bring extensive and far-reaching changes and accelerate the development of the entire industry chain. Complex interfacial reactions and nucleation deposition phenomena show that various factors in lithium battery systems, with complex relationships and effects, jointly affect the working performance of the lithium negative electrode. Therefore, it is necessary to use new methods to achieve overall optimization. For example, the lithiation process greatly affects the structure and performance of the matrix material, which requires the identification of the structure of the matrix material. Further, using the structure prediction method, the lithium ion migration ability and the dendrite inhibition ability are used as key indicators to systematically study related problems, so as to realize the theoretical optimization of the matrix material in the future, and hope to realize the design and discovery of the next generation of high-energy-density energy systems guided by theory and driven by data.
[0003] In order to better match the model with the actual problem, researchers often need to add a base phase or an impurity phase (such as a separator in a lithium battery) in the calculation domain. At present, the method for adding these phases is usually to directly establish a volume grid based on a surface grid or a base / impurity phase of a regular subdomain that is easy to combine with the grid in the model. The biggest limitation of this method is that the geometry of the domain to be solved must be simple enough. However, in the actual engineering field, the geometry of the added base or impurity phase is very complex, which makes it necessary to perform complex and time-consuming structure grid division on the target region during calculation, resulting in low calculation efficiency and inability to meet the accuracy requirements. SUMMARY
[0004] The technical problem to be solved by the present application is to provide a calculation method of adding an irregular subdomain with a flux-free boundary condition in a calculation domain, which simplifies the model building, improves the operation efficiency and accuracy, and enables the model to better match the actual situation, thereby better predicting the actual physical phenomenon.
[0005] To solve the above technical problems, the technical scheme adopted by the present application is: a calculation method for adding an irregular subdomain with a flux-free boundary condition in a calculation domain, comprising the following steps: Step S1, determining an irregular subdomain of a base / impurity phase added in the calculation domain, and drawing a base / impurity phase shape image; Step S2, processing the base / impurity phase shape image obtained in step S1 into data of a required model size by using computer graphics processing; Step S3, establishing a calculation equation of the calculation method for coupling the irregular subdomain with the flux-free boundary condition added in the calculation domain; Step S4, discretely differentiating the calculation equation established in step S3 by using a finite difference method to obtain a discrete difference format equation; Step S5, applying the discrete difference format equation in step S4 to write a phase field preprocessor; Step S6, importing the data obtained in step S2 into the phase field preprocessor obtained in step S5; Step S7, importing the running result of the phase field preprocessor in step S6 as an initial condition of the base / impurity phase into a phase field model constructed for describing a dendrite growth process to obtain a calculation model coupled with the base / impurity phase with the irregular domain; Step S8, obtaining a numerical solution of the calculation model in step S7 by computer simulation.
[0006] The calculation method for adding an irregular subdomain with a flux-free boundary condition in a calculation domain, the specific process of establishing a calculation equation of the calculation method for coupling the irregular subdomain with the flux-free boundary condition added in the calculation domain in step S3 is as follows: Step S301, deriving an initial condition of a Neumann boundary condition formula; Step S302, deriving a boundary condition of the Neumann boundary condition formula of the Neumann boundary applied to the diffusion interface; Step S303, deriving a contact angle boundary condition formula of the Neumann boundary condition formula; Step S304, introducing an irregular subdomain parameter , and merging the contact angle boundary condition into the original control equation; Step S305, listing the calculation equation.
[0007] When deriving the initial condition of the Neumann boundary condition formula in step S301, an arbitrary function is taken as an example, and the Laplacian of the arbitrary function is multiplied by a subdomain parameter , the differential method of the identity product is used to obtain:
[0008] Thus, the calculation term proportional to is obtained; wherein, represents the gradient sign, is the gradient of the function in space, is the gradient of the function in space.
[0009] In the above-mentioned calculation method of adding the irregular sub-domain of the no-flux boundary condition in the calculation domain, the Neumann boundary condition formula derived in step S302 is applied to the boundary condition of the boundary on the diffusion interface, and the inward unit normal vector of the boundary (pointing to ) is given by , and the boundary condition of the Neumann boundary applied to the boundary on the diffusion interface is expressed by the formula:
[0010] wherein, is the normal direction.
[0011] In the above-mentioned calculation method of adding the irregular sub-domain of the no-flux boundary condition in the calculation domain, the specific process of the contact angle boundary condition formula of the Neumann boundary condition formula in step S303 is as follows: Step S3031, for the A-C and C-H equations in the phase field method, the free energy form is expressed as:
[0012] wherein, is the free energy function, is the phase field order parameter used to define different phases, is the double potential well free energy functional, is the gradient energy coefficient, is the entire calculation region; Step S3032, at the extreme value of the free energy function , the variation derivative of the total free energy is obtained, and thus:
[0013] Step S3033, multiply on both sides of the formula obtained in step S3032 to obtain:
[0014] The contact angle boundary condition formula is obtained as follows: .
[0015] The calculation method of adding irregular sub-domains with no flux boundary condition in the calculation domain, the introduction of the irregular sub-domain parameter in step S304 When the contact angle boundary condition is combined into the original control equation, the specific process is as follows: Step S3041, express the contact angle θ as follows on the whole calculation domain:
[0016] Wherein, the sub-domain parameter satisfies the diffuse reflection boundary, and satisfies , is the unit normal vector of the phase interface (pointing to the region of ); Step S3042, the following equation of the contact angle boundary condition is derived: .
[0017] The calculation method of adding irregular sub-domains with no flux boundary condition in the calculation domain, the specific process of listing the calculation equation in step S305 is as follows: Step S3051, the chemical potential driving the morphology evolution is defined as the variational derivative of the system total free energy:
[0018] Step S3052, apply this calculation method to the chemical potential, multiply it by the sub-domain parameter , and apply the product rule and the boundary condition equation to obtain:
[0019] For a conserved order parameter, its evolution is controlled by the Cahn-Hilliard equation, in which the change rate of the order parameter is equal to the divergence of its flux, which is proportional to the gradient of the chemical potential: ; Wherein, is time, is the mobility coefficient; Step S3053, the calculation method obtains: is the flux of the conserved order parameter; wherein, is the material flux; represents the flux perpendicular to the domain boundary; The calculation method formula of step S3054, the Cahn-Hilliard equation, is written as:
[0020] wherein, is a contact angle boundary condition applied at the three-phase boundary, .
[0021] The above calculation method of adding a no-flux boundary condition irregular subdomain in the calculation domain is for a closed system, .
[0022] Compared with the prior art, the present application has the following advantages: the present application is based on rewriting the original partial differential equation into a simple form, and automatically coupling the no-flux boundary condition on the irregular domain; the method of the present application better establishes the irregular subdomain such as an electrode or a separator in electrochemical energy storage simulation and modeling, does not need to perform complex and time-consuming structure grid division on the target region, avoids the problem that the finite element grid can only construct a volume grid based on a surface grid or easily combines the regular subdomain of the grid; and simplifies the model building, improves the operation efficiency and accuracy, and can better fit the actual situation, thereby better predicting the actual physical phenomenon.
[0023] The technical solutions of the present application will be further described in detail below with reference to the drawings and embodiments. BRIEF DESCRIPTION OF DRAWINGS
[0024] Figure 1 It is a method flow chart of the present application; Figure 2 It is a two-dimensional image of a lithium metal battery separator drawn in a specific embodiment of the present application; Figure 3 The specific embodiment of the present application shows the difference between the sharp interface and the diffusion interface of the separator phase before and after the calculation processing, Figure 3 a is a separator phase with a sharp interface, Figure 3 b is a separator phase with a diffusion interface; Figure 4 It is a result graph after modeling of the electrodeposition model in a specific embodiment of the present application. DETAILED DESCRIPTION
[0025] As Figure 1 shown, the calculation method of adding a no-flux boundary condition irregular subdomain in the calculation domain of the present application includes the following steps: Step S1, determining the irregular subdomain of the added base / impurity phase in the calculation domain, and drawing it into a base / impurity phase shape image; Step S2, using computer graphics processing to process the substrate / impurity phase shape image obtained in step S1 into data of a required model size; Step S3, establishing a calculation equation of a calculation method of coupling an irregular subdomain adding a flux-free boundary condition in a calculation domain; Step S4, using a finite difference method to discretely difference the calculation equation established in step S3 to obtain a discrete difference format equation; Step S5, applying the discrete difference format equation in step S4 to write a phase field preprogram; Step S6, importing the data obtained in step S2 into the phase field preprogram obtained in step S5; Step S7, importing the running result of the phase field preprogram in step S6 as an initial condition of the substrate / impurity phase into the phase field model constructed for describing the dendrite growth process to obtain a calculation model coupling the substrate / impurity phase with the irregular domain; Step S8, obtaining a numerical solution of the calculation model in step S7 through computer simulation.
[0026] In this embodiment, the specific process of establishing the calculation equation of the calculation method of coupling the irregular subdomain adding the flux-free boundary condition in the calculation domain in step S3 is as follows: Step S301, deducing an initial condition of a Neumann boundary condition formula; Step S302, deducing a boundary condition of the Neumann boundary condition formula of the Neumann boundary applied on the diffusion interface; Step S303, deducing a contact angle boundary condition formula of the Neumann boundary condition formula; Step S304, introducing an irregular subdomain parameter , and merging the contact angle boundary condition into the original control equation; Step S305, listing the calculation equation.
[0027] In this embodiment, when the initial condition of the Neumann boundary condition formula is deduced in step S301, an arbitrary function is taken as an example, the Laplacian of which is multiplied by the subdomain parameter , and the differential method of the product of the identity is used to obtain:
[0028] Thus, a calculation term proportional to is obtained; wherein, represents a gradient symbol, is the gradient of the function in space, For sub-domain parameters Gradient in space.
[0029] In this embodiment, when the Neumann boundary condition formula is derived in step S302, the Neumann boundary condition is applied to the boundary condition on the diffusion interface, the inward unit normal vector of the boundary (pointing to the ) is given by , and the boundary condition on the diffusion interface is expressed by the formula:
[0030] wherein, is the normal direction.
[0031] In this embodiment, when the contact angle boundary condition formula is derived in step S303, the specific process is as follows: Step S3031, for the A-C and C-H equations in the phase field method, the free energy form is expressed as:
[0032] wherein, is the free energy function, is the phase field order parameter used to define different phases, is the double potential well free energy functional, is the gradient energy coefficient, is the entire calculation area. Step S3032, at the extreme value of the free energy function , the variation derivative of the total free energy , thus:
[0033] Step S3033, multiply on both sides of the formula obtained in step S3032, to obtain:
[0034] The contact angle boundary condition formula is obtained as: .
[0035] In this embodiment, when the irregular sub-domain parameter is introduced in step S304, and the contact angle boundary condition is combined into the original control equation, the specific process is as follows: Step S3041, the contact angle θ is expressed as:
[0036] where the subdomain parameter satisfies the diffuse reflection boundary, and satisfies , is the unit normal vector of the phase interface (pointing towards the region ); Step S3042, the following equation of the contact angle boundary condition is derived: .
[0037] In this embodiment, the specific process of listing the calculation equation in step S305 is as follows: Step S3051, the chemical potential driving the evolution is defined by the variational derivative of the total free energy of the system:
[0038] Step S3052, the calculation method is applied to the chemical potential, by multiplying it by the subdomain parameter , and the product rule and the boundary condition equation are applied to obtain:
[0039] For the conserved order parameter, its evolution is controlled by the Cahn-Hilliard equation, in which the rate of change of the order parameter is equal to the divergence of its flux, which is proportional to the gradient of the chemical potential: ; where is time, is the mobility coefficient; Step S3053, the calculation method is obtained from: is the flux of the conserved order parameter; where is the material flux; represents the flux perpendicular to the domain boundary; Step S3054, the calculation method formula of the Cahn-Hilliard equation is written as:
[0040] where is the contact angle boundary condition imposed at the three-phase boundary, .
[0041] In this embodiment, for a closed system, .
[0042] In order to verify the technical effects generated by the present application, the computing method of the present application is applied to establish a phase field model of the separator in the lithium metal battery, the method of step S1 is used to determine the required lithium metal battery separator phase shape, which is drawn into a picture as shown in Figure 2 For convenience, it is drawn into a black and white image; after steps S1-S8, the difference between the sharp interface and the diffusion interface of the separator phase before and after the calculation processing is as shown in Figure 3 Figure 3 a is the separator phase with a sharp interface, Figure 3 b is the separator phase with a diffusion interface; as can be seen from the figure, after a period of calculation simulation, the sharp interface of the separator phase (the density of the isochrones at the interface does not change) is completely converted into a diffusion interface (the isochrones on both sides of the interface are sparse, and the middle is dense).
[0043] The simulation results of the phase field preprocessor of step S6 in step S7 are imported into the electrodeposition model based on the smooth interface method as the initial conditions of the lithium metal battery separator phase, and the electrodeposition model coupled with the lithium metal battery separator phase with a diffusion interface is obtained; Figure 4 is the result after simulation of the electrodeposition model obtained by step S7, for the convenience of observation, the separator phase part is left blank. As can be seen from the figure, with the elapse of time, the deposition phase gradually advances to the liquid phase area, when the deposition phase advances to the separator, it will completely bypass the deposition phase to advance, which is consistent with the expectation of the lithium metal battery separator phase with a diffusion interface The above is only a preferred embodiment of the present application, and does not limit the present application, any simple modification, change and equivalent structure change according to the technical essence of the present application to the above embodiment are still within the protection scope of the technical solution of the present application.
Claims
1. A computational method for adding an irregular subdomain without flux boundary conditions within a computational domain, characterized in that, The method includes the following steps: Step S1: Determine the irregular subdomains of the substrate / impurity phase added in the computational domain and plot them as substrate / impurity phase shape images; Step S2: Use computer graphics processing to process the substrate / impurity phase shape image obtained in step S1 into data of the required model size; Step S3: Establish the computational equations for a computational method that couples an irregular subdomain with no flux boundary conditions within the computational domain. Step S4: Discrete the computational equations established in step S3 using the finite difference method to obtain the discrete difference scheme equations. Step S5: Write the phase field preprocessing program using the discrete difference scheme equations from step S4; Step S6: Import the data obtained in step S2 into the phase field pre-processing program obtained in step S5; Step S7: Use the running result of the phase field pre-program in step S6 as the initial condition of the substrate / impurity phase, and import it into the constructed phase field model for describing the dendrite growth process to obtain a calculation model coupled with the substrate / impurity phase having a random domain. Step S8: Obtain the numerical solution of the computational model from step S7 through computer simulation.
2. The computational method for adding an irregular subdomain without flux boundary conditions within the computational domain according to claim 1, characterized in that: The specific process of establishing the computational equations for the computational method described in step S3, which couples an irregular subdomain with added flux-free boundary conditions within the computational domain, is as follows: Step S301: Derive the initial conditions of the Neumann boundary condition formula; Step S302: Derive the Neumann boundary conditions applied to the diffusion interface by the Neumann boundary conditions formula. Step S303: Derive the contact angle boundary condition formula of the Neumann boundary condition formula; Step S304: Introduce irregular subdomain parameters The contact angle boundary conditions are incorporated into the original governing equations; Step S305: List the calculation equations.
3. The computational method for adding an irregular subdomain without flux boundary conditions within the computational domain according to claim 2, characterized in that: In step S301, when deriving the initial conditions of the Neumann boundary condition formula, an arbitrary function is used. For example, its Laplace quantity Multiply by subdomain parameter We can obtain the following using the differential method of the product of identities: Thus, we obtained the same as Proportional calculation items; in, Indicates the gradient sign. For function Gradient in space, For subdomain parameters Gradient in space.
4. The computational method for adding an irregular subdomain without flux boundary conditions within the computational domain according to claim 3, characterized in that: In step S302, when the Neumann boundary condition formula is derived and applied to the diffusion interface, the inward unit normal vector is... The boundary is formed by The boundary conditions imposed by the Neumann boundary on the diffusion interface are given by the following formula: in, The direction of the Dharma.
5. The computational method for adding an irregular subdomain without flux boundary conditions within the computational domain according to claim 4, characterized in that: The specific process for deriving the contact angle boundary condition formula of the Neumann boundary condition formula in step S303 is as follows: Step S3031: For the AC and CH equations in the phase-field method, their free energy form is expressed as: in, Let be the free energy function. To define the phase field sequence parameters for different phases, For the free energy functional of a double potential well, The gradient energy coefficient, For the entire computational region; Step S3032, in the free energy function The variational derivative of the total free energy at the extreme point Therefore, we can conclude that: Step S3033: Multiply both sides of the formula obtained in step S3032 by ,get: The contact angle boundary condition equation is obtained as follows: 。 6. The computational method for adding an irregular subdomain without flux boundary conditions within the computational domain according to claim 5, characterized in that: The introduction of irregular subdomain parameters in step S304 When incorporating the contact angle boundary conditions into the original governing equations, the specific process is as follows: Step S3041: Over the entire computational domain, the contact angle θ is expressed as: Among them, subdomain parameters It satisfies the diffuse reflection boundary and satisfies , It is the unit normal vector of the phase interface; Step S3042: Derive the following equations for the contact angle boundary conditions: 。 7. The computational method for adding an irregular subdomain without flux boundary conditions within the computational domain according to claim 6, characterized in that: The specific process of listing the calculation equations in step S305 is as follows: Step S3051: Define the chemical potential driving morphological evolution as the variational derivative of the system's total free energy: Step S3052: Apply the calculation method to the chemical potential by multiplying it by the subdomain parameter. Applying the product rule and boundary condition equations, we obtain: For the conserved order parameter, its evolution is governed by the Cahn–Hilliard equation, where the rate of change of the order parameter is equal to the divergence of its flux and proportional to the gradient of the chemical potential: ; in, For time, This is the mobility coefficient; Step S3053: The calculation method yields the following: It is the flux that conserves the order parameter; where, For mass flux; This represents the flux perpendicular to the domain boundary; Step S3054, the formula for calculating the Cahn–Hilliard equation is written as follows: in, The contact angle boundary condition applied at the three-phase boundary. .
8. The computational method for adding an irregular subdomain without flux boundary conditions within the computational domain according to claim 7, characterized in that: For closed systems, .