Three-parameter coupled grading early warning method for thermal runaway of energy storage battery

Abstract

The invention discloses a three-parameter coupled grading early warning method for thermal runaway of an energy storage battery, the method is operated based on a monitoring system, the monitoring system comprises a sampling module, an analysis module and a judgment module, and a battery state model is established through coupling of three physical quantities of stress, target component concentration and temperature to describe a battery state. After data collected by the sampling module is transmitted to the analysis module, the analysis module predicts future moment values of three physical quantities based on a battery state model, then a final prediction result is obtained through a PF algorithm, then a risk index is obtained, and then the judgment module judges the current state of the battery according to the risk index. According to the method, the battery state can be described more accurately, the reliability of the prediction result is improved, and the result obtained by the model is more accurate; all prediction results are processed through a PF algorithm, so that the finally obtained battery state is closer to the real situation; the battery state is reflected through the risk index, and a reliable basis can be provided for judgment of the battery state.

Description

Technical Field

[0001] This invention belongs to the field of lithium-ion battery safety monitoring technology, specifically relating to a graded early warning method for thermal runaway of energy storage batteries based on the coupling of stress-gas-temperature three parameters. Background Technology

[0002] With the widespread application of lithium-ion batteries in new energy vehicles, energy storage devices, and other fields, their safety has become an increasingly important concern. Under extreme conditions such as high temperature, overcharging, and external impact, lithium-ion batteries may experience thermal runaway, leading to fires or even explosions, causing serious property damage and personal injury. Existing thermal runaway early warning technologies mainly rely on monitoring a single parameter, such as using temperature sensors to detect the temperature rise of the battery surface, strain gauges to sense the expansion and deformation of the battery, or gas sensors to detect the concentration of flammable gases released from side reactions inside the battery. These monitoring methods can achieve certain results in relatively simple scenarios, but they have significant limitations under complex conditions. First, the monitored parameters may have a response lag. Second, single parameters are easily affected by environmental interference or local defects. For example, the battery surface temperature usually only rises after thermal runaway has occurred, and the battery surface temperature may also rise under prolonged sunlight exposure. Therefore, monitoring a single parameter may result in false alarms or missed alarms. Furthermore, the change of a single parameter over time is insufficient to reflect the complex physical-chemical evolution process inside the battery, making it difficult to detect risks in the early stages of thermal runaway. The lack of information about the physical-chemical evolution process inside the battery also poses significant challenges to subsequent accident cause analysis.

[0003] With the development of multiphysics coupling theory in recent years, researchers in this field have gradually recognized that thermal runaway in batteries is a holistic process involving stress changes, gas release, and heat accumulation. A single parameter is insufficient to reflect the true state of the battery, and multi-parameter synchronous sensing and fusion modeling has become a more forward-looking early warning strategy. However, research on multi-parameter fusion early warning methods is still in its early stages, facing the following challenges: first, the complex nonlinear coupling relationships between different parameters make modeling the internal physical-chemical evolution process of the battery difficult; second, multi-source data fusion algorithms are not yet mature, making it difficult to balance the real-time performance and stability of the results; and third, there is no unified risk warning standard. Therefore, it is necessary to construct a thermal runaway early warning method based on a multi-parameter coupling mechanism, possessing stable and efficient multi-level risk identification capabilities, in order to more accurately predict abnormal states within the battery. Summary of the Invention

[0004] The purpose of this invention is to provide a three-parameter coupled method for graded early warning of thermal runaway in energy storage batteries. This method establishes a model that describes the thermal runaway process of the battery using three physical quantities: stress, target component concentration, and temperature. In this model, the three physical quantities are coupled. Then, the predicted values ​​of the three physical quantities at future times are obtained through the model, and the risk index J(t) used to assess the current state of the battery is derived based on the predicted values.

[0005] The technical solution adopted in this invention is as follows:

[0006] A three-parameter coupled method for graded early warning of thermal runaway in energy storage batteries, specifically including the following steps:

[0007] Step 1. Build a monitoring system. The monitoring system includes a sampling module, an analysis module, and a judgment module. The sampling module and the analysis module are connected by a signal, and the analysis module and the judgment module are also connected by a signal. The sampling module includes three types of sensors: a stress sensor, a gas concentration sensor, and a temperature sensor. All sensors sample simultaneously and the time interval between two adjacent samplings is equal.

[0008] Step 2. Establish a battery state model, including the following steps:

[0009] Step 2.1. Establish a model for the time-varying nature of three physical quantities related to the battery: stress, target component concentration, and temperature;

[0010] Step 2.2. Determine the coupling relationship between the above three physical quantities in the physical-chemical evolution process of the battery;

[0011] Step 3. The analysis module predicts the future values ​​of the three physical quantities, including the following steps:

[0012] Step 3.1. The data measured by all sensors in the sampling module are transmitted to the analysis module in real time. The analysis module substitutes the received data into the battery state model to obtain the predicted values ​​of three physical quantities. There are multiple predicted values ​​for each physical quantity.

[0013] Step 3.2. Construct multiple sample vectors based on the predicted values. Then, the analysis module processes all sample vectors using the particle filtering algorithm, i.e., the PF algorithm, to obtain a prediction vector. The prediction vector contains the final predicted values ​​of the three physical quantities.

[0014] Step 4. Based on the final predicted values ​​of the three physical quantities, the analysis module derives the current risk index J(t) of the battery.

[0015] Step 5. The analysis module transmits the J(t) value to the judgment module in real time. The judgment module makes a judgment on the current state of the battery based on the J(t) value and outputs the judgment result.

[0016] In existing technologies, battery thermal runaway processes are mostly described based on a single physical quantity. This makes it difficult to accurately describe the internal physicochemical evolution of the battery, and thus impossible to make accurate predictions about the battery state. The method of this invention couples three physical quantities: stress, target component concentration, and temperature, providing a more comprehensive and accurate description of the battery thermal runaway process. Besides causing temperature increases, battery thermal runaway can also induce side reactions such as electrolyte volatilization and SEI film rupture, releasing gases such as CO2 and CH4. Simultaneously, thermal runaway causes battery expansion, and due to electrolyte decomposition and lithium plating, the battery casing will experience significant volumetric stress. Therefore, stress, target component concentration, and temperature are chosen as monitoring indicators. To minimize the deviation between predicted and actual results, this invention uses a battery thermal runaway process model to derive multiple prediction results, and then processes all prediction results using the PF algorithm to make the final result closer to the actual result. Based on the final result, the current battery risk index J(t) is derived, and combined with pre-set judgment criteria, the current state of the battery is judged, thus standardizing the risk assessment of battery thermal runaway.

[0017] Further optimization, Step 2.1 specifically includes the following steps:

[0018] Step 2.1.1. Establish a model of temperature change over time, which is expressed by the following formula.

[0019]

[0020] Q ohmic (t)=I 2 (t)R int (T)

[0021]

[0022] Q gas (t)=ΔH g ·f g (t)

[0023] Q loss (t)=h·A·[T(t)-T env ]

[0024] Where C p Q represents the specific heat capacity of the battery casing material. ohmic (t) represents ohmic heat, Q side (t) represents the heat released by the side reaction, Q. gas(t) represents the heat released by the gas, Q loss (t) represents the heat dissipation; I(t) represents the current inside the battery at time t; R int (T) represents the battery's internal resistance, ΔH rxn For the heat of reaction, k rxn E is the reaction coefficient. a Let R be the gas internal energy, γ1 be the gas constant, and ΔH be the gas constant. g Let h be the gas friction coefficient, h be the thermal conductivity coefficient of the battery casing, A be the surface area of ​​the battery casing, and T be the total surface area of ​​the battery casing. env Ambient temperature;

[0025] Step 2.1.2. Establish a model for the change of target component concentration over time. This model is expressed by the following formula.

[0026] C g (t)+λ g ∫C g (t)dt=∫f g (t)dt+C0 (3)

[0027]

[0028] Where t is time, C g (t) represents the concentration of the target component at time t, C0 represents the initial concentration of the target component, and λ g f is the gas diffusion coefficient. g (t) represents the gas release rate, where k1, k2, k3, β1, and β2 are constants;

[0029] The target component concentration C is obtained from formula (3). g (t) represents the rate of change at time t, which is given by the following formula.

[0030]

[0031] Step 2.1.3. Establish a model for the stress variation of the battery casing over time. This model is expressed by the following formula.

[0032]

[0033] Where σ(t) is the stress on the battery casing, E is the elastic modulus of the battery casing material, L0 is the original length of the monitored portion of the battery casing, ΔL(t) is the elongation of the monitored portion of the battery casing, and α T Let α be the temperature-induced strain coefficient, T0 be the initial temperature of the battery, T(t) be the battery temperature, and α be the temperature-induced strain coefficient. g The strain coupling coefficient caused by gas expansion;

[0034] Based on formula (1), the rate of change of stress σ(t) at time t is obtained, i.e. It can be expressed as the following formula

[0035]

[0036] Where ε(t) is the strain rate, The value is obtained directly from the sensor;

[0037] The model constructed in this step is not a model that describes the battery state using a single physical quantity, but rather a model that describes the physical quantity based on other physical quantities. In these models, there are linear relationships between different physical quantities, and there are also linear relationships between the derivatives of different physical quantities. This is part of the three-parameter coupling relationship.

[0038] Further optimization, Step 2.2 specifically includes the following steps:

[0039] Step 2.2.1. Construct a low-order vector X(t) and a high-order vector X(t) each containing five elements. 2 The two vectors are represented by the following formulas.

[0040]

[0041] in Let σ(t) be the first derivative of σ(t) with respect to time t, that is, the rate of change of σ(t) at time t. Let T(t) be the first derivative of T(t) with respect to time t, that is, the rate of change of T(t) at time t.

[0042] Step 2.2.2. Express the first derivative of the low-order vector X(t) with respect to time t as follows:

[0043]

[0044] Where A is the state coupling matrix, B is the input control matrix, G is the generalized coupling matrix, u(t) is the external disturbance vector, and w(t) is the system noise; let b be the number of elements in u(t), where b is a positive integer, then matrix B is a 5xb matrix. Matrices A and G are expressed in the following forms.

[0045]

[0046] Where a ij With g ij Let i and j be elements in matrices A and G, respectively, where i and j represent the row and column numbers of the element in their respective matrices, i,j∈{1,2,3,4,5}. All elements in A and G are constants.

[0047] Step 1.2.3. From formula (8), the derivative of X(t) with respect to time t is:

[0048]

[0049] in C g The first derivative of (t) with respect to time t, i.e., C g (t) is the rate of change of time t. Let σ(t) be the second derivative with respect to time t. Let T(t) be the second derivative of T(t) with respect to time t;

[0050] Substituting the matrix determinants in formulas (8), (9), (11), (12), and (13) into formula (10), we obtain the following relation.

[0051]

[0052] Where μ1, μ2, μ3, μ4, and μ5 are all error terms, and the error terms are jointly determined by matrix B, u(t), and w(t), δ i For formula (10)G·X(t) 2 The element in the i-th row of the determinant.

[0053] In the coupled model, the three physical quantities are not only linearly related, but also nonlinearly related. The nonlinear relationship is obtained by constructing a relational expression through higher-order equations. As can be seen from formula (14), the quadratic term, the first derivative term, or the product of two physical quantities of a certain physical quantity are linearly related to the second derivative term of a certain physical quantity. By describing this linear relationship through a linear relational expression, a nonlinear relational expression between a certain physical quantity and other physical quantities can be constructed.

[0054] Further optimization involves the following steps in Step 3.1 to derive the predicted value:

[0055] Step 3.1.1. In formula (14), the terms on the left side of the equal sign are called lower-order terms, and the terms on the right side of the equal sign are called higher-order terms. All lower-order terms in any equation are considered to be linearly related to the higher-order terms in the equation. Select multiple lower-order terms and establish a linear relationship between each lower-order term and its corresponding higher-order term.

[0056] Step 3.1.2. The derivative in the linear relationship can be expressed using the finite difference method as follows:

[0057]

[0058] Where U is the physical quantity to be calculated, U∈{σ,C} g ,T},U t U is the value of this physical quantity at a certain moment. t+1 This represents the value at the next moment, where Δt is the time interval between two consecutive records. This is the first derivative of the physical quantity at a certain moment. The first derivative at the next time step. The second derivative at a certain moment;

[0059] Step 3.1.3. Record experimental data n times during the battery thermal runaway experiment, where n is a positive integer. The experimental data includes the values ​​of all time-related variables in formulas (1) to (7), which are Q. ohmic (t), Q side (t), Q gas (t), Q loss (t), I(t), T(t), ΔL(t) The time interval between two consecutive records is equal; the experimental data are sorted according to the recording time, and the data sequence number is denoted as p, where p is a positive integer and 1≤p≤n. All values ​​in the same record have the same sequence number.

[0060] Step 3.1.4. Substitute the time-related variables in the experimental data into the corresponding formulas in formulas (1) to (7) in order to obtain the stress, target component concentration and temperature at different times, as well as the first derivatives of the three physical quantities with respect to t, and then obtain the values ​​of the low-order and high-order terms in the linear relationship at different times; then fit the coefficients in the linear relationship by the least squares method to determine the linear relationship.

[0061] Step 3.1.5. The stress, target component concentration, and temperature at the current moment are denoted as σ. s C g,s and T s , will σ s C g,s and T s Substituting these values ​​into the linear equation, we obtain the predicted values ​​of stress, target component concentration, and temperature at future times, denoted as σ. fc C g,fc and T fc All three physical quantities have multiple predicted values.

[0062] Further optimization involves the following steps in Step 3.2 to derive the prediction vector:

[0063] Step 3.2.1. Denote the sample vector as S. q,r Where q is the index and r is the iteration number, both q and r are positive integers. Based on the predicted values, m initial sample vectors S are constructed. q,0 S q,0 =[σ fc C g,fc ,T fc ] TEach initial sample vector corresponds to a unique index and 1≤q≤m. Any two S... q,0 All are unequal; at the same time, the iteration threshold D is determined based on the number of initial sample vectors;

[0064] Step 3.2.2. Obtain the observation vector corresponding to the current sample vector using the following formula.

[0065] Z q,r =h(S q,r )+v q,r (16)

[0066] Z q,r Let h be the observation vector, and v be the observation function. q,r For measuring noise;

[0067] Step 3.2.3. Calculate the weight value corresponding to each sample vector using the following formula.

[0068] W q,r =W q,r-1 ·p(Z q,r |S q,r (17)

[0069] Among them W q,r For S q,r The corresponding weight values, where p is a probability function, p(A|B) represents the probability of event A occurring given that event B has occurred, and its value is determined by the observation model itself; initial weight values.

[0070] Step 3.2.4. Generate a new sample vector by substituting the current sample vector into the following formula.

[0071] S q,r =f(S) q,r-1 ,u' q,r )+w' q,r (18)

[0072] Where u' q,r Due to external disturbances, w' q,r This refers to system process noise.

[0073] Compare r in the newly generated sample vector with D. If the condition r = D is met, proceed to Step 3.2.5. If not, return to Step 3.2.2 and use the newly generated sample vector as the current sample vector.

[0074] Step 3.2.5. Obtain the prediction vector F by weighted summation of all sample vectors. F is expressed as the following formula.

[0075]

[0076] In Step 3, generating multiple sets of predicted values ​​and then processing them using the Power Forward (PF) algorithm is crucial because generating only one set of predicted values, i.e., a single sample vector, would compromise the accuracy of that prediction. By taking multiple sample vectors and using the PF algorithm to generate even more sample vectors, and simultaneously calculating the weight value corresponding to each sample vector—the closer the sample vector is to the actual situation, the larger its corresponding weight value—and then generating a certain number of sample vectors before weighted summing of all the sample vectors, the final result will be closer to the actual situation.

[0077] Further optimization is achieved by deriving the battery risk index J(t) at the current moment in Step 4 as follows:

[0078] Step 4.1. Construct the risk function for calculating the risk index J(t). The risk function is expressed as follows:

[0079] J(t) = w σ ·σ'(t)+w g ·C g '(t)+w T ·T'(t) (20)

[0080] Where σ'(t) and C g T'(t) and T'(t) are the predicted values ​​of stress, target component concentration, and temperature at time t, respectively. The predicted value is the ratio of the final predicted value of a physical quantity to its experimentally recorded maximum value. σ w g w T These are the weight values ​​corresponding to the three physical quantities: stress, target component concentration, and temperature.

[0081] Step 4.2. Determine w based on actual working conditions σ w g w T The specific values ​​of the three components, and the three weight values ​​satisfy w σ +w g +w T =1 condition;

[0082] Step 4.3. Based on the final predicted value at the current time, calculate the difference between the current time and the predicted value σ'(t), C. g Substitute the three predicted values, T'(t) and T'(t), into formula (20) to obtain the current time J(t).

[0083] Further optimization is achieved by determining the current battery state in Step 5 as follows:

[0084] Step 5.1. Based on the battery operating conditions, take multiple critical values ​​within the range of [0,1], and divide the range of [0,1] into multiple intervals through these critical values. Each interval corresponds to a warning state.

[0085] Step 5.2. Compare J(t) with the magnitude of all critical values ​​to determine the interval in which J(t) is located, and then determine the warning status at the current moment.

[0086] The J(t) obtained through the above process has a value in the range [0,1]. Thus, the value of J(t) can be used to measure the risk of the battery in its current state. The closer J(t) is to 0, the lower the risk; the closer it is to 1, the higher the risk. Based on this characteristic, the range [0,1] is divided into multiple intervals. The interval with the median value closer to 1 represents a higher risk level. When J(t) falls into an interval with the median value closer to 1, the judgment module should output a higher level of warning.

[0087] Further optimization involves determining the values ​​of constants k1, k2, k3, β1, and β2 in Step 1.1.1 through the following steps:

[0088] I. Set the target component concentration C g The rate of change over time is expressed by the finite difference method, and by combining formulas (2) and (3), we can obtain the relationship between σ and σ. p C g,p and T p formula

[0089]

[0090] II. All C g,p Substituting these values ​​into formula (21) in sequence, we obtain the values ​​on the left side of the equation (21), which are denoted as Ω. p Then, the right side of formula (21) is regarded as an algebraic expression with k1, k2, k3, β1, and β2 as unknowns, and all σ p C g,p and T p Substituting the values ​​into the equation in order, we construct the following fitting function.

[0091]

[0092] Where MSE is the sum of squared residuals, and k1', k2', k3', β1', and β2' are the variables corresponding to k1, k2, k3, β1, and β2, respectively. Changing the values ​​of these variables yields different MSE values. Each MSE value and the k1', k2', k3', β1', and β2' at which that value is obtained are recorded. Then, the minimum MSE value is found, and the k1', k2', k3', β1', and β2' at which the minimum MSE value is obtained are taken as the final determined values ​​of k1, k2, k3, β1, and β2.

[0093] Since formulas describing the concentration of a component in a gas are often not derived precisely from physical principles, but rather based on fitting statistical data, the accuracy of the fitting directly affects the reliability of the formula. In the method of this invention, the coefficients in the formula for the concentration of the target component are obtained through nonlinear least squares fitting. When the sum of squared residuals (MSE) is minimized, the curve of the fitted function has the highest degree of overlap with the actual data distribution, and therefore the function obtained when the MSE is minimized is the most accurate.

[0094] The beneficial effects of the method of the present invention are as follows:

[0095] 1. The method of this invention establishes a multi-parameter coupled model, which can more accurately describe the battery state compared with the traditional single-parameter model, and increases the reliability of the prediction results;

[0096] 2. The battery state model describes the relationship between different physical quantities in two parts: linear and nonlinear, making the results obtained from the model more accurate and reasonable;

[0097] 3. The analysis module generates multiple sample vectors, i.e., multiple prediction results. Then, the PF algorithm processes all prediction results to make the final battery state more closely reflect the real situation.

[0098] 4. By using the risk index J(t) to reflect the risk status of the battery, the risk of the battery can be quantified and observed, enabling the method of this invention to judge the battery status based on reliable evidence. Attached Figure Description

[0099] Figure 1 Schematic diagram of the overall structure of the monitoring system;

[0100] Figure 2 A schematic diagram of the overall process of the method of this invention. Detailed Implementation

[0101] To make the objectives, technical solutions, and advantages of the present invention clearer, the technical solutions of the present invention will be clearly and completely described below through specific embodiments. Obviously, the described embodiments are only some embodiments of the present invention, 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.

[0102] Example 1:

[0103] A three-parameter coupled method for graded early warning of thermal runaway in energy storage batteries includes the following steps, and the overall monitoring process of this method is as follows: Figure 2 As shown:

[0104] Step 1. Construct a monitoring system. The monitoring system includes a sampling module, an analysis module, and a judgment module. The sampling module contains three sensors. In this embodiment, it is a warning response experiment. The monitored battery is a 48V / 100Ah lithium iron phosphate battery, operating under 1.5 times overcharge rate conditions, with a constant ambient temperature of 45℃. The sensors are a piezoelectric element for monitoring stress changes on the battery casing, an SGA-504-SKCM four-in-one gas detector installed at the battery vent, monitoring CH4 as the target component, and a K-type thermocouple for monitoring the internal temperature of the battery. All three sensors sample simultaneously, with a 1-second interval between adjacent samplings. The analysis module includes a processor, and all three sensors are signal-connected to the processor. The judgment module includes a decision unit, and the processor is also signal-connected to the decision unit. The overall structure of the monitoring system is as follows: Figure 1 As shown in the figure, 1 is the monitored battery, 2 is the mounting bracket, 3 is the piezoelectric element, 4 is the gas detector, 5 is the thermocouple, 6 is the processor, and 7 is the decision device.

[0105] Step 2. Establish a battery state model, including the following steps:

[0106] Step 2.1. Establish a model describing the changes of three physical quantities in the battery—stress, target component concentration, and temperature—over time. This includes the following steps:

[0107] Step 2.1.1. Establish a model for the change of target component concentration over time. This model is expressed by the following formula.

[0108] C g (t)+λ g ∫C g (t)dt=∫f g (t)dt+C0 (1)

[0109]

[0110] Where t is time, Cg (t) represents the concentration of the target component at time t, C0 represents the initial concentration of the target component, and λ g f is the gas diffusion coefficient. g (t) represents the gas release rate, where k1, k2, k3, β1, and β2 are constants;

[0111] The target component concentration C is obtained from formula (1). g (t) represents the rate of change at time t, which is given by the following formula.

[0112]

[0113] The constants k1, k2, k3, β1, and β2 in Step 2.1.1 depend on various factors such as battery materials, operating environment, and sensor sampling accuracy. In this embodiment, their values ​​are obtained by fitting experimental data, and the specific steps are as follows:

[0114] I. Set the target component concentration C g The rate of change over time is expressed by the finite difference method, and by combining formulas (2) and (3), we can obtain the relationship between σ and σ. p C g,p and T p formula

[0115]

[0116] II. All C g,p Substituting these values ​​into formula (21) in sequence, we obtain the values ​​on the left side of the equation (21), which are denoted as Ω. p Then, the right side of formula (21) is regarded as an algebraic expression with k1, k2, k3, β1, and β2 as unknowns, and all σ p C g,p and T p Substituting the values ​​into the equation in order, we construct the following fitting function.

[0117]

[0118] Where MSE is the sum of squared residuals, and k1', k2', k3', β1', and β2' are the variables corresponding to k1, k2, k3, β1, and β2, respectively. Changing the values ​​of these variables yields different MSE values. Each MSE value and the k1', k2', k3', β1', and β2' values ​​at which that value is obtained are recorded. The Levenberg-Marquard algorithm is used to find the k1', k2', k3', β1', and β2' values ​​corresponding to the minimum MSE value. The obtained k1', k2', k3', β1', and β2' values ​​are used as the final determined k1, k2, k3, β1, and β2 values.

[0119] Step 2.1.2. Establish a model of temperature variation over time, which is expressed by the following formula.

[0120]

[0121] Q ohmic (t)=I 2 (t)R int (T)

[0122]

[0123] Q gas (t)=ΔH g ·f g (t)

[0124] Q loss (t)=h·A·[T(t)-T env ]

[0125] Where C p Q represents the specific heat capacity of the battery casing material. ohmic (t) represents ohmic heat, Q side (t) represents the heat released by the side reaction, Q. gas (t) represents the heat released by the gas, Q loss (t) represents the heat dissipation; I(t) represents the current inside the battery at time t; R int (T) represents the battery's internal resistance, ΔH rxn For the heat of reaction, k rxn E is the reaction coefficient. a Let R be the gas internal energy, γ1 be the gas constant, and ΔH be the gas constant. g Let h be the gas friction coefficient, h be the thermal conductivity coefficient of the battery casing, A be the surface area of ​​the battery casing, and T be the total surface area of ​​the battery casing. env Ambient temperature;

[0126] Step 2.1.3. Establish a model for the stress variation of the battery casing over time. This model is expressed by the following formula.

[0127]

[0128] Where σ(t) is the stress on the battery casing, E is the elastic modulus of the battery casing material, L0 is the original length of the monitored portion of the battery casing, ΔL(t) is the elongation of the monitored portion of the battery casing, and α T Let α be the temperature-induced strain coefficient, T0 be the initial temperature of the battery, T(t) be the battery temperature, and α be the temperature-induced strain coefficient. g The strain coupling coefficient caused by gas expansion;

[0129] Based on formula (1), the rate of change of stress σ(t) at time t is obtained, i.e. It can be expressed as the following formula

[0130]

[0131] Where ε(t) is the strain rate, The value is obtained directly from the sensor;

[0132] Step 2.2. Determine the coupling relationship between the above three physical quantities in the battery's physical-chemical evolution process, specifically including the following steps:

[0133] Step 2.2.1. Construct a low-order vector X(t) and a high-order vector X(t) each containing five elements. 2 The two vectors are represented by the following formulas.

[0134]

[0135] in Let σ(t) be the first derivative of σ(t) with respect to time t, that is, the rate of change of σ(t) at time t. Let T(t) be the first derivative of T(t) with respect to time t, that is, the rate of change of T(t) at time t.

[0136] Step 2.2.2. Express the first derivative of the low-order vector X(t) with respect to time t as follows:

[0137]

[0138] Where A is the state coupling matrix, B is the input control matrix, G is the generalized coupling matrix, u(t) is the external disturbance vector, and w(t) is the system noise, which is Gaussian white noise in this embodiment. In this embodiment, u(t) can be expressed as a matrix determinant of the following form.

[0139]

[0140] Where F mech u(t) represents the mechanical load, u(t) has 3 elements, matrix B is a 5x3 matrix, and matrices A and G are represented as follows:

[0141]

[0142]

[0143] Where a ij With g ij Let i and j be elements in matrices A and G, respectively, where i and j represent the row and column numbers of the element in their respective matrices, i,j∈{1,2,3,4,5}. All elements in A and G are constants.

[0144] Step 2.2.3. From formula (8), the derivative of X(t) with respect to time t is:

[0145]

[0146] in C g The first derivative of (t) with respect to time t, i.e., C g (t) is the rate of change of time t. Let σ(t) be the second derivative with respect to time t. Let T(t) be the second derivative of T(t) with respect to time t;

[0147] Substituting the matrix determinants in formulas (8), (9), (11), (12), and (13) into formula (10), we obtain the following relation.

[0148]

[0149] Where μ1, μ2, μ3, μ4, and μ5 are all error terms, and the error terms are jointly determined by matrix B, u(t), and w(t), δ i For formula (10)G·X(t) 2 The element in the i-th row of the determinant.

[0150] Step 3. The analysis module predicts the future values ​​of the three physical quantities, including the following steps:

[0151] Step 3.1. Data measured by all sensors in the sampling module is transmitted to the processor in the analysis module in real time. The processor substitutes the received data into the battery state model to obtain predicted values ​​for three physical quantities. There are multiple predicted values ​​for each physical quantity. The process of obtaining the predicted values ​​includes the following steps:

[0152] Step 3.1.1. In formula (14), the terms on the left side of the equal sign are called lower-order terms, and the terms on the right side of the equal sign are called higher-order terms. All lower-order terms in any equation are considered to be linearly related to the higher-order terms in that equation. In order from top to bottom, select the first and second equations T(t), the third equation σ(t), and the fourth equation and the fifth equation Five lower-order terms and establish the relationships between the five lower-order terms and their corresponding terms. The linear relationships of higher-order terms are as follows: The five linear relationships are as follows:

[0153]

[0154] Where ψ1, ψ2, ψ3, ψ4, and ψ5 are all constant terms in the linear relationship;

[0155] Step 3.1.2. All derivatives in the linear relationship are expressed using the finite difference method as follows:

[0156]

[0157] Where U is the physical quantity to be calculated, U∈{σ,C} g ,T},U t U is the value of this physical quantity at a certain moment. t+1 The next time value is given, where Δt is the time interval between two consecutive records. This is the first derivative of the physical quantity at a certain moment. The first derivative at the next time step. The second derivative at a certain moment;

[0158] Step 3.1.3. Record experimental data n times during the battery thermal runaway experiment, where n is a positive integer. The experimental data includes the values ​​of all time-related variables in formulas (1) to (7), which are Q. ohmic (t), Q side (t), Q gas (t), Q loss (t), I(t), T(t), ΔL(t) The time interval between two consecutive records is equal; the recorded experimental data are sorted according to the recording time, and the data sequence number is denoted as p, where p is a positive integer and 1≤p≤n. Data from the same record have the same sequence number.

[0159] Step 3.1.4. Substitute the time-related variables from the experimental data into the corresponding formulas in formulas (1) to (7) in sequence to obtain the stress, target component concentration, and temperature at different times, as well as the first derivatives of the three physical quantities with respect to t. Then, obtain the values ​​of the five lower-order terms and the five higher-order terms at different times. Finally, fit the coefficients in the linear relationship using the least squares method to determine the linear relationship. The fitting of the coefficients in the five linear relationships is shown in the following formula.

[0160]

[0161] Step 3.1.5. Denote the current stress, target component concentration, and temperature as σ. s C g,s and T s , will σ s C g,s and T s Substituting these values ​​into the linear equation, we obtain the predicted values ​​of stress, target component concentration, and temperature at future times, denoted as σ. fc C g,f c and T fc ,

[0162]

[0163] Where σ s+1 C g,s+1 and T s+1 These are the values ​​of stress, target component concentration, and temperature after 1 second, respectively. σ s+2 and T s+2 The stress and temperature values ​​are then obtained after 2 seconds, and Δt0 is the sensor sampling time interval.

[0164] In this embodiment, the predicted values ​​of the three physical quantities are all values ​​two seconds later, i.e., σ. s+2 C g,s+2 and T s+2 The first derivatives of the three physical quantities can be obtained from formula (15), and the first derivatives of the three physical quantities can also be obtained from formulas (3), (6), and (9). Therefore, the first derivatives of the three physical quantities can all be given two values, σ. fc and T fc Therefore, there are two possible values; after knowing σ and T at a certain moment, C at that moment can be obtained from formulas (1) and (2). g Therefore, C g,fc There are four possible values.

[0165] Step 3.2. Construct sample vectors based on the predicted values; then the post-processor processes all sample vectors using the particle filter algorithm, i.e., the PF algorithm, to obtain a prediction vector. The prediction vector contains the final predicted values ​​of the three physical quantities. The process of obtaining the prediction vector includes the following steps:

[0166] Step 3.2.1. Denote the sample vector as S. q,r Where q is the index and r is the iteration number, both q and r are positive integers. Eight distinct initial sample vectors S are constructed based on the predicted values. q,0 S q,0 =[σ fc C g,fc ,T fc ] T Each initial sample vector corresponds to a unique index and 1≤q≤8; at the same time, the iteration threshold D is determined according to the number of initial sample vectors;

[0167] Step 3.2.2. Obtain the observation vector corresponding to the current sample vector using the following formula.

[0168] Z q,r =h(S q,r )+v q,r (twenty two)

[0169] Z q,r Let h be the observation vector, and v be the observation function.q,r For measuring noise;

[0170] Step 3.2.3. Calculate the weight value corresponding to each sample vector using the following formula.

[0171] W q,r =W q,r-1 ·p(Z q,r |S q,r ) (twenty three)

[0172] Among them W q,r For S q,r The corresponding weight values, where p is a probability function, p(A|B) represents the probability of event A occurring given that event B has occurred, and its value is determined by the observation model itself; initial weight values.

[0173] Step 3.2.4. Generate a new sample vector by substituting the current sample vector into the following formula.

[0174] S q,r =f(S) q,r-1 ,u' q,r )+w' q,r (twenty four)

[0175] Where u' q,r Due to external disturbances, w' q,r This refers to system process noise.

[0176] Compare r in the newly generated sample vector with D. If the condition r = D is met, proceed to Step 3.2.5. If not, return to Step 3.2.2 and use the newly generated sample vector as the current sample vector.

[0177] Step 3.2.5. Obtain the prediction vector F by weighted summation of all sample vectors. F is expressed as the following formula.

[0178]

[0179] Step 4. Based on the final predicted values ​​of the three physical quantities, the processor derives the current battery risk index J(t), as follows:

[0180] Step 4.1. Construct the risk function for calculating the risk index J(t). The risk function is expressed as follows:

[0181] J(t) = w σ ·σ'(t)+w g ·C g '(t)+w T ·T'(t) (26)

[0182] Where σ'(t) and C g T'(t) and T'(t) are the predicted values ​​of stress, target component concentration, and temperature at time t, respectively. The predicted value is the ratio of the final predicted value of a physical quantity to its experimentally recorded maximum value. σ w g w T These are the weight values ​​corresponding to the three physical quantities: stress, target component concentration, and temperature.

[0183] Step 4.2. Determine w based on actual working conditions σ w g w T The specific values ​​of the three components, and the three weight values ​​satisfy w σ +w g +w T In this embodiment, the condition is w = 1. σ =0.3, w g =0.4, w T =0.3;

[0184] Step 4.3. Based on the final predicted value at the current time, calculate the difference between the current time and the predicted value σ'(t), C. g Substitute the three predicted values, T'(t) and T'(t), into formula (20) to obtain the current time J(t).

[0185] Step 5. The processor transmits the value of J(t) to the decision-making module in real time. The decision-making module makes a judgment on the current state of the battery based on the value of J(t) and displays the judgment result on the screen. The specific judgment process is as follows:

[0186] Step 5.1. Based on the battery operating conditions, take two critical values ​​in the range [0,1], namely 0.4 and 0.7. Divide the range [0,1] into three intervals: [0,0.4], [0.4,0.7], and [0.7,1]. The interval [0,0.4] corresponds to a low-risk warning state, the interval [0.4,0.7] corresponds to a medium-risk warning state, and the interval [0.7,1] corresponds to a high-risk warning state.

[0187] Step 5.2. Compare J(t) with the magnitude of all critical values ​​to determine the interval in which J(t) is located, and then determine the warning state at the current moment. In this embodiment, after the experiment starts, the internal temperature of the battery only fluctuates within a small range, the CH4 concentration gradually increases, and the stress on the outer shell increases slightly. At the 28th minute, the risk index J(t) = 0.55, and the judgment module outputs a medium risk warning state. Before this, the judgment module has been outputting a low risk warning state.

Claims

1. A three-parameter coupled method for graded early warning of thermal runaway in energy storage batteries, characterized in that, Specifically, the following steps are included: Step 1. Build a monitoring system. The monitoring system includes a sampling module, an analysis module, and a judgment module. The sampling module and the analysis module are connected by a signal, and the analysis module and the judgment module are also connected by a signal. The sampling module includes three types of sensors: a stress sensor, a gas concentration sensor, and a temperature sensor. All sensors sample simultaneously and the time interval between two adjacent samplings is equal. Step 2. Establish a battery state model, including the following steps: Step 2.

1. Establish a model for the time-varying nature of three physical quantities related to the battery: stress, target component concentration, and temperature; Step 2.

2. Determine the coupling relationship between the above three physical quantities in the physical-chemical evolution process of the battery; Step 3. The analysis module predicts the future values ​​of the three physical quantities, including the following steps: Step 3.

1. The data measured by all sensors in the sampling module are transmitted to the analysis module in real time. The analysis module substitutes the received data into the battery state model to obtain the predicted values ​​of three physical quantities. There are multiple predicted values ​​for each physical quantity. Step 3.

2. Construct multiple sample vectors based on the predicted values. Then, the analysis module processes all sample vectors using the particle filtering algorithm, i.e., the PF algorithm, to obtain a prediction vector. The prediction vector contains the final predicted values ​​of the three physical quantities. Step 4. Based on the final predicted values ​​of the three physical quantities, the analysis module derives the current risk index J(t) of the battery. Step 5. The analysis module transmits the J(t) value to the judgment module in real time. The judgment module makes a judgment on the current state of the battery based on the J(t) value and outputs the judgment result.

2. The three-parameter coupled method for graded early warning of thermal runaway in energy storage batteries as described in claim 1, characterized in that, Step 2.1 specifically includes the following steps: Step 2.1.

1. Establish a model of temperature change over time, which is expressed by the following formula. Q ohmic (t)=I 2 (t)R int (T) Q gas (t)=ΔH g ·f g (t) Q loss (t)=h·A·[T(t)-T env ] Where C p Q represents the specific heat capacity of the battery casing material. ohmic (t) represents ohmic heat, Q side (t) represents the heat released by the side reaction, Q. gas (t) represents the heat released by the gas, Q loss (t) represents the heat dissipation; I(t) represents the current inside the battery at time t; R int (T) represents the battery's internal resistance, ΔH rxn For the heat of reaction, k rxn E is the reaction coefficient. a Let R be the gas internal energy, γ1 be the gas constant, and ΔH be the gas constant. g Let be the gas friction coefficient, h be the thermal conductivity coefficient of the battery casing, A be the surface area of ​​the battery casing, and T(t) be the battery temperature at time t. env Ambient temperature; Step 2.1.

2. Establish a model for the change of target component concentration over time. This model is expressed by the following formula. C g (t)+λ g ∫C g (t)dt=∫f g (t)dt+C0 (3) Where t is time, C g (t) represents the concentration of the target component at time t, C0 represents the initial concentration of the target component, and λ g f is the gas diffusion coefficient. g (t) represents the gas release rate, where k1, k2, k3, β1, and β2 are constants; The target component concentration C is obtained from formula (3). g (t) represents the rate of change at time t, which is given by the following formula. Step 2.1.

3. Establish a model for the stress variation of the battery casing over time. This model is expressed by the following formula. Where σ(t) is the stress on the battery casing, E is the elastic modulus of the battery casing material, L0 is the original length of the monitored portion of the battery casing, ΔL(t) is the elongation of the monitored portion of the battery casing, and α T Let α be the temperature-induced strain coefficient, T0 be the initial temperature of the battery, T(t) be the battery temperature, and α be the temperature-induced strain coefficient. g The strain coupling coefficient caused by gas expansion; Based on formula (1), the rate of change of stress σ(t) at time t is obtained, i.e. It can be expressed as the following formula Where ε(t) is the strain rate, The value is obtained directly from the sensor.

3. The three-parameter coupled method for graded early warning of thermal runaway in energy storage batteries as described in claim 2, characterized in that, Step 2.2 specifically includes the following steps: Step 2.2.

1. Construct a low-order vector X(t) and a high-order vector X(t) each containing five elements. 2 The two vectors are represented by the following formulas. in Let σ(t) be the first derivative of σ(t) with respect to time t, that is, the rate of change of σ(t) at time t. Let T(t) be the first derivative of T(t) with respect to time t, that is, the rate of change of T(t) at time t. Step 2.2.

2. Express the first derivative of the low-order vector X(t) with respect to time t as follows: Where A is the state coupling matrix, B is the input control matrix, G is the generalized coupling matrix, u(t) is the external disturbance vector, and w(t) is the system noise; let b be the number of elements in u(t), where b is a positive integer, then matrix B is a 5xb matrix. Matrices A and G are expressed in the following forms. Where a ij With g ij Let i and j be elements in matrices A and G, respectively, where i and j represent the row and column numbers of the element in their respective matrices, i,j∈{1,2,3,4,5}. All elements in A and G are constants. Step 2.2.

3. From formula (8), the derivative of X(t) with respect to time t is: in C g The first derivative of (t) with respect to time t, i.e., C g (t) is the rate of change of time t. Let σ(t) be the second derivative with respect to time t. Let T(t) be the second derivative of T(t) with respect to time t; Substituting the matrix determinants in formulas (8), (9), (11), (12), and (13) into formula (10), we obtain the following relation. Where μ1, μ2, μ3, μ4, and μ5 are all error terms, and the error terms are jointly determined by matrix B, u(t), and w(t), δ i For formula (10)G·X(t) 2 The element in the i-th row of the determinant.

4. The three-parameter coupled method for graded early warning of thermal runaway in energy storage batteries as described in claim 3, characterized in that, The process of obtaining the predicted value in Step 3.1 specifically includes the following steps: Step 3.1.

1. In formula (14), the terms on the left side of the equal sign are called lower-order terms, and the terms on the right side of the equal sign are called higher-order terms. All lower-order terms in any equation are considered to be linearly related to the higher-order terms in the equation. Select multiple lower-order terms and establish a linear relationship between each lower-order term and its corresponding higher-order term. Step 3.1.

2. The derivative in the linear relationship can be expressed using the finite difference method as follows: Where U is the physical quantity to be calculated, U∈{σ,C} g ,T},U t U is the value of this physical quantity at a certain moment. t+1 This represents the value at the next moment, where Δt is the time interval between two consecutive records. This is the first derivative of the physical quantity at a certain moment. The first derivative at the next time step. The second derivative at a certain moment; Step 3.1.

3. Record experimental data n times during the battery thermal runaway experiment, where n is a positive integer. The experimental data includes the values ​​of all time-related variables in formulas (1) to (7). The time interval between two adjacent records is equal. Sort the experimental data according to the order of recording time, and record the data sequence number as p, where p is a positive integer and 1≤p≤n. All values ​​recorded in the same time have the same sequence number. Step 3.1.

4. Substitute the time-related variables in the experimental data into the corresponding formulas in formulas (1) to (7) in order to obtain the stress, target component concentration and temperature at different times, as well as the first derivatives of the three physical quantities with respect to t, and then obtain the values ​​of the low-order and high-order terms in the linear relationship at different times; then fit the coefficients in the linear relationship by the least squares method to determine the linear relationship. Step 3.1.

5. The stress, target component concentration, and temperature at the current moment are denoted as σ. s C g,s and T s , will σ s C g,s and T s Substituting these values ​​into the linear equation, we obtain the predicted values ​​of stress, target component concentration, and temperature at future times, denoted as σ. fc C g,fc and T fc All three physical quantities have multiple predicted values.

5. The three-parameter coupled method for graded early warning of thermal runaway in energy storage batteries as described in claim 4, characterized in that, The process of obtaining the prediction vector in Step 3.2 specifically includes the following steps: Step 3.2.

1. Denote the sample vector as S. q,r Where q is the index and r is the iteration number, both q and r are positive integers. Based on the predicted values, m initial sample vectors S are constructed. q,0 S q,0 =[σ fc C g,fc ,T fc ] T Each initial sample vector corresponds to a unique index and 1≤q≤m. Any two S... q,0 All are unequal; at the same time, the iteration threshold D is determined based on the number of initial sample vectors; Step 3.2.

2. Obtain the observation vector corresponding to the current sample vector using the following formula. From q,r =h(S q,r )+v q,r (16) Z q,r Let h be the observation vector, and v be the observation function. q,r For measuring noise; Step 3.2.

3. Calculate the weight value corresponding to each sample vector using the following formula. W q,r =W q,r-1 ·p(Z q,r |S q,r ) (17) Among them W q,r For S q,r The corresponding weight values, where p is a probability function, p(A|B) represents the probability of event A occurring given that event B has occurred, and its value is determined by the observation model itself; initial weight values. Step 3.2.

4. Generate a new sample vector by substituting the current sample vector into the following formula. S q,r =f(S q,r-1 ,u’ q,r )+w’ q,r (18) Where u' q,r Due to external disturbances, w' q,r This refers to system process noise. Compare r in the newly generated sample vector with D. If the condition r = D is met, proceed to Step 3.2.

5. If not, return to Step 3.2.2 and use the newly generated sample vector as the current sample vector. Step 3.2.

5. Obtain the prediction vector F by weighted summation of all sample vectors. F is expressed as the following formula.

6. The three-parameter coupled method for graded early warning of thermal runaway in energy storage batteries as described in claim 5, characterized in that, The specific process for obtaining the battery risk index J(t) at the current moment in Step 4 is as follows: Step 4.

1. Construct the risk function for calculating the risk index J(t). The risk function is expressed as follows: J(t)=w σ ·σ’(t)+w g ·C g ’(t)+w T ·T’(t) (20) Where σ'(t) and C g T'(t) and T'(t) are the predicted values ​​of stress, target component concentration, and temperature at time t, respectively. The predicted value is the ratio of the final predicted value of a physical quantity to its experimentally recorded maximum value. σ w g w T These are the weight values ​​corresponding to the three physical quantities: stress, target component concentration, and temperature. Step 4.

2. Determine w based on actual working conditions σ w g w T The specific values ​​of the three components, and the three weight values ​​satisfy w σ +w g +w T =1 condition; Step 4.

3. Based on the final predicted value at the current time, calculate the difference between the current time and the predicted value σ'(t), C. g Substitute the three predicted values, T'(t) and T'(t), into formula (20) to obtain the current time J(t).

7. The three-parameter coupled method for graded early warning of thermal runaway in energy storage batteries as described in claim 6, characterized in that, The process for determining the current state of the battery in Step 5 is as follows: Step 5.

1. Based on the battery operating conditions, take multiple critical values ​​within the range of [0,1], and divide the range of [0,1] into multiple intervals through these critical values. Each interval corresponds to a warning state. Step 5.

2. Compare J(t) with the magnitude of all critical values ​​to determine the interval in which J(t) is located, and then determine the warning status at the current moment.

8. The three-parameter coupled method for graded early warning of thermal runaway in energy storage batteries as described in claim 7, characterized in that, The values ​​of constants k1, k2, k3, β1, and β2 in Step 2.1.1 are determined through the following steps: I. Set the target component concentration C g The rate of change over time is expressed by the finite difference method, and by combining formulas (4) and (5), we can obtain the relationship between σ and σ. p C g,p and T p formula II. All C g,p Substituting these values ​​into formula (21) in sequence, we obtain the values ​​on the left side of the equation (21), which are denoted as Ω. p Then, the right side of formula (21) is regarded as an algebraic expression with k1, k2, k3, β1, and β2 as unknowns, and all σ p C g,p and T p Substituting the values ​​into the equation in order, we construct the following fitting function. Where MSE is the sum of squared residuals, and k1', k2', k3', β1', and β2' are the variables corresponding to k1, k2, k3, β1, and β2, respectively. Changing the values ​​of these variables yields different MSE values. Each MSE value and the k1', k2', k3', β1', and β2' at which that value is obtained are recorded. Then, the minimum MSE value is found, and the k1', k2', k3', β1', and β2' at which the minimum MSE value is obtained are taken as the final determined values ​​of k1, k2, k3, β1, and β2.

Images