Real-time evolution method of battery electrochemical characteristics based on digital twin model

Abstract

This invention relates to the field of battery testing technology, specifically a real-time evolution method for the electrochemical characteristics of batteries based on a digital twin model. The method includes simultaneously applying high-frequency polarization and ultrasonic pulses to the battery under test, acquiring acoustic-electric intermodulation distortion signals, surface thermal gradients, and ultrasonic attenuation delays; after signal separation, calculating the spatial gradient using a multidimensional physical boundary constraint tensor, and dynamically reconstructing a three-dimensional heterogeneous spatial domain; generating a virtual battery entity by mapping conductivity and mass transport parameters, driving transient acceleration derivation after loading a limit load spectrum, and outputting the virtual hotspot diffusion trajectory and capacity decay evolution curve; finally, extracting thermal runaway and aging inflection point features, generating a tiered degradation matching code, and executing it. This invention achieves precise battery sorting through acoustic-electric twin evolution.

Description

Technical Field

[0001] This invention relates to the field of battery testing technology, specifically a method for real-time evolution of battery electrochemical characteristics based on a digital twin model. Background Technology

[0002] With the rapid development of new energy vehicles and large-scale energy storage industries, the full life cycle management of batteries and the assessment of the residual utilization value of retired batteries have become key concerns in the industry. In order to ensure the safe operation of battery systems and accurately assess their remaining utilization value, it is necessary to effectively assess the internal health status and aging degree of batteries.

[0003] Currently, in practical engineering applications, conventional methods for detecting and evaluating the state of batteries typically focus on collecting macroscopic external physical characteristics such as terminal voltage, charging and discharging current, and casing surface temperature during battery operation using sensor systems. These macroscopic parameters are then substituted into pre-defined electrochemical empirical formulas or equivalent circuit models to estimate the battery's current overall capacity decay state or internal resistance change trend.

[0004] However, the operation inside a battery is essentially a complex process involving highly coupled multi-physics fields, encompassing ion mass transfer, charge transfer at the solid-liquid interface, and mechanical strain of porous electrodes. Under complex and variable dynamic operating conditions or extreme loads, the microscopic pore structure and electrochemical reaction interface inside the battery often exhibit significant spatial heterogeneity and highly nonlinear dynamic changes. Existing methods based on macroscopic external signals and traditional lumped parameter models, while reflecting the overall macroscopic performance of the battery, are limited by their technical principles and struggle to deeply capture and map the three-dimensional spatial distribution differences of microstructural defects and electrochemical kinetics within the battery, nor can they continuously characterize the dynamic accumulation and evolution of these microscopic features under varying operating conditions.

[0005] Therefore, how to achieve spatial characterization of the continuous dynamic evolution of the internal electrochemical characteristics and physical state of batteries under complex operating conditions has become a pressing technical problem to be solved in the field of in-depth battery condition assessment.

[0006] To address this, a method for real-time evolution of battery electrochemical characteristics based on a digital twin model is proposed. Summary of the Invention

[0007] The purpose of this invention is to provide a real-time evolution method for the electrochemical characteristics of batteries based on a digital twin model, achieving precise battery sorting through acoustic-electric twin evolution. This invention relates to the field of battery testing technology, specifically a real-time evolution method for the electrochemical characteristics of batteries based on a digital twin model. The method includes simultaneously applying high-frequency polarization and ultrasonic pulses to the battery under test, acquiring acoustic-electric intermodulation distortion signals, surface thermal gradients, and ultrasonic attenuation delays; after signal separation, calculating the spatial gradient using a multidimensional physical boundary constraint tensor, and dynamically reconstructing a three-dimensional heterogeneous spatial domain; generating a virtual battery entity by initializing with mapped conductivity and mass transport parameters, driving transient acceleration deduction after loading a limit load spectrum, and outputting the virtual hotspot diffusion trajectory and capacity decay evolution curve; finally, extracting thermal runaway and aging inflection point features, generating a tiered degradation matching code, and issuing it for execution.

[0008] To achieve the above objectives, the present invention provides the following technical solution: A real-time evolution method for battery electrochemical characteristics based on digital twin models includes: High-frequency polarization pulses and ultrasonic pulses of the same frequency and phase-locked loop are applied to the battery under test, and the acoustic-electric intermodulation nonlinear distortion signal, surface thermal gradient matrix and ultrasonic penetration attenuation time delay sequence are acquired simultaneously. Wavelet decomposition is applied to process the nonlinear distortion signal of acoustic-electric intermodulation, separating the high-frequency polarization component and the low-frequency acoustic-electric component; based on the preset three-dimensional finite element mesh, the surface thermal gradient matrix and the ultrasonic penetration attenuation delay sequence are projected onto the mesh nodes to construct a multidimensional physical boundary constraint tensor. The spatial gradient distribution of the multidimensional physical boundary constraint tensor is calculated, a microscopic solution mesh is generated in the structural anomaly region where the gradient exceeds the limit, and a three-dimensional heterogeneous spatial domain is reconstructed. The high-frequency polarization component, the low-frequency acoustic-electric component and the multidimensional physical boundary constraint tensor are substituted into the multiphysics field model established on the three-dimensional heterogeneous spatial domain to initialize and generate a virtual battery entity. Inject the extreme peak frequency modulation load spectrum into the virtual battery entity to drive the execution of accelerated time step integral derivation, and output the virtual hot spot diffusion trajectory and capacity decay evolution curve; extract the thermal runaway critical time of the virtual hot spot diffusion trajectory and the aging inflection point cycle number of the capacity decay evolution curve, calculate the tiered degradation matching code based on the two, and send it to the local sorting execution device.

[0009] Preferably, the process of acquiring the acoustic-electric intermodulation nonlinear distortion signal, surface thermal gradient matrix, and ultrasonic penetration attenuation delay sequence includes: extracting the mechanical intrinsic frequency characteristic value of the battery under test, and synchronously configuring the mechanical intrinsic frequency characteristic value as the reference resonant frequency of the electrical excitation circuit and the acoustic wave transmission circuit; controlling the output delay of the electrical excitation circuit and the acoustic wave transmission circuit through a hardware phase-locked loop to lock the zero-phase-difference state, and synchronously injecting high-frequency polarization pulses and ultrasonic pulses into the battery under test; acquiring the original voltage response sequence of the battery under test in the zero-phase-difference state, inputting the original voltage response sequence into the demodulation operator, and using the reference resonant frequency as... Demodulation is performed using a reference carrier; the fundamental frequency response component with the same reference resonant frequency as the original voltage response sequence is filtered out, and the harmonic components modulated by the acoustic-electric coupling effect are extracted. The harmonic components are defined as acoustic-electric intermodulation nonlinear distortion signals; the values ​​of the array temperature measurement nodes arranged on the surface of the battery under test are read, and the temperature change rate of adjacent nodes is calculated by applying a spatial interpolation function to reconstruct the surface thermal gradient matrix; the ultrasonic echo waveform penetrating the battery under test is captured, and the time difference vector is calculated by subtracting the high-frequency polarization pulse excitation time from the arrival time of the ultrasonic echo waveform envelope peak. The time difference vector is defined as the ultrasonic penetration attenuation time delay sequence.

[0010] Preferably, the process of separating the high-frequency polarization component and the low-frequency acoustic-electric component includes: introducing a preset orthogonal wavelet basis function, performing multi-scale discrete wavelet transform decomposition on the acoustic-electric intermodulation nonlinear distortion signal to generate a wavelet signal matrix containing multiple frequency band levels; extracting the highest frequency band detail coefficient sequence from the wavelet signal matrix and defining the highest frequency band detail coefficient sequence as the high-frequency polarization component; extracting the lowest frequency band approximation coefficient sequence from the wavelet signal matrix and defining the lowest frequency band approximation coefficient sequence as the low-frequency acoustic-electric component; calculating the wavelet energy spectral density feature value of the high-frequency polarization component within the time domain window, extracting the waveform envelope peak feature value of the low-frequency acoustic-electric component, and concatenating the wavelet energy spectral density feature value and the waveform envelope peak feature value to construct an evolution feature vector.

[0011] Preferably, the process of constructing the multidimensional physical boundary constraint tensor includes: extracting the spatial node coordinate system of a preset three-dimensional finite element mesh, and separating the surface contour node set and the internal volume node set from the spatial node coordinate system; aligning the surface thermal gradient matrix to the surface contour node set through coordinate interpolation transformation, and assigning values ​​to the surface contour node set to generate thermodynamic anomaly boundary weights; projecting the ultrasonic penetration attenuation delay sequence to the internal volume node set along the ultrasonic propagation and transmission path, calculating the local attenuation coefficient of the corresponding node and assigning values ​​to generate mechanical porosity degradation factors; extracting the topological adjacency matrix of the spatial node coordinate system, and concatenating the thermodynamic anomaly boundary weights and the mechanical porosity degradation factors along the dimensions of the topological adjacency matrix to construct the multidimensional physical boundary constraint tensor.

[0012] Preferably, the process of reconstructing the three-dimensional heterogeneous spatial domain includes: performing differential operations along the spatial coordinate axes of the multidimensional physical boundary constraint tensor, splicing partial derivative values ​​to generate a three-dimensional spatial gradient matrix; comparing the values ​​of the three-dimensional spatial gradient matrix with a preset gradient threshold, marking coordinate intervals with values ​​greater than the preset gradient threshold as structurally abnormal regions, and marking coordinate intervals with values ​​less than the preset gradient threshold as volumetric regions; triggering a mesh subdivision command in the structurally abnormal regions, cutting the preset three-dimensional finite element mesh into a high-node-density micro-solution mesh according to a preset geometric ratio; triggering a mesh coarsening command in the volumetric regions, merging the preset three-dimensional finite element mesh according to topological connectivity to generate a low-node-density volumetric macro-mesh; extracting the suspended boundary nodes between the micro-solution mesh and the volumetric macro-mesh, aligning the suspended boundary nodes through spatial coordinate interpolation, and stitching the boundary closure to generate a three-dimensional heterogeneous spatial domain.

[0013] Preferably, the initialization process for generating the virtual battery entity includes: extracting high-frequency polarization components and mapping and updating the charge transfer internal resistance reference value and interface double-layer capacitance value of each grid node in the three-dimensional heterogeneous spatial domain; extracting low-frequency acoustic-electric components and mapping and updating the solid-phase diffusion time constant and maximum lithium insertion / extraction concentration limit value of each grid node in the three-dimensional heterogeneous spatial domain; analyzing the multidimensional physical boundary constraint tensor and converting the thermodynamic anomaly boundary weights into the local thermal conductivity attenuation parameters of the corresponding grid nodes; and converting the mechanical porosity degradation factor in the multidimensional physical boundary constraint tensor into a parameter for... The effective reaction surface area reduction parameter of the grid nodes is used; the reference value of charge transfer internal resistance, the value of interface double layer capacitance, the solid-phase diffusion time constant, the maximum lithium insertion / extraction concentration limit, the local thermal conductivity decay parameter and the effective reaction surface area reduction parameter are extracted, and the initial grid node parameter matrix of the multiphysics model is loaded; the multiphysics model loaded with the initial grid node parameter matrix is ​​driven to perform the initial steady-state solution, and a multiphysics simulation instance object containing the initial state variables of the nodes in the three-dimensional heterogeneous spatial domain is constructed, and the multiphysics simulation instance object is defined as a virtual battery entity.

[0014] Preferably, the process of outputting the virtual hotspot diffusion trajectory and capacity decay evolution curve includes: analyzing the limiting peak frequency modulated load spectrum and extracting the time-domain transient current boundary excitation sequence; loading the time-domain transient current boundary excitation sequence onto the tab topology boundary node of the virtual battery entity, triggering the transient multiphysics coupling solver to perform adaptive time-step partial differential integral operation; within each time step of the adaptive time-step partial differential integral operation, traversing and extracting the three-dimensional spatial coordinates corresponding to the absolute extreme values ​​of the internal temperature scalar field of the virtual battery entity, and splicing the three-dimensional spatial coordinates in the time evolution order to construct the virtual hotspot diffusion trajectory; within each time step of the adaptive time-step partial differential integral operation, performing volume integral calculation on the global active lithium ion concentration scalar field of the virtual battery entity, extracting the global effective integrated capacity value, and fitting the global effective integrated capacity value sequence along the time axis to generate the capacity decay evolution curve.

[0015] Preferably, the process of sending the data to the local sorting execution device includes: performing a time-domain first-order differential operation on the virtual hotspot diffusion trajectory, extracting the simulation timestamp when the temperature rise slope first reaches a preset safety limit value, and defining it as the thermal runaway critical time; calculating the second derivative distribution of the capacity decay evolution curve, identifying the peak point with the largest absolute value in the second derivative sequence, and defining the evolution cycle number corresponding to the peak point as the aging inflection point cycle number; substituting the thermal runaway critical time and the aging inflection point cycle number as input variables into a preset tiered utilization classification logic table, and retrieving and matching the corresponding battery health level identifier; converting the battery health level identifier into hexadecimal machine code containing sorting fork action instructions, and defining it as the tiered degradation matching code; transmitting the tiered degradation matching code to the logic controller of the local sorting execution device through the communication bus, and driving the execution mechanism to perform sorting actions according to the physical channel corresponding to the tiered degradation matching code.

[0016] Compared with the prior art, the beneficial effects of the present invention are as follows: 1. By synchronously applying hardware phase-locked electrical and acoustic dual-physics field excitation to the battery under test, and using wavelet decomposition to process the acoustic-electric intermodulation distortion signal, this invention establishes a unified time zero point and initial phase reference interface by strictly implementing dual-field phase-locking at the injection boundary. This ensures that the phase lag generated by the sound wave during its propagation in the three-dimensional space inside the battery is no longer an interference, but is transformed into a stable physical characteristic characterizing internal microstructural defects and mechanical porosity degradation.

[0017] 2. By constructing a multidimensional physical boundary constraint tensor and performing adaptive mesh reconstruction based on spatial gradients, micro-mesh subdivision is applied to densely defective anomalous regions to capture evolutionary details, based on the distribution of thermal and acoustic physical tensors. Mesh coarsening is performed on stable volumetric regions, and dangling nodes are closed using interpolation functions. This three-dimensional heterogeneous spatial domain significantly reduces redundant solution nodes while fully preserving the dynamic evolution characteristics of local micro-defects, and reasonably balances the computational burden of multiphysics simulation with the ability to restore spatial states.

[0018] 3. By injecting extreme load spectrum into customized virtual entities for accelerated simulation, and converting the evolution results into local device execution instructions, the system outputs the critical time of thermal runaway reflecting the safety margin and the number of aging inflection point cycles reflecting lifespan mutations. Based on the above-mentioned forward-looking simulation indicators, the system automatically matches the hierarchical logic, compiles and generates hexadecimal tiered degradation matching codes, and directly sends them to the underlying controller to drive the physical fork action, thus shortening the business flow cycle of the identification of retired battery tiered utilization. Attached Figure Description

[0019] Figure 1 A schematic diagram of the real-time evolution method for battery electrochemical characteristics based on a digital twin model provided by this invention; Figure 2 This invention provides a schematic diagram of the process for reconstructing a three-dimensional heterogeneous spatial domain. Figure 3 This is a schematic diagram of the process for initializing and generating a virtual battery entity according to an embodiment of the present invention. Detailed Implementation

[0020] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0021] Please see Figures 1 to 3 This invention provides a method for real-time evolution of battery electrochemical characteristics based on a digital twin model, the technical solution of which is as follows: A method for real-time evolution of battery electrochemical characteristics based on digital twin models, the specific process of which is as follows: Figure 1 As shown, it includes: High-frequency polarization pulses and ultrasonic pulses of the same frequency and phase-locked loop are applied to the battery under test, and the acoustic-electric intermodulation nonlinear distortion signal, surface thermal gradient matrix and ultrasonic penetration attenuation time delay sequence are acquired simultaneously. Wavelet decomposition is applied to process the nonlinear distortion signal of acoustic-electric intermodulation, separating the high-frequency polarization component and the low-frequency acoustic-electric component; based on the preset three-dimensional finite element mesh, the surface thermal gradient matrix and the ultrasonic penetration attenuation delay sequence are projected onto the mesh nodes to construct a multidimensional physical boundary constraint tensor. The spatial gradient distribution of the multidimensional physical boundary constraint tensor is calculated, a microscopic solution mesh is generated in the structural anomaly region where the gradient exceeds the limit, and a three-dimensional heterogeneous spatial domain is reconstructed. The high-frequency polarization component, the low-frequency acoustic-electric component and the multidimensional physical boundary constraint tensor are substituted into the multiphysics field model established on the three-dimensional heterogeneous spatial domain to initialize and generate a virtual battery entity. Inject the extreme peak frequency modulation load spectrum into the virtual battery entity to drive the execution of accelerated time step integral derivation, and output the virtual hot spot diffusion trajectory and capacity decay evolution curve; extract the thermal runaway critical time of the virtual hot spot diffusion trajectory and the aging inflection point cycle number of the capacity decay evolution curve, calculate the tiered degradation matching code based on the two, and send it to the local sorting execution device.

[0022] Example 1: Furthermore, the process of acquiring the acoustic-electric intermodulation nonlinear distortion signal, surface thermal gradient matrix, and ultrasonic penetration attenuation delay sequence includes: extracting the mechanical intrinsic frequency characteristic value of the battery under test, and synchronously configuring the mechanical intrinsic frequency characteristic value as the reference resonant frequency of the electrical excitation circuit and the acoustic wave transmission circuit; controlling the output delay of the electrical excitation circuit and the acoustic wave transmission circuit through a hardware phase-locked loop to lock the zero-phase-difference state, and synchronously injecting high-frequency polarization pulses and ultrasonic pulses into the battery under test; acquiring the original voltage response sequence of the battery under test in the zero-phase-difference state, inputting the original voltage response sequence into the demodulation operator, and using the reference resonant frequency as... Demodulation is performed using a reference carrier; the fundamental frequency response component with the same reference resonant frequency as the original voltage response sequence is filtered out, and the harmonic components modulated by the acoustic-electric coupling effect are extracted. The harmonic components are defined as acoustic-electric intermodulation nonlinear distortion signals; the values ​​of the array temperature measurement nodes arranged on the surface of the battery under test are read, and the temperature change rate of adjacent nodes is calculated by applying a spatial interpolation function to reconstruct the surface thermal gradient matrix; the ultrasonic echo waveform penetrating the battery under test is captured, and the time difference vector is calculated by subtracting the high-frequency polarization pulse excitation time from the arrival time of the ultrasonic echo waveform envelope peak. The time difference vector is defined as the ultrasonic penetration attenuation time delay sequence.

[0023] First, to extract the mechanical intrinsic frequency characteristic values ​​of the battery under test, the system uses a wideband ultrasonic transducer to send a swept-frequency excitation signal to the battery (the swept-frequency excitation power remains constant, specifically set to a reference power of 5W to 15W). Multi-point triaxial vibration sensors located at the geometric center of the battery casing and the root of the tabs are used to collect the battery's resonant response. For cases with multiple resonant modes, the system extracts each resonant peak with an amplitude greater than a preset noise floor threshold (e.g., greater than 3dB) through frequency domain analysis, and identifies the fundamental frequency point with the lowest frequency and largest response amplitude as the mechanical intrinsic frequency characteristic value. Specifically, the frequency sweep range of the swept-frequency excitation signal is set to 20kHz to 500kHz, and the sweep step size is set to 100Hz to 500Hz. Subsequently, this characteristic value parameter is input into the Direct Digital Synthesizer (DDS) module of the Field Programmable Gate Array (FPGA). This module synchronously generates the reference resonant frequency carrier signal required for both the electrical excitation circuit and the acoustic wave transmission circuit, completing the basic configuration of the excitation source.

[0024] Next, synchronous injection of phase-locked loop pulses at the same frequency is performed. Using the digital phase-locked loop (DPLL) integrated within the FPGA or an external high-precision hardware phase-locked loop chip, phase discrimination is performed on the electrical polarization excitation signal and the ultrasonic emission signal, respectively. To obtain the true phase feedback signal, a wideband voltage follower is connected in parallel at the battery terminal tabs to acquire the actual phase of the high-frequency polarization voltage at the electrodes. Simultaneously, a miniature piezoelectric ceramic probe is embedded at the physical coupling interface between the ultrasonic transducer and the battery surface to acquire the actual acoustic wave vibration phase penetrating the battery surface. The two acquired actual feedback signals are shaped into square waves by a high-speed comparator and then input to the phase detector inside the FPGA for phase comparison. The system monitors the arrival of the two signals at the battery in real time. The true phase state of the electrode-surface coupling point is determined by a digital phase-locked loop dynamically adjusting the output delay clock of the internal digital-to-analog converter channel. This forces the electrical high-frequency polarization pulse and the ultrasonic pulse to maintain a zero phase difference at the physical injection boundary of the battery (i.e., the electrode terminal and the surface acoustic coupling point). Here, both the electrical excitation frequency and the acoustic wave frequency are uniformly set to the mechanical eigenfrequency. The physical mechanism is as follows: utilizing the acoustic-electric resonance enhancement effect, the local strain field generated by mechanical vibration modulates the ion migration process of the electrochemical double layer to the greatest extent. In this state of synchronous frequency, the periodic fluctuation of the dielectric density caused by the acoustic wave and the change in polarization intensity caused by the electrical pulse produce nonlinear interference, thereby enabling the extraction of intermodulation components reflecting internal microscopic defects with an extremely high signal-to-noise ratio. Specifically, the synchronous frequency refers to the fundamental frequency of the electrical high-frequency polarization pulse being consistent with the transmission center frequency of the ultrasonic transducer. Although the ultrasonic wave will have an inherent phase lag due to the difference in sound velocity in the dielectric when penetrating the microscopic interface inside the battery, boundary phase-locking ensures the consistency of this lag law. In this phase-locked state, the system continuously injects dual-physical field excitation into the battery through a high-frequency constant current source and an ultrasonic emission drive circuit.

[0025] In the signal acquisition and acoustic-electric intermodulation distortion component separation stage, the system utilizes a high-speed analog-to-digital converter (ADC) at a sampling rate far exceeding the reference resonant frequency to synchronously acquire the original voltage response sequence across the battery terminals under zero phase difference conditions. This sequence is transmitted to the digital signal processor (DSP). The DSP calls a digital quadrature demodulation operator to perform time-domain multiplication and low-pass filtering operations on the original voltage response sequence and a pre-stored sine / cosine reference carrier at the reference resonant frequency, completing the digital frequency shift of the signal. Subsequently, the DSP calls a digital band-stop filter to filter out the fundamental frequency response component with a frequency equal to the reference resonant frequency and extracts the harmonic components generated by nonlinear modulation due to internal acoustic-electric coupling effects. The extracted pure harmonic components are then processed. The sequence is stored and defined as an acoustic-electric intermodulation nonlinear distortion signal. The harmonic components are specifically set as the second harmonic (2f0) and third harmonic (3f0) components of the reference resonant frequency. In the filter settings, the center frequency of the digital bandstop filter is set to the reference resonant frequency f0, the stopband bandwidth is set to ±0.05f0, and a digital bandpass filter with a cutoff frequency of 1.5f0 to 3.5f0 is cascaded in the subsequent stage. The physical mechanism for extracting the 2f0 and 3f0 components is that when ultrasonic waves propagate through the electrode / electrolyte interface, the alternating mechanical stress generated by the sound pressure induces a periodic change in the thickness of the double layer. This mechanical modulation effect is nonlinearly superimposed with the electrochemical response generated by the high-frequency polarization current. As the battery ages (such as SEI film thickening or local lithium plating), the nonlinear stiffness of the interface increases, which leads to a significant enhancement of the second and third harmonic energy generated by acoustic-electric coupling. Therefore, by monitoring the amplitude changes of the 2f0 and 3f0 components, the degree of physical degradation of the internal interface of the battery can be quantitatively identified.

[0026] To reconstruct the surface thermal gradient matrix, a multi-channel flexible thin-film thermocouple array or a miniature thermistor matrix is ​​attached to the surface of the battery casing under test. The data acquisition card scans and reads the absolute temperature values ​​of each physical temperature measurement node according to a preset cycle. After receiving the discrete temperature node data, the host computer software calls a two-dimensional bicubic interpolation or Kriging spatial interpolation algorithm to generate virtual encrypted nodes between the physical temperature measurement nodes and calculates the temperature change rate of any two adjacent grid nodes in the spatial coordinate direction. Finally, the surface thermal gradient matrix covering the entire battery surface is generated according to the spatial topological relationship. The temperature change rate mentioned here refers to the spatial temperature gradient. The specific process is as follows: along the X-axis and Y-axis directions of the preset three-dimensional orthogonal coordinate system, the temperature scalar values ​​of the virtual encrypted nodes are subjected to first-order spatial difference operations to obtain the partial derivative matrix in the X-direction and the partial derivative matrix in the Y-direction, and synthesizes the surface thermal gradient matrix representing the distribution of heat flow potential energy.

[0027] Finally, during the acquisition of the ultrasonic penetration attenuation delay sequence, the receiving transducer positioned on the opposite side of the ultrasonic transmitter of the battery captures the ultrasonic echo waveform after penetrating the internal medium of the battery. The system calls envelope extraction algorithms such as Hilbert transform to obtain the low-frequency envelope of the high-frequency echo waveform, and uses a peak detection algorithm to accurately locate the absolute arrival timestamp corresponding to the point of maximum envelope amplitude. The high-precision digital timer inside the system's main control board synchronously records the initial excitation timestamp of the rising edge triggered by the high-frequency polarization pulse. By subtracting the excitation timestamp from the arrival timestamp of the envelope peak, the time difference vector of a single excitation cycle is obtained. The time difference vectors recorded in multiple consecutive test cycles are spliced ​​together in chronological order along the time axis to form a complete ultrasonic penetration attenuation delay sequence.

[0028] By introducing a hardware phase-locked loop to force the acoustic and electrical dual-physics fields to operate synchronously in a zero-phase-difference state, the physical consistency of the internal excitation conditions is ensured. Demodulation filtering techniques are used to extract the modulated harmonic components, effectively separating the nonlinear characteristics of acoustic-electric coupling reflecting the complex internal interface of the battery. Simultaneously, by combining the spatially interpolated reconstructed thermal gradient matrix with the ultrasonic penetration attenuation delay sequence, discrete single-point monitoring data are expanded into continuous spatially distributed parameters.

[0029] Furthermore, the process of separating the high-frequency polarization component and the low-frequency acoustic-electric component includes: introducing a preset orthogonal wavelet basis function, performing multi-scale discrete wavelet transform decomposition on the acoustic-electric intermodulation nonlinear distortion signal to generate a wavelet signal matrix containing multiple frequency band levels; extracting the highest frequency band detail coefficient sequence from the wavelet signal matrix and defining the highest frequency band detail coefficient sequence as the high-frequency polarization component; extracting the lowest frequency band approximation coefficient sequence from the wavelet signal matrix and defining the lowest frequency band approximation coefficient sequence as the low-frequency acoustic-electric component; calculating the wavelet energy spectral density feature value of the high-frequency polarization component within the time domain window, extracting the waveform envelope peak feature value of the low-frequency acoustic-electric component, and concatenating the wavelet energy spectral density feature value and the waveform envelope peak feature value to construct an evolution feature vector.

[0030] First, Daubechies wavelet basis functions with tight support and orthogonality are preloaded, for example, the db6 wavelet basis is selected. The discrete sequence of the acoustic-electric intermodulation nonlinear distortion signal is read from the buffer memory, and the standard Mallat tower fast algorithm is called. Through cascaded low-pass and high-pass digital filter banks, the discrete sequence is subjected to stepwise downsampling and orthogonal filtering operations of a preset depth (6 layers in this embodiment), thereby generating and storing a wavelet signal matrix containing different frequency band resolution characteristics from high to low in the system memory. Specifically, the sampling frequency of the discrete sequence of the acoustic-electric intermodulation nonlinear distortion signal is set to a preset value fs. According to the Nyquist sampling theorem and wavelet frequency band division rules, the frequency band range after the first layer of high-pass filtering decomposition corresponds to fs / 4 to fs / 2, covering the frequency band of high-frequency polarization pulses; at the same time, the frequency band range corresponding to the output of the 6th layer of low-pass filtering is 0 to fs / 128, covering the low-frequency band modulated by acoustic-electric coupling.

[0031] Subsequently, the generated wavelet signal matrix is ​​subjected to target-level separation and extraction. The output node after the first high-pass filter decomposition is directly located, and the detailed coefficient sequence containing the highest frequency fluctuation information of the original signal under this node is extracted. This sequence is then peeled off and defined as a high-frequency polarization component. At the same time, the output node of the deepest decomposition layer (i.e., the 6th layer) is located, and the approximate coefficient sequence representing the global slow change trend of the signal under this node is extracted. This sequence is then independently defined as a low-frequency acoustic-electric component.

[0032] In the high-frequency feature extraction step, a sliding time-domain window containing a fixed number of data points is set along the time axis of the discrete sequence for the separated high-frequency polarization components. Within each translation period of the time-domain window, the amplitude of all discrete data points covered by the window is squared, and all squared values ​​within the window are accumulated. This sum is then divided by the total number of data points within the window. By calculating the mean square value within this time period, the wavelet energy spectral density feature value of the corresponding high-frequency polarization component within the time-domain window is obtained. The specific formula is as follows: in, This represents the wavelet energy spectral density eigenvalue of the high-frequency polarization component obtained within the k-th sliding time-domain window; This indicates the length of the fixed number of data points contained in the set sliding time-domain window; This represents the amplitude of the nth discrete data point in the discrete sequence of the separated high-frequency polarization components; Indicates the starting data point position of the current sliding time domain window; This means squaring the magnitudes of all discrete data points covered by the window and summing all the squared values ​​within the window.

[0033] Subsequently, the generated wavelet signal matrix is ​​subjected to target-level separation and extraction. The output node after the first high-pass filter decomposition is located, and the detailed coefficient sequence containing the highest frequency fluctuation information of the original signal under this node is extracted. This sequence is then peeled off and defined as a high-frequency polarization component. At the same time, the output node of the deepest decomposition layer (i.e., the 6th layer) is located, and the approximate coefficient sequence representing the global slow change trend of the signal under this node is extracted. This sequence is then independently defined as a low-frequency acoustic-electric component.

[0034] In the low-frequency feature extraction step, for the separated low-frequency acoustic-electric components, a discrete Hilbert transform is performed on the approximate coefficient sequence to construct the corresponding analytical signal sequence. The instantaneous magnitude of the analytical signal is calculated point by point, and a smooth envelope sequence reflecting the low-frequency fluctuation state of the signal is fitted to generate a smooth envelope sequence. In order to eliminate the boundary truncation effect and sudden noise interference generated at both ends of the finite-length data sequence by the discrete Hilbert transform, the smooth envelope sequence is subjected to edge removal and filtering: data points at the beginning and end of the envelope sequence with a preset proportion (5% in this embodiment) are removed, and the moving average filter operator is applied to smooth the truncated envelope sequence twice. Subsequently, a global extremum search algorithm is used to traverse the entire envelope sequence, identify and lock the absolute convex point with the largest amplitude in the sequence, and extract the value corresponding to the point as the waveform envelope peak feature value of the low-frequency acoustic-electric components.

[0035] Finally, the combination and splicing of the feature parameters are performed. The calculated high-frequency polarization component wavelet energy spectral density feature value is used as the first element of the one-dimensional vector, and the extracted low-frequency acoustic-electric component waveform envelope peak feature value is used as the second element. The two are directly spliced ​​along the data dimension to construct a one-dimensional feature data set, which is used as the output evolution feature vector.

[0036] By introducing wavelet multi-scale decomposition at a specific level, the method effectively decouples complex acoustic-electric nonlinear distortion signals into physically meaningful high-frequency polarization and low-frequency acoustic-electric components. Combining this with sliding time-domain window quantization of high-frequency energy spectral density and employing endpoint removal and sliding smoothing mechanisms when extracting the low-frequency envelope, the method significantly reduces data redundancy while preserving core dynamic features, providing physically clear parameter inputs for subsequent multiphysics mapping in digital twin models.

[0037] Furthermore, the process of constructing the multidimensional physical boundary constraint tensor includes: extracting the spatial node coordinate system of a preset three-dimensional finite element mesh, and separating the surface contour node set and the internal volume node set from the spatial node coordinate system; aligning the surface thermal gradient matrix to the surface contour node set through coordinate interpolation transformation, and assigning values ​​to the surface contour node set to generate thermodynamic anomaly boundary weights; projecting the ultrasonic penetration attenuation delay sequence to the internal volume node set along the ultrasonic propagation and transmission path, calculating the local attenuation coefficient of the corresponding node and assigning values ​​to generate mechanical porosity degradation factors; extracting the topological adjacency matrix of the spatial node coordinate system, and concatenating the thermodynamic anomaly boundary weights and the mechanical porosity degradation factors along the dimensions of the topological adjacency matrix to construct the multidimensional physical boundary constraint tensor.

[0038] A three-dimensional geometric solid model of the battery to be tested is pre-established. The model is geometrically discretized using a hexahedral or tetrahedral meshing strategy. A three-dimensional finite element mesh is constructed, and a global three-dimensional orthogonal coordinate system is established. The spatial coordinate data of all geometric nodes in the mesh are extracted. A boundary recognition algorithm is executed through the facet sharing rules of the element topology: nodes contained in facets belonging only to the surface of a single solid element are extracted and divided into a surface contour node set; the remaining nodes that are shared by multiple solid elements are extracted and divided into an internal volume phase node set.

[0039] The surface thermal gradient matrix reconstructed in the previous steps and its corresponding two-dimensional measurement point coordinate system are obtained. Through a three-dimensional spatial coordinate affine transformation, this two-dimensional measurement point coordinate system is rigidly registered and aligned to the surface contour node set of the three-dimensional finite element mesh. Specifically, a surface unfolding algorithm based on the geometric topology of the battery casing under test is introduced to unfold the surface contour node set of the three-dimensional finite element mesh and map it to the global two-dimensional UV coordinate system. Then, using the physical origin of the two-dimensional measurement point coordinate system and the center point of the reference plane of the UV coordinate system as alignment anchor points, in-plane rotation and translation affine transformations are performed to make the discrete measurement points coincide with the unfolded mesh surface, completing the transformation from two-dimensional measurement points to three-dimensional... The reverse spatial localization of the surface involves setting a spatial search radius centered on the coordinates of each target node in the surface contour node set. Inverse distance weighted interpolation or spatial bilinear interpolation algorithms are used to smoothly map the measurement point values ​​in the thermal gradient matrix to the target mesh node. Subsequently, the absolute value of the mapped node's thermal gradient is compared with a preset baseline healthy thermal gradient threshold to generate a dimensionless thermodynamic anomaly boundary weight, which is then assigned as an independent attribute parameter to the corresponding surface contour node. The baseline healthy thermal gradient threshold is not a fixed constant but is adjusted based on the current injected load current of the battery under test and the environmental conditions. Temperature is a benchmark value dynamically retrieved by consulting a pre-stored MAP (Metrological Map of Healthy New Batteries) in the system. Specifically, the construction of the MAP is an offline calibration process, which includes the following steps: acquiring several brand-new, healthy battery samples of the same batch and model as the battery under test; placing the battery samples in a high-low temperature alternating controllable environment test chamber and setting multiple discrete environmental temperature test nodes covering the operating range of the battery under test; and at each set environmental temperature node, sequentially injecting multiple preset gradient load currents into the battery samples using a battery charge-discharge test system. After the battery sample reaches thermal equilibrium, the absolute temperature of its surface is simultaneously collected using a multi-channel temperature measurement array, and the steady-state extreme value of the surface thermal gradient is calculated. The discrete data point matrix obtained from the test is input into the host computer with the ambient temperature test node as the X-axis scalar, the load current ratio as the Y-axis scalar, and the corresponding steady-state extreme value of the surface thermal gradient as the Z-axis scalar. Two-dimensional bicubic interpolation is used to perform continuous smoothing processing on the discrete data point matrix to generate a three-dimensional feature surface data table within the full operating range. This data table is compiled and solidified into the non-volatile memory of the system and defined as the thermodynamic standard feature map of the healthy new battery.

[0040] The corresponding structural variation degree is derived from the pre-calibrated sound velocity-porosity attenuation mapping curve. The topological adjacency matrix representing the spatial connectivity between all nodes in the three-dimensional finite element mesh is extracted. The sound velocity-porosity attenuation mapping curve is established through previous offline calibration experiments: ultrasonic time-of-flight tests and X-ray computed tomography porosity analysis are performed simultaneously on multiple sets of battery samples from the same batch at different aging stages. The nonlinear discrete corresponding scatter points of sound velocity and microporosity at a specific ultrasonic frequency are extracted, and a continuous mapping curve is generated using a polynomial least squares fitting algorithm. The node sequence number of the topological adjacency matrix is ​​used as the reference dimension. The three-dimensional coordinate scalar of each node, the newly generated thermodynamic anomaly boundary weight scalar, and the mechanical porosity degradation factor scalar are used as different physical feature channels. The above physical feature channels are matrix-stitched along the longitudinal structure according to the node sequence number, and the single topological matrix is ​​upgraded to a multi-feature collection containing node spatial position, surface thermodynamic constraints, and internal mechanical degradation state, thus constructing a multi-dimensional physical boundary constraint tensor.

[0041] By deeply fusing three-dimensional geometric topology with multi-source sensing data, the spatial mapping of discrete surface thermal field and internal acoustic attenuation characteristics to finite element nodes is realized. Using offline calibrated feature pulse spectrum and mapping curve, the physical consistency of thermodynamic weights and mechanical degradation factors is ensured. This process transforms heterogeneous physical signals into a unified spatial constraint tensor, providing high-fidelity boundary conditions and data foundation for the subsequent multi-physics coupling evolution of digital twin models in the three-dimensional heterogeneous spatial domain.

[0042] Further, the process of reconstructing the three-dimensional heterogeneous spatial domain includes: performing differential operations along the spatial coordinate axes of the multidimensional physical boundary constraint tensor, and splicing partial derivative values ​​to generate a three-dimensional spatial gradient matrix; comparing the values ​​of the three-dimensional spatial gradient matrix with a preset gradient threshold, marking coordinate intervals with values ​​greater than the preset gradient threshold as structurally abnormal regions, and marking coordinate intervals with values ​​less than the preset gradient threshold as volumetric regions; triggering a mesh subdivision command in the structurally abnormal regions, cutting the preset three-dimensional finite element mesh into a high-node-density micro-solution mesh according to a preset geometric ratio; triggering a mesh coarsening command in the volumetric regions, merging the preset three-dimensional finite element mesh according to topological connectivity to generate a low-node-density volumetric macro-mesh; extracting the dangling boundary nodes between the micro-solution mesh and the volumetric macro-mesh, aligning the dangling boundary nodes through spatial coordinate interpolation, and stitching the boundary closure to generate a three-dimensional heterogeneous spatial domain. The specific process is as follows: Figure 2 As shown.

[0043] For the constructed multidimensional physical boundary constraint tensor, discrete central difference operations are performed along the three axes of the defined three-dimensional orthogonal spatial coordinate system. For any internal node in the tensor field, the physical constraint coefficient values ​​of its adjacent nodes in both positive and negative directions of a single coordinate axis are obtained. By calculating the difference between the eigenvalues ​​of these two adjacent nodes and dividing it by the spatial physical distance between them, the first-order partial derivative value of the point in the direction of the coordinate axis is obtained. After completing the difference calculation in all directions by traversing the three spatial coordinate axes in turn, the partial derivatives in these three directions are orthogonally vector synthesized according to the coordinate system vector operation rules to obtain the absolute value of the synthesized spatial gradient at the node. The above operation is performed on all nodes in the tensor space, and the absolute values ​​of the spatial gradients calculated at each point are rearranged according to the original spatial topology to construct a complete three-dimensional spatial gradient matrix.

[0044] After obtaining the generated three-dimensional spatial gradient matrix, a pre-set global gradient threshold parameter is introduced to calibrate the entire matrix point by point. The pre-setting method of the global gradient threshold parameter is as follows: calculate the statistical average and standard deviation of the absolute values ​​of the spatial gradients of all nodes in the three-dimensional spatial gradient matrix, and add the average value to the standard deviation by a preset multiple (such as 1.5 to 3 times) as the global gradient threshold parameter under the specific working condition to adaptively match the overall characteristic fluctuation level of the current tensor field. During the traversal comparison process, the relationship between the absolute value of the spatial gradient of each node and the preset gradient threshold is judged. If the absolute value of the gradient of a node is greater than the given preset gradient threshold, the node and the finite element geometric unit to which it belongs are marked as a structurally abnormal region on the topology map. If the absolute value of the gradient of a node does not exceed the preset gradient threshold, the corresponding unit is marked as a volumetric region. Through the region marking operation, the three-dimensional finite element mesh space is divided into two discrete coordinate interval sets: one with drastic feature changes and the other with stable feature evolution.

[0045] After marking the regional attributes, a local mesh subdivision command is triggered for the set of mesh units identified as structurally anomalous regions. A geometrically adaptive cutting algorithm is used to reduce the dimensionality of the original geometric units within these anomalous regions according to a pre-set subdivision ratio. A single hexahedral parent unit is divided into eight smaller sub-units along its three central axes, or a tetrahedral unit is cut along the midpoints of its edges. During the cutting process, new spatial nodes are inserted on the existing geometric edges and faces, and the physical coordinates of the new nodes are generated by linear interpolation based on the coordinate positions of the surrounding original nodes. During the local subdivision iteration, the gradient range (i.e., the maximum difference in the absolute values ​​of the gradients of all nodes within the sub-unit) or the geometric feature dimensions of the sub-unit are monitored in real time. When the gradient range within the sub-unit is less than a preset smoothing threshold, or when the volume of the sub-unit reaches a preset lower limit for the micro-mesh volume, the subdivision iteration for that unit is terminated. Through local subdivision iteration, a micro-solution mesh set with a spatial scale smaller than the initial mesh and a dense node distribution is constructed within the structurally anomalous region.

[0046] For the set of mesh elements marked as volumetric regions, a region mesh coarsening command is triggered. An element merging operation is performed along the connectivity path of the geometric topology. Multiple basic elements adjacent to each other within these regions and sharing the same common geometric surface are retrieved. By removing their internal common surfaces and related center nodes, discrete small-volume elements are wrapped and merged into a large-volume super-element. During the merging operation, the retained vertex coordinates are re-smoothed and repositioned according to the volume weights of adjacent meshes. Specifically, the original centroid coordinates and corresponding element volumes of each basic element participating in the merging are extracted. The centroid coordinates of each basic element are multiplied by its own volume to obtain a volume moment vector. Then, the volume moment vectors of all basic elements are vector-accumulated and divided by the total volume of the basic elements participating in the merging. The calculated three-dimensional coordinate values ​​are used as the vertex coordinates of the re-smoothed super-element. Through a geometric simplification process based on topological clustering, a low-node-density macroscopic mesh set with large individual volume and sparse spatial node distribution is generated in the volumetric region.

[0047] Because the structurally anomalous region and the volumetric region underwent mesh reconstruction operations at different scales, the entire mesh topology connectivity table was traversed to identify and extract suspended boundary nodes that were connected to the micro-mesh on only one side and lost the direct topology of the macro-mesh on the other side. Using the geometric patch of the volumetric macro-mesh located at the interface as the interpolation reference plane, the vertex coordinates and shape functions of the reference plane were extracted. Using the area coordinate interpolation algorithm, the spatial position of the suspended node was projected and registered onto the adjacent macro-scale geometric surface, establishing a multi-point constraint topology binding relationship across scales. Specifically, the binding relationship is as follows: the set of macro-master nodes adjacent to the suspended node on the reference plane is extracted, and the physical state variables (such as temperature and electric potential) of the suspended node are forcibly defined as a linear combination of the physical state variables of these macro-master nodes. The weight coefficient of the linear combination is strictly equal to the area interpolation shape function value of each macro-master node when the suspended node is projected onto the reference plane. In this way, the degrees of freedom of the suspended node are eliminated from the global solution matrix. After completing the geometric alignment and topological constraint definition of all suspended boundary nodes, the transition boundaries with gaps are stitched together to construct a three-dimensional heterogeneous spatial domain with closed boundaries.

[0048] By accurately identifying anomalous regions where the internal structural characteristics of the battery undergo drastic changes, the mesh is refined in these anomalous regions to ensure the accuracy of local multiphysics solutions, while the mesh is coarsened in stable regions to significantly reduce computational redundancy. Furthermore, multi-point constraint binding is applied to suspended boundary nodes, directly eliminating additional degrees of freedom in the solution process. This achieves a reasonable balance between simulation accuracy and computational efficiency in the digital twin model while ensuring the geometric closure and physical continuity of the cross-scale mesh.

[0049] Furthermore, the initialization process for generating the virtual battery entity includes: extracting high-frequency polarization components and mapping and updating the charge transfer internal resistance reference value and interface double-layer capacitance value of each grid node in the three-dimensional heterogeneous spatial domain; extracting low-frequency acoustic-electric components and mapping and updating the solid-phase diffusion time constant and maximum lithium insertion / extraction concentration limit value of each grid node in the three-dimensional heterogeneous spatial domain; analyzing the multidimensional physical boundary constraint tensor and converting the thermodynamic anomaly boundary weights into the local thermal conductivity attenuation parameters of the corresponding grid nodes; and converting the mechanical porosity degradation factor in the multidimensional physical boundary constraint tensor into the corresponding grid node's... The effective reactive surface area reduction parameter of the lattice nodes is calculated; the baseline value of charge transfer internal resistance, the value of interface double layer capacitance, the solid-phase diffusion time constant, the maximum lithium insertion / extraction concentration limit, the local thermal conductivity decay parameter, and the effective reactive surface area reduction parameter are extracted; the initial grid node parameter matrix of the multiphysics model is loaded; the multiphysics model loaded with the initial grid node parameter matrix is ​​driven to perform initial steady-state solution; a multiphysics simulation instance object containing the initial state variables of all nodes in the three-dimensional heterogeneous spatial domain is constructed; the multiphysics simulation instance object is defined as a virtual battery entity, and the specific process is as follows: Figure 3 As shown.

[0050] First, the extracted high-frequency polarization components are used as regulating factors for local electrochemical activity. Based on a pre-defined electrochemical parameter benchmark library of the battery system under test, the initial charge transfer resistance and interfacial double-layer capacitance benchmark values ​​corresponding to each grid node are retrieved. Utilizing the intensity differences of the high-frequency polarization components in the time and frequency domains, weight compensation calculations are performed on the benchmark values ​​of each node in the three-dimensional heterogeneous spatial domain. For grid regions marked as structurally abnormal, the charge transfer resistance of nodes in that region is amplified based on the proportion of high-frequency characteristic values ​​deviating from the normal range, and the double-layer capacitance is dynamically reduced. The process involves mapping the microscopic polarization features to the conductivity parameters of the grid nodes. The specific mapping process is as follows: the ratio of the extracted high-frequency polarization component feature value to the pre-stored reference high-frequency feature value of a healthy battery of the same type is calculated, and this ratio is defined as the polarization offset coefficient. For grid regions marked as structurally abnormal, the initial charge transfer internal resistance reference value of the node is directly multiplied by the polarization offset coefficient to implement gain amplification. At the same time, the initial interface double layer capacitance reference value of the node is divided by the polarization offset coefficient to implement dynamic scaling, thereby completing the quantitative assignment of conductivity parameters.

[0051] Next, the mass transport parameters of the nodes are updated using the eigenvalues ​​of the low-frequency acoustic-electric components. The waveform envelope peak value in the low-frequency acoustic-electric components is extracted and correlated with the solid-phase lithium-ion diffusion kinetics. Based on the material interface state reflected by the acoustic-electric coupling effect, the solid-phase diffusion time constant of each grid node in the three-dimensional heterogeneous spatial domain is mapped and adjusted. Simultaneously, based on the propagation attenuation characteristics of the acoustic signal in the porous electrode, the effective lithium intercalation space of the active material is inverted, and the maximum lithium intercalation / deintercalation concentration limit corresponding to each grid node is updated. This process ensures that the low-frequency acoustic information can accurately characterize the ion transport obstacles and energy within the battery. The specific process of correlation adjustment and inversion is as follows: A mapping correlation table between low-frequency envelope peak values ​​and solid-phase diffusion time constants, pre-established through offline electrochemical impedance spectroscopy calibration, is called. The currently extracted waveform envelope peak characteristics are input, and the corresponding solid-phase diffusion time constant is directly obtained through linear interpolation. Simultaneously, the increase ratio of the ultrasonic penetration attenuation delay of the current grid node relative to the baseline healthy attenuation delay is calculated. 1 is subtracted from this increase ratio to obtain the effective volume reduction factor. The initial maximum lithium intercalation / deintercalation concentration limit is multiplied by this effective volume reduction factor to obtain the updated maximum lithium intercalation / deintercalation concentration limit.

[0052] Specifically, the offline construction process of the mapping relationship table between the low-frequency envelope peak and the solid-phase diffusion time constant, which was previously established through offline electrochemical impedance spectroscopy calibration, is as follows: Several sample batteries of the same model and batch as the battery under test, but in different health states (covering the range from brand new to scrap extreme values), are selected; under a preset standard ambient temperature, electrochemical impedance spectroscopy tests are sequentially performed on each sample battery using an electrochemical workstation (test frequency range set to 10mHz to 10kHz, AC disturbance voltage amplitude set to 5mV, no DC bias); based on the equivalent circuit model, the low-frequency Weber impedance data in the Nyquist plot obtained from the test is fitted, and the current true solid-phase diffusion time constant of each sample battery is extracted; under the same test environment... Simultaneously apply high-frequency polarized pulses and ultrasonic pulses of the same frequency and phase-locked loop, with parameters identical to those of the online detection (including the amplitude and pulse width of the polarization pulse, the sound pressure intensity of the ultrasonic pulse, etc.), to each of the above sample batteries. Acquire the nonlinear distortion signal of acoustic-electric intermodulation and separate the low-frequency acoustic-electric components, and extract the corresponding waveform envelope peak feature values. Use the extracted waveform envelope peak feature values ​​as independent variables and the corresponding solid-phase diffusion time constant obtained analytically as dependent variables to construct a two-dimensional data scatter array. Perform polynomial curve fitting on the two-dimensional data scatter array using the least squares method to generate a continuous mathematical mapping function, or generate a discrete data interpolation lookup table according to a preset step size. Solidify it into the non-volatile memory of the local controller and define it as a mapping association table.

[0053] Subsequently, the multidimensional physical boundary constraint tensor is transformed into physical parameters. The thermodynamic anomaly boundary weights stored in the tensor are analyzed and converted into local thermal conductivity attenuation parameters of the corresponding grid nodes according to an inverse proportional mapping relationship. This allows the spatial distribution of the surface thermal gradient to be quantified as anisotropic degradation of the internal thermal conductivity of the battery. Simultaneously, the mechanical porosity degradation factor in the tensor is analyzed. Based on the theory of effective medium in porous electrodes, the evolution of porosity is transformed into the effective reactive surface area reduction parameter of the corresponding grid nodes. This characterizes the loss of the chemical reaction interface due to mechanical strain or microstructure damage. The specific calculation process for the physical parameter transformation is as follows: First, based on the introduced thermodynamic anomaly boundary weights, the local thermal conductivity attenuation parameter of the corresponding grid nodes is calculated. The specific calculation formula is as follows: in, This represents the calculated local thermal conductivity attenuation parameter. This represents the initial local thermal conductivity reference value for this grid node. This represents the boundary weight of the thermodynamic anomaly obtained through analysis; This represents the preset attenuation penalty factor.

[0054] Wherein, the preset attenuation penalty factor This is an empirical coefficient characterizing the sensitivity of a material's thermal conductivity to temperature gradients. Specifically, it is obtained by first measuring the true local thermal conductivity of sample batteries of the same type under different known thermodynamic anomaly boundary weights using a thermal conductivity meter. The measured values ​​are then substituted into the above formula for inverse least-squares fitting, thereby calibrating the degradation penalty factor specific to that battery model. For different battery systems and aging states, this degradation penalty factor... The rules for determining the value are as follows: For lithium iron phosphate battery systems: due to their relatively high thermal stability, their thermal conductivity is moderately sensitive to microstructural distortions. The preferred value range is 0.25 to 0.45; for ternary lithium battery systems: due to their high energy density and high risk of thermal runaway, the thermal conductivity of the material decreases more drastically with microstructural failure (such as electrolyte oxidation and decomposition). The preferred value range is 0.50 to 0.75. Under online testing conditions, the system performs adaptive adjustment based on the estimated state of health (SOH) of the battery under test: when the battery is in the early stage of its life (SOH > 80%), the lower limit of the corresponding system range is taken to reflect the relative integrity of the material structure; when the battery enters the late stage of aging (SOH < 80%), the value is linearly increased as SOH decreases. Up to the upper limit of the corresponding system range.

[0055] Simultaneously, the current abnormal porosity characterized by the mechanical porosity degradation factor is extracted, and the effective reaction surface area reduction parameter of the corresponding grid node is calculated based on the effective medium theory of porous electrodes. The specific calculation formula is as follows: in, This represents the calculated effective reaction surface area reduction parameter; This represents the initial effective reaction surface area of ​​the grid node; This represents the current anomalous porosity characterized by the mechanical porosity degradation factor; This indicates the pre-calibrated initial healthy porosity; This represents a preset empirical index. In this embodiment, for the NCM ternary lithium battery system, this index... The preset value is 1.5. This value is obtained by least-squares fitting calibration based on the microstructure reconstruction results of multiple sets of batteries with different aging levels using FIB-SEM (Focused Ion Beam Scanning Electron Microscopy), combined with measured electrochemical impedance spectroscopy data. In practical applications, the tortuosity of the battery cathode material will vary. The value of n is usually between 1.2 and 1.8. By selecting n=1.5, the computational load and the characterization accuracy of specific surface area loss can be well balanced.

[0056] Finally, the multiphysics model is initialized and entities are generated. In this embodiment, the preset multiphysics model is a mechanistic model composed of multiple sets of governing equations that follow the laws of battery electrochemistry and thermodynamics. This model integrates the charge conservation mechanism describing the potential spatial distribution, the diffusion mechanism describing lithium-ion mass transfer, and the heat conduction mechanism describing energy fluctuations. The initial mesh node parameter matrix and multidimensional physical boundary constraint tensor constructed in the preceding steps are specifically used as spatially variable parameter terms of the above governing equations, and the underlying coupling solution is performed according to the following logic: based on the node index of the three-dimensional heterogeneous spatial domain, the locally obtained thermal conductivity attenuation parameter is analyzed. The effective thermal conductivity of the corresponding mesh nodes is directly substituted into the conduction terms of the three-dimensional heat conduction equation to characterize the anisotropic heat transfer degradation caused by thermodynamically anomalous boundaries. The analytically obtained effective reaction surface area reduction parameter is then applied. Substituting the source term of the electrochemical kinetics constructed based on the Butler-Volmer law into the charge conservation equation, the effective electrochemical reaction activity at this node is corrected in real time, thereby characterizing the loss of the chemical reaction interface caused by mechanical strain or microstructure damage. Simultaneously, the values ​​of charge transfer resistance and double-layer capacitance mapped by the high-frequency polarization component, and the solid-phase diffusion time constant and concentration limit value mapped by the low-frequency acoustic-electric component are extracted. The initial grid node parameter matrix of the model is overloaded, driving the simulation engine to perform initial steady-state iterative solutions under the current ambient temperature and open-circuit voltage boundary conditions, calculating the initial state variables such as potential, ion concentration, and temperature distribution of the global grid nodes in equilibrium. This simulation object, containing a complete physical parameter matrix and an initial state variable field, is defined as a virtual battery entity.

[0057] Furthermore, the process of outputting the virtual hotspot diffusion trajectory and capacity decay evolution curve includes: analyzing the limiting peak frequency modulated load spectrum and extracting the time-domain transient current boundary excitation sequence; loading the time-domain transient current boundary excitation sequence onto the tab topology boundary node of the virtual battery entity, triggering the transient multiphysics coupling solver to perform adaptive time-step partial differential integral operation; within each time step of the adaptive time-step partial differential integral operation, traversing and extracting the three-dimensional spatial coordinates corresponding to the absolute extreme values ​​of the internal temperature scalar field of the virtual battery entity, and splicing the three-dimensional spatial coordinates in the time evolution order to construct the virtual hotspot diffusion trajectory; within each time step of the adaptive time-step partial differential integral operation, performing volume integral calculation on the global active lithium ion concentration scalar field of the virtual battery entity, extracting the global effective integrated capacity value, and fitting the global effective integrated capacity value sequence along the time axis to generate the capacity decay evolution curve.

[0058] First, the input extreme peak frequency-modulated load spectrum is analyzed and serialized. This extreme peak frequency-modulated load spectrum is a test boundary set based on the nominal capacity of the battery under test. The current rate range it includes is set to be greater than or equal to the battery's rated maximum discharge rate, and the direction and amplitude of the charging and discharging current are switched according to a preset deterministic pseudo-random jump sequence. This pseudo-random jump sequence is pre-generated by a specified pseudo-random number generator (such as a Mason slew algorithm with a fixed random seed) to ensure the determinism and repeatability of the load spectrum in each digital twin simulation. In one specific embodiment, the frequency band is divided into a low-frequency high-current (e.g., 3C to 5C discharge) range of 0.1Hz to 1Hz and a high-frequency low-current (e.g., 1C charge / discharge) range of 1Hz to 10Hz, used to simulate the temperature rise and polarization evolution trend of the battery under the worst operating conditions. A preset load test condition file containing different frequency bands and extreme current amplitudes is read, and the current change node data and corresponding timestamps are extracted. By using piecewise linear interpolation or spline interpolation algorithms, discrete node data is resampled into a continuous discrete data sequence with a fixed minimum time resolution, thereby constructing a time-domain transient current boundary excitation sequence. This provides a continuous external excitation input for solving transient multiphysics. The fixed minimum time resolution is determined based on the Nyquist sampling theorem: the highest frequency component in the limiting peak frequency modulated load spectrum is pre-extracted, and the fixed minimum time resolution is strictly set to be less than or equal to half the period corresponding to the highest frequency component. This ensures that the resampled boundary excitation sequence can retain all high-frequency transient impact characteristics of the original load spectrum without aliasing.

[0059] Next, the boundary loading and transient derivation of the multiphysics solver are executed. The set of topological boundary nodes for the positive and negative tabs in the 3D mesh of the virtual battery entity is identified. The constructed time-domain transient current boundary excitation sequence is dynamically loaded onto the aforementioned tab node set as the Newman boundary condition (i.e., current flux boundary) over time. The transient multiphysics coupled solver is started, and the adaptive step-size control algorithm is activated. During the derivation process, the solver monitors the rate of change of local partial derivatives of the temperature and potential fields within the mesh in real time. When physical quantities undergo drastic changes, the time step is automatically reduced to ensure the convergence of the nonlinear partial differential equation system. When the evolution of physical quantities is stable, the time step is automatically increased to advance the integral solution throughout the entire load cycle. Specifically, the step-size update logic of the adaptive time-step control algorithm is as follows: after each current time step is solved, the relative truncation error of the global temperature and potential fields within the mesh relative to the previous time step is calculated; a preset relative error tolerance threshold (e.g., 10) is set. -4The relative truncation error is compared with a preset error tolerance threshold. If the relative truncation error is greater than the tolerance threshold, it is determined that the physical quantity has undergone a drastic change. The current calculation result is discarded, and the current time step is multiplied by a preset reduction attenuation coefficient (with a value range of 0.1 to 0.5) before iteratively solving again. If the relative truncation error is less than the tolerance threshold, it is determined that the physical quantity has evolved steadily. The current calculation result is accepted, and the step size of the next time step is multiplied by a preset amplification gain coefficient (with a value range of 1.2 to 2.0) to dynamically adjust the integration step size.

[0060] To meet the real-time cycle time requirements of industrial production line sorting, a model order reduction module based on intrinsic orthogonal decomposition was introduced before executing the transient multiphysics coupled solver operation. The specific process is as follows: First, under typical working conditions, the steady-state and partial transient solutions of state variables such as temperature, concentration and electric potential in the three-dimensional heterogeneous spatial domain are extracted, and a snapshot matrix is ​​constructed. To address the issue of inconsistent node counts in heterogeneous meshes, the system performs spatial reference alignment before constructing the snapshot matrix: the solution results of all heterogeneous meshes are uniformly projected onto a set of preset reference reference meshes (usually using the highest density microgrid as the reference), and missing nodes are filled in using an interpolation algorithm, thereby constructing a snapshot matrix with consistent dimensions; then, singular value decomposition is performed on the snapshot matrix to extract the top N orthogonal vectors with energy proportions exceeding a preset threshold (e.g., 99%), which serve as the principal component characteristic basis of the state variables; finally, the Galerkin projection method is used to project the complex system of partial differential equations onto this low-dimensional subspace, constructing a reduced-order multiphysics model with a small number of degrees of freedom. This reduced-order model runs on a high-performance industrial host computer, significantly reducing the dimensionality and computational power requirements of the matrix solution. At the same time, the final result obtained from the derivation is compiled into hexadecimal machine code containing only channel addresses and toggle instructions. The local sorting and execution device (such as a PLC) does not need to undertake any model calculation tasks, but only needs to parse the machine code to execute the action.

[0061] At the end of each convergence time step in the simulation, the virtual hotspot diffusion trajectory is constructed. The spatial extremum search algorithm is invoked to traverse the temperature scalar field data of the global grid nodes of the virtual battery entity at the current time step. The target grid node whose temperature value reaches the absolute maximum value is compared and locked. The three-dimensional spatial coordinates (X, Y, Z axis coordinate values) of the target grid node in the global coordinate system are extracted and bound to the specific simulation timestamp of the current time step. As the time step continues to advance, the spatial coordinates of the absolute extremum recorded at each time step are spliced ​​together in chronological order to generate a discrete geometric path that records the dynamic transfer of the temperature extremum point in the three-dimensional space inside the battery. This path is then output as the virtual hotspot diffusion trajectory.

[0062] Synchronously, at the end of each convergence time step, the capacity decay evolution curve is extracted and fitted. For the current time step, the global active lithium-ion concentration scalar field within the active material region of the virtual battery entity's positive and negative electrodes is extracted. Using the volume of each grid cell as an integral element, the local lithium-ion concentration within each grid cell is multiplied by its corresponding cell volume, and then summed over the entire active phase domain to complete the volume integral calculation and obtain the total global available active lithium reserve at the current moment. Based on the mapping relationship between the Faraday constant and the theoretical specific capacity of the electrode material, this total lithium reserve is converted into a macroscopic global effective integral capacity value. Specifically, the physical equation for the conversion calculation is as follows: in, This represents the global effective integration capacity value obtained from the transformation. Denotes Faraday's constant. This represents the global effective volume integral domain of the positive and negative electrode active materials. This represents the scalar field of local lithium-ion concentration in three-dimensional space within each grid cell at the current time step.

[0063] The effective integral capacity values ​​extracted at each time step throughout the entire simulation period are arranged along the time axis to form the original discrete sequence. A moving average filter is then applied to remove high-frequency numerical fluctuations in the sequence, ultimately generating a smooth capacity decay evolution curve that reflects the irreversible decline in the battery's usable capacity over time.

[0064] Furthermore, the process of sending data to the local sorting execution device involves: performing a first-order time-domain differential operation on the virtual hotspot diffusion trajectory to extract the simulation timestamp when the temperature rise slope first reaches a preset safety limit value, defining it as the thermal runaway critical time; calculating the second derivative distribution of the capacity decay evolution curve, identifying the peak point with the largest absolute value in the second derivative sequence, and defining the evolution cycle number corresponding to the peak point as the aging inflection point cycle number; substituting the thermal runaway critical time and the aging inflection point cycle number as input variables into a preset tiered utilization classification logic table to retrieve and match the corresponding battery health level identifier; converting the battery health level identifier into hexadecimal machine code containing sorting fork action instructions, defining it as the tiered degradation matching code; and transmitting the tiered degradation matching code to the logic controller of the local sorting execution device via a communication bus to drive the execution mechanism to perform sorting actions according to the physical channel corresponding to the tiered degradation matching code.

[0065] First, for the extracted virtual hotspot diffusion trajectory, since this trajectory data is represented as a discrete temperature extreme value sequence on the time axis, a discrete-time difference algorithm is used to perform a first-order time-domain differential operation on it. Specifically, the discrete sequence is traversed, the temperature extreme value difference between two adjacent time steps is calculated, and it is divided by the corresponding time step size to generate a temperature rise slope sequence representing the rate of temperature change. Subsequently, a preset safety limit value is introduced. This limit value is a critical temperature rise rate threshold calibrated by conducting destructive thermal runaway experiments on the same type of battery in advance. The values ​​in the temperature rise slope sequence are compared with this safety limit value one by one in chronological order. The simulation timestamp corresponding to the first time when the temperature rise slope is greater than or equal to the safety limit value is extracted, and this timestamp is established as the critical time of thermal runaway of the battery under test, quantifying the battery's safety margin.

[0066] Next, for the output capacity decay evolution curve, a smoothing filtering algorithm is first used to denoise the original capacity sequence. Specifically, the smoothing filtering algorithm adopts the Savitzky-Golay polynomial smoothing filtering algorithm (SG filtering), setting the sliding window length to a preset odd number of cycles (such as 5 to 11 cycles), and the polynomial fitting order to be second or third order. This is to filter out high-frequency random measurement noise while retaining the true local shape features and inflection point positions of the capacity decay curve without phase shift. Subsequently, two discrete difference operations are performed on the smoothed capacity sequence to calculate the second derivative values ​​corresponding to each evolution cycle, forming a second derivative distribution sequence. This sequence reflects the acceleration of the battery capacity decay rate. The extreme value search logic is invoked to traverse the entire second derivative distribution sequence to find and lock the value point with the largest absolute value. The number of evolution cycles corresponding to this largest absolute value point is extracted and defined as the aging inflection point cycle number, which characterizes the expected lifespan node where a malignant change occurs in the battery's internal aging mechanism.

[0067] After obtaining the two key evaluation parameters mentioned above, a health level matching based on a two-dimensional lookup table is performed to construct a preset tiered utilization classification logic table. This table uses multiple time threshold intervals of thermal runaway critical time as the horizontal evaluation dimension and multiple cycle number threshold intervals of aging inflection point cycle number as the vertical evaluation dimension. Each intersecting grid in the table has a pre-fixed corresponding battery health level identifier (e.g., tiered utilization level A, level B, or scrap level). In a specific embodiment, for a certain type of lithium iron phosphate battery, the preset safety reference time is, for example, 600 seconds, and the preset cycle life lower limit is, for example, 1500 cycles. The threshold cross-judgment rules of the logic table are configured as follows: if the thermal runaway critical time is greater than the preset safety reference time (i.e., >600s) and the number of aging inflection point cycles is greater than the preset lower limit of cycle life (i.e., >1500 times), then the corresponding cross grid is solidified as Grade A for secondary use (e.g., it can be used in energy storage power stations); if the thermal runaway critical time is greater than the safety reference time (i.e., >600s), but the number of aging inflection point cycles is less than or equal to the lower limit of cycle life (i.e., ≤1500 times), then the corresponding cross grid is solidified as Grade B for secondary use (e.g., it can be used in low-speed electric vehicles); if the thermal runaway critical time is less than or equal to the safety reference time (i.e., ≤600s), then regardless of the number of aging inflection point cycles, the cross grid it falls into is forcibly solidified as scrap level (limited to physical dismantling after safe discharge). The currently calculated thermal runaway critical time and the number of aging inflection point cycles are used as joint index variables and input into the logic table to determine the cross interval into which these two variables fall, and retrieve the battery health level identifier corresponding to the interval.

[0068] Finally, the encapsulation of control commands and the issuance of physical sorting actions involve calling the pre-built hardware control protocol mapping library within the system to translate the retrieved battery health level identifier into a control message recognizable by the underlying hardware. This message is encapsulated into a standard hexadecimal machine code data frame. The data frame strictly includes a start bit, the physical sorting channel address corresponding to the health level, the push-pull action command code for the sorting fork, and a cyclic redundancy check code for verification. This complete data frame is defined as the tiered degradation matching code. After construction, this hexadecimal tiered degradation matching code is sent in real-time to the programmable logic controller (PLC) on the physical sorting line via an industrial field communication bus (such as CAN bus or Modbus bus). The PLC decodes and verifies the received machine code and outputs the corresponding level signal to the corresponding pneumatic or electric actuator to drive the sorting fork to push the battery under inspection into the physical tiered utilization track that matches its health level.

[0069] Example 2: At the automated testing station, the system first controls a broadband ultrasonic transducer to send a sweep frequency signal to the battery under test, extracting its mechanical intrinsic frequency characteristic value under the current aging state. Subsequently, this characteristic value is synchronously configured as the resonant fundamental frequency of the electrical and acoustic excitation circuits. Through hardware digital phase-locked loop dynamic compensation for transmission delay, the electrical high-frequency polarization pulse and ultrasonic pulse are forced to maintain zero phase difference and synchronous injection at the physical injection boundary outside the battery, thereby eliminating the initial system error of spatial lag calculation. In this state, the high-speed data acquisition card acquires the original voltage response sequence at the tab end. After digital quadrature demodulation and band-stop filtering, the acoustic-electric intermodulation harmonic components reflecting the complex interface state inside the battery are extracted. At the same time, the system synchronously reads the array temperature measurement node data attached to the outer shell surface to reconstruct the surface thermal gradient matrix and captures the ultrasonic echo waveform to calculate the ultrasonic penetration attenuation time delay sequence.

[0070] The system calls the preset orthogonal wavelet basis functions to perform multi-scale discrete wavelet decomposition on the extracted acoustic-electric intermodulation nonlinear distortion signal. The detailed coefficient sequence of the highest frequency band is extracted and the wavelet energy spectral density within the sliding time-domain window is calculated and defined as the high-frequency polarization component. The approximate coefficient sequence of the lowest frequency band is extracted, and after edge removal and smoothing filtering, the absolute convex point is locked and defined as the low-frequency acoustic-electric component. Subsequently, the system loads the preset three-dimensional finite element mesh corresponding to this battery model. The surface thermal gradient matrix is ​​rigidly aligned to the surface of the outer shell mesh through bilinear coordinate interpolation to generate thermodynamic anomaly boundary weights. The ultrasonic attenuation delay is projected along the transmission path to the interior of the bulk mesh to generate the mechanical porosity degradation factor. The two are then spliced ​​along the mesh topological connectivity matrix dimension to construct a multi-dimensional physical boundary constraint tensor.

[0071] The system performs discrete difference operations along the three orthogonal coordinate axes of the multidimensional physical boundary constraint tensor to synthesize a global three-dimensional spatial gradient matrix. The system dynamically calculates the statistical mean and standard deviation of the current gradient matrix to generate an adaptive global gradient threshold. Through comparison, the system identifies a region within the battery electrode roll where the absolute gradient value exceeds the limit, indicating a microstructural anomaly (such as local lithium plating or active material shedding). The system then triggers a mesh subdivision command in the anomalous region, adaptively dividing it into a high-node-density micro-solution mesh based on the local range smoothing criterion; simultaneously, it performs mesh coarsening on stable volumetric regions where the gradient does not exceed the limit. For dangling nodes at the macro-micro mesh boundary, the area interpolation shape function of the macro-surface is extracted, and multi-point constraint equations are established to eliminate redundant degrees of freedom, ultimately stitching together a closed three-dimensional heterogeneous spatial domain.

[0072] The system performs discrete difference operations along the three orthogonal coordinate axes of the multidimensional physical boundary constraint tensor to synthesize a global three-dimensional spatial gradient matrix. The system dynamically calculates the statistical mean and standard deviation of the current gradient matrix to generate an adaptive global gradient threshold. Through comparison, the system identifies a region within the battery electrode roll where the absolute gradient value exceeds the limit, indicating a microstructural anomaly (such as local lithium plating or active material shedding). The system then triggers a mesh subdivision command in the anomalous region, adaptively dividing it into a high-node-density micro-solution mesh based on the local range smoothing criterion; simultaneously, it performs mesh coarsening on stable volumetric regions where the gradient does not exceed the limit. For dangling nodes at the macro-micro mesh boundary, the area interpolation shape function of the macro-surface is extracted, and multi-point constraint equations are established to eliminate redundant degrees of freedom, ultimately stitching together a closed three-dimensional heterogeneous spatial domain.

[0073] The system performs a first-order time-domain differential on the virtual hotspot propagation trajectory, compares it with the preset safe temperature rise rate limit, and extracts the simulation timestamp of the first exceeding the limit as the critical time for thermal runaway. Simultaneously, it uses polynomial smoothing filtering to denoise the capacity decay curve, and finds the maximum peak point through two consecutive discrete differences to determine the number of aging inflection point cycles. The system inputs these two key parameters into a tiered utilization hierarchical logic table. Because the critical time for thermal runaway of the abnormal battery is shorter than the safety benchmark threshold set in the logic table, the system triggers a rejection judgment and retrieves the scrap-level health identifier. Finally, the system converts this identifier into a standard hexadecimal machine code data frame containing a fixed frame header, scrap channel address, fork action command, and cyclic redundancy check code. This data frame is then sent to the programmable logic controller (PLC) on the sorting line via the industrial bus, driving the pneumatic actuator to precisely push the battery under inspection into the scrap dismantling physical channel, completing the automated sorting closed loop.

[0074] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A method for real-time evolution of electrochemical characteristics of a battery based on a digital twin model, characterized in that, include: High-frequency polarization pulses and ultrasonic pulses of the same frequency and phase-locked loop are applied to the battery under test, and the acoustic-electric intermodulation nonlinear distortion signal, surface thermal gradient matrix and ultrasonic penetration attenuation time delay sequence are acquired simultaneously. Wavelet decomposition was applied to process the nonlinear distortion signal of acoustic-electric intermodulation, separating the high-frequency polarization component and the low-frequency acoustic-electric component. Using a pre-defined three-dimensional finite element mesh as a reference, the surface thermal gradient matrix and the ultrasonic penetration attenuation delay sequence are projected onto the mesh nodes to construct a multi-dimensional physical boundary constraint tensor. The spatial gradient distribution of the multidimensional physical boundary constraint tensor is calculated, a microscopic solution mesh is generated in the structural anomaly region where the gradient exceeds the limit, and a three-dimensional heterogeneous spatial domain is reconstructed. Substitute the high-frequency polarization component, the low-frequency acoustic-electric component, and the multi-dimensional physical boundary constraint tensor into the multiphysics field model established on the three-dimensional heterogeneous spatial domain to initialize and generate a virtual battery entity. Inject the extreme peak frequency modulation load spectrum into the virtual battery entity to drive the execution of accelerated time step integral derivation, and output the virtual hot spot diffusion trajectory and capacity decay evolution curve; extract the thermal runaway critical time of the virtual hot spot diffusion trajectory and the aging inflection point cycle number of the capacity decay evolution curve, calculate the tiered degradation matching code based on the two, and send it to the local sorting execution device.

2. The method for real-time evolution of battery electrochemical characteristics based on a digital twin model according to claim 1, characterized in that, The process of acquiring the acoustic-electric intermodulation nonlinear distortion signal, surface thermal gradient matrix, and ultrasonic penetration attenuation delay sequence includes: extracting the mechanical intrinsic frequency characteristic value of the battery under test, and synchronously configuring the mechanical intrinsic frequency characteristic value as the reference resonant frequency of the electrical excitation circuit and the acoustic wave transmission circuit; controlling the output delay of the electrical excitation circuit and the acoustic wave transmission circuit through a hardware phase-locked loop, using the actual voltage phase reaching the electrode terminal and the actual acoustic wave phase reaching the shell surface as feedback control nodes, defining the electrode terminal and the shell surface as the physical injection boundary for initial phase alignment, locking the zero phase difference state at the physical injection boundary, and synchronously injecting high-frequency polarization pulses and ultrasonic pulses into the battery under test; and acquiring the original voltage response sequence of the battery under test in the zero phase difference state. The original voltage response sequence is input into the demodulation operator, and demodulation is performed using the reference resonant frequency as the reference carrier. The fundamental frequency response component with the same reference resonant frequency as the original voltage response sequence is filtered out, and the harmonic components modulated by the acoustic-electric coupling effect are extracted. The harmonic components are defined as acoustic-electric intermodulation nonlinear distortion signals. The values ​​of the array temperature measurement nodes arranged on the surface of the battery under test are read under the Joule heat conduction effect generated by the basic DC diagnostic load. The temperature change rate of adjacent nodes is calculated by applying the spatial interpolation function, and the surface thermal gradient matrix is ​​reconstructed. The ultrasonic echo waveform penetrating the battery under test is captured. The time difference vector is calculated by subtracting the high-frequency polarization pulse excitation time from the arrival time of the ultrasonic echo waveform envelope peak. The time difference vector is defined as the ultrasonic penetration attenuation time delay sequence.

3. The method for real-time evolution of battery electrochemical characteristics based on a digital twin model according to claim 1, characterized in that, The process of separating the high-frequency polarization component and the low-frequency acoustic-electric component includes: introducing a preset orthogonal wavelet basis function, performing multi-scale discrete wavelet transform decomposition on the acoustic-electric intermodulation nonlinear distortion signal to generate a wavelet signal matrix containing multiple frequency band levels; extracting the highest frequency band detail coefficient sequence from the wavelet signal matrix and defining the highest frequency band detail coefficient sequence as the high-frequency polarization component; extracting the lowest frequency band approximation coefficient sequence from the wavelet signal matrix and defining the lowest frequency band approximation coefficient sequence as the low-frequency acoustic-electric component; calculating the wavelet energy spectral density feature value of the high-frequency polarization component within the time domain window, extracting the waveform envelope peak feature value of the low-frequency acoustic-electric component, and concatenating the wavelet energy spectral density feature value and the waveform envelope peak feature value to construct an evolution feature vector.

4. The method for real-time evolution of battery electrochemical characteristics based on a digital twin model according to claim 1, characterized in that, The process of constructing the multidimensional physical boundary constraint tensor includes: extracting the spatial node coordinate system of a preset three-dimensional finite element mesh, and separating the surface contour node set and the internal volume phase node set from the spatial node coordinate system; aligning the surface thermal gradient matrix to the surface contour node set through coordinate interpolation transformation, and assigning values ​​to the surface contour node set to generate thermodynamic anomaly boundary weights; projecting the ultrasonic penetration attenuation delay sequence to the internal volume phase node set along the ultrasonic propagation and transmission path, calculating the local attenuation coefficient of the corresponding node and assigning values ​​to generate mechanical porosity degradation factors; extracting the topological adjacency matrix of the spatial node coordinate system, and concatenating the thermodynamic anomaly boundary weights and the mechanical porosity degradation factors along the dimensions of the topological adjacency matrix to construct the multidimensional physical boundary constraint tensor.

5. The method for real-time evolution of battery electrochemical characteristics based on a digital twin model according to claim 1, characterized in that, The process of reconstructing a three-dimensional heterogeneous spatial domain includes: performing differential operations along the spatial coordinate axes of the multidimensional physical boundary constraint tensor, splicing partial derivative values ​​to generate a three-dimensional spatial gradient matrix; comparing the values ​​of the three-dimensional spatial gradient matrix with a preset gradient threshold, marking coordinate intervals with values ​​greater than the preset gradient threshold as structurally abnormal regions, and marking coordinate intervals with values ​​less than the preset gradient threshold as volumetric regions; triggering a mesh subdivision command in the structurally abnormal regions, cutting the preset three-dimensional finite element mesh into a high-node-density micro-solution mesh according to a preset geometric ratio; triggering a mesh coarsening command in the volumetric regions, merging the preset three-dimensional finite element mesh according to topological connectivity to generate a low-node-density volumetric macro-mesh; extracting the suspended boundary nodes between the micro-solution mesh and the volumetric macro-mesh, aligning the suspended boundary nodes through spatial coordinate interpolation, and stitching the boundary closure to generate a three-dimensional heterogeneous spatial domain.

6. The method for real-time evolution of battery electrochemical characteristics based on a digital twin model according to claim 1, characterized in that, The initialization process for generating a virtual battery entity includes: extracting high-frequency polarization components and mapping and updating the charge transfer internal resistance reference value and interface double-layer capacitance value of each grid node in the three-dimensional heterogeneous spatial domain; extracting low-frequency acoustic-electric components and mapping and updating the solid-phase diffusion time constant and maximum lithium insertion / extraction concentration limit value of each grid node in the three-dimensional heterogeneous spatial domain; analyzing the multidimensional physical boundary constraint tensor and converting the thermodynamic anomaly boundary weights into the local thermal conductivity attenuation parameter of the corresponding grid node; converting the mechanical porosity degradation factor in the multidimensional physical boundary constraint tensor into the effective reaction surface area reduction parameter of the corresponding grid node; extracting the charge transfer internal resistance reference value, interface double-layer capacitance value, solid-phase diffusion time constant, maximum lithium insertion / extraction concentration limit value, local thermal conductivity attenuation parameter, and effective reaction surface area reduction parameter, and reloading the initial grid node parameter matrix of the multiphysics model; driving the multiphysics model loaded with the initial grid node parameter matrix to perform initial steady-state solution, constructing a multiphysics simulation instance object containing the initial state variables of all nodes in the three-dimensional heterogeneous spatial domain, and defining the multiphysics simulation instance object as a virtual battery entity.

7. The method for real-time evolution of battery electrochemical characteristics based on a digital twin model according to claim 1, characterized in that, The process of outputting the virtual hotspot diffusion trajectory and capacity decay evolution curve includes: analyzing the limiting peak frequency modulated load spectrum and extracting the time-domain transient current boundary excitation sequence; loading the time-domain transient current boundary excitation sequence onto the tab topology boundary node of the virtual battery entity, triggering the transient multiphysics coupling solver to perform adaptive time-step partial differential integral operation; within each time step of the adaptive time-step partial differential integral operation, traversing and extracting the three-dimensional spatial coordinates corresponding to the absolute extreme values ​​of the internal temperature scalar field of the virtual battery entity, and splicing the three-dimensional spatial coordinates in the time evolution order to construct the virtual hotspot diffusion trajectory; within each time step of the adaptive time-step partial differential integral operation, performing volume integral calculation on the global active lithium ion concentration scalar field of the virtual battery entity, extracting the global effective integrated capacity value, and fitting the global effective integrated capacity value sequence along the time axis to generate the capacity decay evolution curve.

8. The method for real-time evolution of battery electrochemical characteristics based on a digital twin model according to claim 1, characterized in that, The process of sending the data to the local sorting execution device includes: performing a first-order time-domain differential operation on the virtual hotspot diffusion trajectory, extracting the simulation timestamp when the temperature rise slope first reaches the preset safety limit value, and defining it as the thermal runaway critical time; calculating the second derivative distribution of the capacity decay evolution curve, identifying the peak point with the largest absolute value in the second derivative sequence, and defining the evolution cycle number corresponding to the peak point as the aging inflection point cycle number; substituting the thermal runaway critical time and the aging inflection point cycle number as input variables into a preset tiered utilization classification logic table, and retrieving and matching the corresponding battery health level identifier; converting the battery health level identifier into hexadecimal machine code containing sorting fork action instructions, and defining it as the tiered degradation matching code; transmitting the tiered degradation matching code to the logic controller of the local sorting execution device through the communication bus, and driving the execution mechanism to perform sorting actions according to the physical channel corresponding to the tiered degradation matching code.

Images