Next Article in Journal
Optimization of Coal Production Based on the Modeling of the Jig Operation
Previous Article in Journal
Class Thresholds Pre-Definition by Clustering Techniques for Applications of ELECTRE TRI Method
Previous Article in Special Issue
Atmospheric Ecology Modeling for the Sustainable Development of the Urban Environment
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Generalized Method of Mathematical Prototyping of Energy Processes for Digital Twins Development

1
Department of Electrical Engineering and Aviation Electrical Equipment, Moscow State Technical University of Civil Aviation (MSTU CA), 125993 Moscow, Russia
2
Department of Chemical Engineering, Ariel University, Ariel 4070000, Israel
*
Author to whom correspondence should be addressed.
Energies 2023, 16(4), 1933; https://doi.org/10.3390/en16041933
Submission received: 13 December 2022 / Revised: 20 January 2023 / Accepted: 11 February 2023 / Published: 15 February 2023
(This article belongs to the Special Issue Smart Energy and Sustainable Environment)

Abstract

:
The use of digital twins in smart power systems at the stages of the life cycle is promising. The dynamics of such systems (smart energy renewable sources, smart energy hydrogen systems, etc.), are determined mainly by the physical and chemical processes occurring inside the systems. The basis for developing digital twins is reliable mathematical models of the systems. In the present paper, the authors present a method of energy processes mathematical prototyping—an overall approach to modeling processes of various physical and chemical natures based on modern non-equilibrium thermodynamics, mechanics, and electrodynamics. Controlled parameters are connected with measured ones by developing a theoretically correct system of process dynamics equations with accuracy up to the experimentally studied properties of substances and processes. Subsequent transformation into particular mathematical models of a specific class of systems makes this approach widely applicable. The properties of substances and processes are given in the form of functional dependencies on the state of the system up to experimentally determined constant coefficients. The authors consider algorithms for identifying the constant coefficients of the functions of substances and processes properties, which complement the proposed unified approach of designing models of various physical and chemical nature systems.

1. Introduction

Utilization of renewable energy sources is primarily associated with the requirements of controlling a distributed system with unstable parameters and implies the need to automate the management of such systems. The complexity of the system’s management appears from the flow of nonlinear physical and chemical processes in the system, which significantly depend not only on operating modes but also on external conditions (such as ambient temperature, wind, and solar insolation). Nowadays, promising solutions are smart control systems based on the creation and implementation of digital twins, in which the parameters of the objects are monitored, and their mathematical models are reconstructed. The basis for digital twin development is reliable mathematical models of the systems.
For the functioning of digital twins, it is necessary to use mathematical models in the form of dependencies of the controlled parameters on the measured ones. The measured parameters are input to the model, and the controlled parameters are returned at the output, according to which practical decisions are made [1,2,3,4,5,6,7,8,9]. In particular, probabilistic characteristics can be obtained at the output, for example, in problems related to the reliability and safety of technical systems [2,3,9].
Methods of identification theory [10,11], methods of machine learning [12,13,14,15,16,17,18,19,20], including deep learning based on neural networks [14], and symbolic regression [15,16,17,18,19,20,21] are often used to obtain mathematical models of renewable energy systems (RES). However, all these modeling methods belong to the class of simulation models that do not consider actual physical and chemical processes and therefore do not guarantee their correctness in the entire range of operating conditions. To obtain correct (not contradicting the general physical laws) models of systems, it is necessary to set some restrictions arising from the general physical laws and the system functioning mechanism [10,11,12,13,14,15,16,17,18,19,20,21,22,23].
Correct mathematical models are formed based on physical laws and the study of the characteristics of the object of study: the structure of the object and the physical and chemical processes occurring in it. That is why the approach based on modeling the dynamics of physical and chemical processes in the system has found wide practical application, in contrast to the simulation approach, which considers the system as a “black box” [4,8,9,10,22,23]. Another aspect of resorting to this mathematical model approach is interpretability and explainability for subsequent decision-making.
In the general case, the physical laws, based on which the system model is built, are conservation laws or the connection between the causes of processes (internal disturbances, internal forces) and the rates of these processes (flows) [10,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39]. This makes it possible to set a class of models for any RES in the form of a connection between internal forces and flows, as well as a connection between flows and the rates of change in its state coordinates based on conservation laws [10,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39]. These models are generally built up to the experimentally determined properties of substances and processes, which generally depend on the state of the system [24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39]. Thus, the model of the particular RES instance should be a system of equations: differential equations for the dynamics of physical and chemical processes, as well as equations for measured and controlled parameters.
The features of building a system model described above make it possible to develop a unified approach to build models of systems of various physical and chemical nature [39,40], based on mechanics [29,30], including continuum mechanics [28,31], electrodynamics [28,32], the theory of electric and magnetic circuits [32], on modern non-equilibrium thermodynamics [33,34,35,36,37,38] and incorporating methods of identification theory, machine learning methods, including deep learning based on neural networks, and on symbolic regression [15,16,17,18,19,20,21]. Such models take measured parameters as input and return controlled parameters that have practical value as output [1,2,3,4,5,6,7,8,9]. This work is devoted to the development of this approach.

2. Description of the Energy Processes Mathematical Prototyping Method

As noted above, the general fundamental laws are conservation laws and the link between system internal disturbances and the speed of physical and chemical processes inside the system [10,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39]. The main conservation law in mechanical systems is the momentum conservation law [28,30,31]; in the theory of electric and magnetic circuits, as well as in electrodynamics—the electric charge and magnetic flux [28,32] conservation law; in modern non-equilibrium thermodynamics—the energy conservation law (the first law of thermodynamics), the mass and stoichiometric ratios conservation law [28,31,33,34,35,36,37,38]. The energy conservation law is a general physical conservation law [28,30,31,32,33,34,35,36,37,38].
From the point of view of modern non-equilibrium thermodynamics, the cause and necessary condition for the occurrence of physicochemical processes in an arbitrary system are thermodynamic forces [33,35,36,37], which are internal disturbances [10,24,25,26,27,28,29,33,35,36,37]. Examples of thermodynamic forces are [33,35,36,37]: normalized temperature difference (reciprocal temperature difference), chemical potential difference, and chemical affinity. In mechanics, internal disturbances are the difference in velocities, which causes the friction force (momentum flux due to friction), as well as potential forces [10,28,29,30,31]. In electrodynamics, the theory of electrical and magnetic circuits, such internal disturbances are the difference in electrical potentials that carry an electric charge, and currents through inductive elements that create magnetic fluxes [10,28,29,32].
Flows in physical and chemical systems are the rates of physical and chemical processes [33,35,36,37]: in mechanics—mechanical forces (due to Newton's second law—momentum flows) and speeds [10,28,29,30,31]; in electrodynamics, theories of electric and magnetic circuits—electric currents that carry an electric charge and EMF of electromagnetic induction [10,28,29,32]. Examples of flows in modern nonequilibrium thermodynamics are [33,35,36,37]: heat flow, substance flow (diffusion flow, phase transition rate), and chemical reaction rate. These internal perturbations do not unambiguously determine the flows, the flows are also determined by the properties of the systems that do not depend on the perturbations. [10,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39].
In the general case of processes of different physical and chemical nature, the authors proposed a unified approach to their description within the framework of mechanics [29,30], including continuum mechanics [28,31], electrodynamics [28,32], the theory of electrical and magnetic circuits [32], modern nonequilibrium thermodynamics [33,34,35,36,37,38]—it is a method of mathematical prototyping that does not contradict the conservation laws and the second law of thermodynamics [39]. The equations of this method can be written as [39]:
d x t d t = B x t , U t δ Δ x t d t + d x t d t ,   δ Δ x t d t = A x t , U t Δ F x t , U t ,
Δ F ( x , U ) = B T ( x , U ) F ( x , U ) ,   F ( x , U ) = x W ( x , U ) ,
y t = g y x t , U t ,   z t = g z x t , U t ,
where x —the coordinates of the system state; Δ x —coordinates of processes in the system; δ Δ x t / d t —speed of physical and chemical processes in the system; B x , U —system topology matrix; d x t / d t is the component of the change rate of the system state parameters, due to its interaction with external systems (it also additively includes a random component of these external flows [36,37,38]); Δ F x , U —dynamic forces (internal disturbances), which are the cause and necessary condition for the flow of physical and chemical processes in the system; F x , U —partial derivatives of the free energy W x , U by the coordinates of the state x , taken with the sign “-”; A x , U —positively defined (in particular cases, if there is inertia in the system, non-degenerated, non-negative defined) dissipative matrix; U —system parameters that do not change as a result of the processes in the system, but change as a result of external influences; y —measured parameters of the system; z —controlled parameters of the system. The measured and controlled parameters of the system (y and z) can be both a function of the system state and functionals of the system dynamics (dynamics x t and U t ). To implement the system of Equations (1)–(3) in numerical form, it is necessary to have [39]:
  • Topology matrix B x , U ;
  • The expression for the free energy W x , U , expressed in terms of x , or its partial derivatives by the state coordinates x , taken with the “-” sign F x , U , satisfying the total differential condition;
  • Positively defined dissipative matrix A x , U ;
  • Functions g y x , U and g z x , U , obtained from the definition of the measured y and controlled z parameters of the system, respectively.
If we denote a part of the parameters y and z as y ¯ and z ¯ , which are functions of the state x :
y ¯ = g ¯ y x , U ,   z ¯ = g ¯ z x , U ,
then if the right-hand sides of (1) and (2) directly depend on y ¯ and z ¯ , we represent the quantities included in (1) and (2) as compositions of functions [38]:
B x , U B y ¯ , z ¯ , U B g ¯ y x , U , g ¯ z x , U , U
W x , U W y ¯ , z ¯ , U W g ¯ y x , U , g ¯ z x , U , U
F x , U F y ¯ , z ¯ , U F g ¯ y x , U , g ¯ z x , U , U
A x , U A y ¯ , z ¯ , U A g ¯ y x , U , g ¯ z x , U , U
y ˜ t = g ˜ y x t , U t g ˜ y y ¯ t , z ¯ t , U t g ˜ y g ¯ y x t , U t , g ¯ z x t , U t , U t
z ˜ t = g ˜ z x t , U t g ˜ z y ¯ t , z ¯ t , U t g ˜ z g ¯ y x t , U t , g ¯ z x t , U t , U t
where y ˜ and z ˜ are another part of the parameters y and z , respectively:
y = y ¯ T y ˜ T T ,   z = z ¯ T z ˜ T T .
In the vast majority of cases, the absolute values of the coordinates of the state x are not determined (for example, the internal energy [38,41]), so, as can be seen from (4), it is advisable to set the quantities y ¯ and z ¯ in differential form [38]:
d y ¯ = C y , x y ¯ , z ¯ , U d x + C y , U y ¯ , z ¯ , U d U ,   d z ¯ = C z , x y ¯ , z ¯ , U d x + C z , U y ¯ , z ¯ , U d U ,
where C y , x y ¯ , z ¯ , U , C y , U y ¯ , z ¯ , U , C z , x y ¯ , z ¯ , U , C z , U y ¯ , z ¯ , U are the Jacobians of the right-hand sides of (4) with respect to x and U . Hence, based on (1) and (12) we obtain the following expressions [38]:
d y ¯ t d t = B ˜ y y ¯ t , z ¯ t , U t δ Δ x t d t + d y ¯ t d t ,   d z ¯ t d t = B ˜ z y ¯ t , z ¯ t , U t δ Δ x t d t + d z ¯ t d t ,
where [38]:
d y ¯ t d t = C y , x y ¯ t , z ¯ t , U t d x t d t + C y , U y ¯ t , z ¯ t , U t d U t d t ,
d z ¯ t d t = C z , x y ¯ t , z ¯ t , U t d x t d t + C z , U y ¯ t , z ¯ t , U t d U t d t ,
B ˜ y y ¯ , z ¯ , U = C y , x y ¯ , z ¯ , U B y ¯ , z ¯ , U ,   B ˜ z y ¯ , z ¯ , U = C z , x y ¯ , z ¯ , U B y ¯ , z ¯ , U .
By virtue of (5)–(10), system (1)–(3) takes the form [38]:
δ Δ x t d t = A y ¯ t , z ¯ t , U t Δ F y ¯ t , z ¯ t , U t ,   Δ F y ¯ , z ¯ , U = B T y ¯ , z ¯ , U F y ¯ , z ¯ , U ,
y ˜ t = g ˜ y y ¯ t , z ¯ t , U t ,   z ˜ t = g ˜ z y ¯ t , z ¯ t , U t .
The system of Equations (13)–(18) is more convenient for practical applications [38].
Thus, in order to obtain the transformed equations of the mathematical prototyping method (13)–(18), it is necessary to obtain functional dependences of the system state from the experimental data for substances and processes properties and system topology [38,39]:
  • Topology matrix B y ¯ , z ¯ , U ;
  • Positively defined dissipative matrix A y ¯ , z ¯ , U ;
  • Jacobi matrices C y , x y ¯ , z ¯ , U , C y , U y ¯ , z ¯ , U and C z , x y ¯ , z ¯ , U , C z , U y ¯ , z ¯ , U of measured y ¯ and controlled z ¯ parameters of the system respectively, satisfying the conditions of the total differential [38]:
    k = 1 n y ¯ C y , x , i y ¯ , z ¯ , U y ¯ k C y , x , k , j y ¯ , z ¯ , U + k = 1 n z ¯ C z , x , i y ¯ , z ¯ , U z ¯ k C z , x , k , j y ¯ , z ¯ , U k = 1 n y ¯ C y , x , j y ¯ , z ¯ , U y ¯ k C y , x , k , i y ¯ , z ¯ , U + k = 1 n z ¯ C z , x , j y ¯ , z ¯ , U z ¯ k C z , x , k , i y ¯ , z ¯ , U ,        j = 1 , i 1 ,    i = 2 , n x ,    n x = dim x ,    n y ¯ = dim y ¯ ,    n z ¯ = dim z ¯ ,    n x = n y ¯ + n z ¯ ,
    k = 1 n y ¯ C y , x , i y ¯ , z ¯ , U y ¯ k C y , U , k , j y ¯ , z ¯ , U + k = 1 n z ¯ C z , x , i y ¯ , z ¯ , U z ¯ k C z , U , k , j y ¯ , z ¯ , U + C y , x , i y ¯ , z ¯ , U U j       k = 1 n y ¯ C y , U , j y ¯ , z ¯ , U y ¯ k C y , x , k , i y ¯ , z ¯ , U + k = 1 n z ¯ C z , U , j y ¯ , z ¯ , U z ¯ k C z , x , k , i y ¯ , z ¯ , U ,    j = 1 , n U ,    i = 1 , n x ,        n U = dim U ,
            k = 1 n y ¯ C y , U , i y ¯ , z ¯ , U y ¯ k C y , U , k , j y ¯ , z ¯ , U + k = 1 n z ¯ C z , U , i y ¯ , z ¯ , U z ¯ k C z , U , k , j y ¯ , z ¯ , U + C y , U , i y ¯ , z ¯ , U U j k = 1 n y ¯ C y , U , j y ¯ , z ¯ , U y ¯ k C y , U , k , i y ¯ , z ¯ , U + k = 1 n z ¯ C z , U , j y ¯ , z ¯ , U z ¯ k C z , U , k , i y ¯ , z ¯ , U + C z , U , j y ¯ , z ¯ , U U i , j = 1 , i 1 ,    i = 2 , n U ,
    where:
    C y , x y ¯ , z ¯ , U = C y , x , 1 y ¯ , z ¯ , U C y , x , n x y ¯ , z ¯ , U ,
    C z , x y ¯ , z ¯ , U = C z , x , 1 y ¯ , z ¯ , U C z , x , n x y ¯ , z ¯ , U ,
    C y , U y ¯ , z ¯ , U = C y , U , 1 y ¯ , z ¯ , U C y , U , n x y ¯ , z ¯ , U ,
    C z , U y ¯ , z ¯ , U = ( C z , U , 1 y ¯ , z ¯ , U C z , U , n x y ¯ , z ¯ , U ) ,
    C y , x y ¯ , z ¯ , U = C y , x , 1 , 1 y ¯ , z ¯ , U C y , x , 1 , n x y ¯ , z ¯ , U C y , x , n y ¯ , 1 y ¯ , z ¯ , U C y , x , n y ¯ , n x y ¯ , z ¯ , U ,
    C z , x y ¯ , z ¯ , U = C z , x , 1 , 1 y ¯ , z ¯ , U C z , x , 1 , n x y ¯ , z ¯ , U C z , x , n y ¯ , 1 y ¯ , z ¯ , U C z , x , n y ¯ , n x y ¯ , z ¯ , U ,
    C y , U y ¯ , z ¯ , U = C y , U , 1 , 1 y ¯ , z ¯ , U C y , U , 1 , n x y ¯ , z ¯ , U C y , U , n y ¯ , 1 y ¯ , z ¯ , U C y , U , n y ¯ , n x y ¯ , z ¯ , U ,
    C z , U y ¯ , z ¯ , U = C z , U , 1 , 1 y ¯ , z ¯ , U C z , U , 1 , n x y ¯ , z ¯ , U C z , U , n y ¯ , 1 y ¯ , z ¯ , U C z , U , n y ¯ , n x y ¯ , z ¯ , U ;
  • Partial derivatives F y ¯ , z ¯ , U by the coordinates of the state x of free energy W x , U , satisfying the conditions of the total differential [38]:
    k = 1 n y ¯ F i y ¯ , z ¯ , U y ¯ k C y , x , k , j y ¯ , z ¯ , U + k = 1 n z ¯ F i y ¯ , z ¯ , U z ¯ k C z , x , k , j y ¯ , z ¯ , U k = 1 n y ¯ F j y ¯ , z ¯ , U y ¯ k C y , x , k , i y ¯ , z ¯ , U + k = 1 n z ¯ F j y ¯ , z ¯ , U z ¯ k C z , x , k , i y ¯ , z ¯ , U , j = 1 , i 1 ,    i = 2 , n x ;
  • Functions of measured g ˜ y y ¯ , z ¯ , U and controlled g ˜ z y ¯ , z ¯ , U parameters (which can be functionals of y ¯ t , z ¯ t , and U t dynamics).
Further, for the dissipative matrix, partial derivatives of the free energy, the Jacobi matrices of the measured and controlled parameters specified in differential form, the topology matrix, and functions for the remaining measured and controlled parameters, it is necessary to specify a class of analytical expressions up to constant coefficients. In this class, for any desired function, there must necessarily be an analytic expression that approximates the desired function with any given accuracy. Such classes, for example, are:
  • Power polynomials, whose convergence is guaranteed by the Weierstrass theorem on the uniform approximation of functions by polynomials [21];
  • Classes of inductively generating functions [15] obtained by symbolic regression methods; convergence, in this case, is confirmed by the universal approximation theorem (Cybenko’s theorem) [42,43];
  • Classes of interpolation expressions (linear, cubic splines, Lagrange interpolation polynomials, etc.) [44].
The values of constant coefficients, in the general case, are determined in such a way that the measured parameters y determined from (1)–(3) or from (13)–(18) coincide with the corresponding values of these parameters y E , obtained by direct measurement [10,11,37]:
y t j = y E t j ,   j = 1 , N t ,
where t j , j = 1 , N t —discrete moments of time; N t —number of discrete moments of time. The coefficients of approximating analytical expressions for the properties of substances and processes that satisfy the relevant restrictions can be determined from experimental data using (13)–(18), (23), by minimizing the objective function Q y t [10,11,37]:
Q y t = 1 2 i = 1 n j = 1 N t , i y i t j y i E t j T L ˜ i y i t j y i E t j ,
where L ˜ i , i = 1 , n is a positively defined symmetric matrix; N t , i , i = 1 , n —the number of discrete moments of time in each i -th operation mode; n is the number of modes in which experimental data are taken y i E t j , j = 1 , N t , i , i = 1 , n [10,11,37]. Minimization of the objective function Q , defined by virtue of (13)–(18), (23), (24) can be carried out in different ways [45], for example, in parts [45]. Hence, the analytical expressions for the state function of the dissipative matrix should be given in the form [37]:
A y ¯ , z ¯ , U = A h a y ¯ , z ¯ , U , p a + j = 0 m a A j i = 1 n x h a , i y ¯ , z ¯ , U , p a n i , j n i , j ! , h a y ¯ , z ¯ , U , p = h a , 1 y ¯ , z ¯ , U , p a h a , n x y ¯ , z ¯ , U , p a T ,
where A h a is the positively defined (non-negative-defined in the case of inertia in the system) basic dissipative matrix; h a , i = h a , i y ¯ , z ¯ , U , p a > 0 , j = 0 , m a are some variables determined by the state of the system; m a is the number of basic additive functions; A j , j = 0 , m a are constant non-negative defined matrices; p a are the parameters by which, along with A j , j = 0 , m a the objective function is optimized Q . For independent components of other properties of substances and processes (partial derivatives of free energy F y ¯ , z ¯ , U , Jacobi matrices C y , x y ¯ , z ¯ , U , C y , U y ¯ , z ¯ , U and C z , x y ¯ , z ¯ , U , C z , U y ¯ , z ¯ , U observed y ¯ and controlled z ¯ parameters, respectively, which are taken in differential form, topology matrices B y ¯ , z ¯ , U ), analytical expressions are given in the form [37]:
H y ¯ , z ¯ , U = H h y ¯ , z ¯ , U , p + j = 0 m c j i = 1 n x h i y ¯ , z ¯ , U , p n i , j n i , j ! , h y ¯ , z ¯ , U , p = h 1 y ¯ , z ¯ , U , p h n x y ¯ , z ¯ , U , p T ,
where H h is the basic component; c j are constant coefficients, j = 0 , m ; h i = h i y ¯ , z ¯ , U , p , i = 1 , m -some variables determined by the state of the system; p are the parameters by which, along with c j the objective function is optimized Q .
A positively defined dissipative matrix A y ¯ , z ¯ , U , in the case of physical and chemical processes, is a kinetic matrix constructed through its reversible and positive irreversible components [36,37], which are given in the form (26) [37]. The positivity of the irreversible components guarantees the positive definiteness of the kinetic matrix [36,37]. In the case of mechanics, system (1)–(3) transforms into the Hamiltonian equations [30,39], in this case, the dissipative matrix is constructed based on the friction coefficients and transfer functions [39]. The positivity of the friction coefficients, also given in the form (26), guarantees the non-negative definiteness of the dissipative matrix of mechanical systems (included in the Hamilton equations) [30,39]. In the case of electrodynamics, in terms of the transfer of electric charge through an anisotropic crystal, in terms of the Hall effect, it is also expedient to build a dissipative matrix through reversible and irreversible components [35,36]. In the case of electric and magnetic circuits, the dissipative matrix is built based on the resistances in the electric circuit and the way they are connected [39]. For physical and chemical processes, the existing model of the specific nature process can be converted to the form (1)–(3) or to the form (13)–(18), in this case the analytical expression of the kinetic matrix is constructed in accordance with (25) [36].
Thus, building a model of physical and chemical processes in the system is reduced to optimizing the objective function Q y τ given in the form (13)–(18), (23)–(26), according to the constant parameters included in (25) and (26). Such a procedure, as it is easy to see from (13)–(18), (23)–(26), is reduced to solving a system of ordinary differential equations. Currently, in order to solve a system of ordinary differential equations, the following methods [46] can be used:
  • Step methods based on calculating the values of a dynamic variable at subsequent time points from the previous values of this variable;
  • Special methods based on the approximation of solutions to a system of differential equations by analytical expressions.
The main advantage of step methods is their universality for any system of ordinary differential equations because as the integration step tends to zero, the approximate solution uniformly converges to the exact one [46]. However, these methods have a disadvantage—the complexity of calculations [46]. Special methods based on specifying the solution by an analytic expression followed by the search for the constant parameters of this expression are free from this shortcoming [46]. This approach is less labor intensive than step methods. However, to apply this approach, it is necessary to take into account the qualitative nature of the solutions to the system of ordinary differential equations [46]. As an approximate solution, one can take a biunique function of the general analytical solution of a simplified system of differential equations [46]. Moreover, in the case of an autonomous system of differential equations, such a solution will have the property of a group and be a solution of an autonomous system of differential equations [47]. Such an approximate solution can be composed with analytical solutions of local simplifications of the system of ordinary differential equations being solved [46] by taking a biunique function from this solution composed of pieces [46].
It should be noted that in the absence of external flows, the system tends to an equilibrium state, and from any initial state (due to the zeroth law of thermodynamics) [34,36,37]. In the case of the presence of external flows, the system can evolve either into a stable stationary state or into an oscillatory regime (self-oscillations, forced oscillations, dynamic chaos, etc.) [48,49]. Oscillations in the system and their occurrence in the considered case are explained by the tendency of the system to a stationary state, which can change as a result of feedback [49]. Local simplifications of the equation’s parameters of the mathematical prototyping method can be piecewise constant or piecewise polynomial [50]. Moreover, these simplifications can be chosen by taking the piecewise constant dissipative matrix, balance matrix, quantities U , and external flows [48,49]. In such local areas, the system tends to a stationary state, which can also change as a result of a change in the balance parameters (for example, the total mass of the system, the total energy of the system, the total momentum of the system, etc.) [49]. Yet upon transition to another region, the stationary state can also change as a result of changes in the properties of the system [49]. Thus, an oscillatory motion of the system arises [49]. The general analytical solution of the system of equations of the mathematical prototyping method is obtained by stitching (without the refinement of the stitching method) simplified analytical solutions in local areas [46], which, in the limit, as the size of the areas tends to zero, converges to the analytical solution of the main system of Equations (1) and (2) [46].
The simplified analytic solutions x ˜ t = x ˜ t x ˜ 0 , t have the group property [47]:
x ˜ 0 = x ˜ t x ˜ 0 , 0 ,   x ˜ t x ˜ t x ˜ 0 , τ , t = x ˜ t x ˜ 0 , t + τ .
Approximate general solution of Equations (1) and (2) x ¯ t = x ¯ t x ¯ 0 , t 0 , t , where x ¯ 0 = x ¯ t x ¯ 0 , t 0 , t 0 , considering (27) is represented as:
x ¯ t x ¯ 0 , t 0 , t x ˜ t x ˜ i 1 0 , t t t i 1 t 0 , t i t 0 ,   x ˜ i 0 = x ˜ t x ˜ i 1 0 , t i ,   x ˜ 0 0 = x ¯ 0 ,   i = 1 , ,
where x ˜ i 0 , i = 1 , additively includes a random component.
If, in the space of system state coordinates, we find such a basis under which the expressions for the parameters of the mathematical prototyping method equations system (functions of the dissipative matrix, topology matrix, the scalar function of free energy or its partial derivatives) are simplified. In the new space of system state coordinates, the number of local areas into which the entire space should be divided decreases.
To transform the state coordinates, we introduce a reversible (if possible) by x ¯ function r x x ¯ , w , γ :
x = r x x ¯ , w , γ ,   w t = W U t , d x t d t ,
where W U t , d x t d t —is a functional that satisfies the condition:
U t c o n s t     t d x d t 0 w t = W U t , d x d t c o n s t
In this case, the expressions for the approximate solution of Equations (1) and (2) up to the coefficients γ found from these Equations (1) and (2) are represented as:
x t = r x x ¯ t x ¯ 0 , t 0 , t , w t , γ ,   w t = W U t , d x t d t .
Let us show that the analytical solution defined by (27)–(31) is the solution of Equations (1) and (2). Indeed, since x ¯ t = x ¯ t x ¯ 0 , t 0 , t is a solution to Equations (1) and (2) written for the topology matrix B ¯ x ¯ , U ¯ , dissipative matrix A ¯ x ¯ , U ¯ , and energy W ¯ x ¯ , U ¯ , as well as piecewise constant external flows d x ¯ t / d t , i.e.:
d x ¯ t d t = B ¯ x ¯ t , U ¯ t δ Δ x ¯ t d t + d x ¯ t d t ,   δ Δ x ¯ t d t = A ¯ x ¯ t , U ¯ t Δ F ¯ x ¯ t , U ¯ t , Δ F ¯ x ¯ , U ¯ = B ¯ T x ¯ , U ¯ F ¯ x ¯ , U ¯ ,   F ¯ x ¯ , U ¯ = x ¯ W ¯ x ¯ , U ¯ ,
and taking into account that:
x W ¯ x ¯ , U ¯ = x W ¯ r x 1 x , w , γ , U ¯ = J r , x ¯ T x ¯ , w , γ 1 x ¯ W ¯ x ¯ , U ¯
where J r , x ¯ x ¯ , w , γ is the Jacobian of matrix functions r x x ¯ , w , γ with respect to x ¯ , by virtue of (29) and (31) we have:
d x t d t = J r , x ¯ x ¯ t , w t , γ B ¯ x ¯ t , U ¯ t δ Δ x ¯ t d t + J r , x ¯ x ¯ t , w t , γ d x ¯ t d t + J r , w x ¯ t , w t , γ d w t d t   δ Δ x ¯ t d t = A ¯ x ¯ t , U ¯ t Δ F ¯ x ¯ t , U ¯ t Δ F ¯ x ¯ , U ¯ = B ¯ T x ¯ , U ¯ F ¯ x ¯ , U ¯ ,   F ¯ x ¯ , U ¯ = J r , x ¯ T x ¯ , w , γ x W ¯ r x 1 x , w , γ , U ¯ ,
where J r , w x ¯ , w , γ is the Jacobian matrix of the function r x x ¯ , w , γ with respect to w ; hence, taking into account that the parameters U ¯ t are piecewise constant U t , and hence U ¯ is a function U , having introduced the topology matrix B x , w , U :
B x , w , U = J r r x 1 x , w , γ , w , γ B ¯ r x 1 x , w , γ , U ,
external flows d x t / d t :
d x t d t = J r , x ¯ x ¯ t , w t , γ d x ¯ t d t + J r , w x ¯ t , w t , γ d w t d t ,
the free energy W x , w , U = W ¯ r x 1 x , w , γ , U ¯ and its partial derivatives F x , w , U by state coordinates x :
F x , w , U = x W ¯ r x 1 x , w , γ , U ¯ = x W x , w , U ,
as well as the dissipation matrix A x , w , U :
A x , w , U = A ¯ r x 1 x , w , γ , U ¯ ,
we have:
d x t d t = B x t , w t , U t δ Δ x ¯ t d t + d x t d t ,   δ Δ x ¯ t d t = A x t , w t , U t Δ F x t , w t , U t ,   Δ F x , w , U = Δ F ¯ r x 1 x , w , γ , U ¯ = B T x , w , U F x , w , U ,   F x , w , U = x W x , w , U .
This directly implies that the analytical solution given by (27)–(31) is the solution of Equations (1) and (2) of the mathematical prototyping method (because from (30) it can be seen that w can also be attributed to the parameters of U ). Thus it is possible to set such a positively defined (in the case of inertia—non-degenerate and non-negatively defined) dissipative matrix A x , w , U and such partial derivatives F x , w , U of the free energy W x , w , U (satisfying the condition of the total differential), as well as the topology matrix B x , w , U , that this general solution given by (27)–(31) will be the general solution of Equations (1) and (2) obtained for the mentioned quantities.
External energy flows in the system, in the general case, can either change the system balance parameters (for example, the internal energy and mass of the entire system), or not.
If conditions (30) are met:
  • The replacement of the state coordinates does not affect the fact that the balance parameters of the system are changed by external flows, it follows from Equations (32) and (33);
  • If the analytical solution in the old coordinate system tends to the stationary state from any initial state, then the solution in the new coordinate system also tends to the stationary state from any initial state, it follows from Equation (31);
  • If the analytical solution in the old coordinate system satisfies the group condition, then in the new coordinate system the solution also satisfies the group condition [47].
This implies the correctness of the general analytical solution (27)–(31) of the equations of the mathematical prototyping method.
Let us replace the state coordinates y ¯ and z ¯ by ξ in (13)–(17) in the same way as the replacement of the state coordinates in Equations (1) and (2) was carried out:
y ¯ ξ 0 , t 0 , t , γ = r y ¯ ξ ξ 0 , t 0 , t , γ ,   z ¯ ξ 0 , t 0 , t , γ = r z ¯ ξ ξ 0 , t 0 , t , γ ,
where the system of functions r y ¯ ξ , γ and r z ¯ ξ , γ is also resolvable with respect to ξ ; ξ ξ 0 , t 0 , t , being a general analytical solution of a piecewise simplified system of Equations (13)–(17), is defined similarly to (27) and (28):
ξ ξ 0 , t 0 , t ξ ˜ ξ ˜ i 1 0 , t t t i 1 t 0 , t i t 0 ,   ξ ˜ i 0 = ξ ˜ ξ ˜ i 1 0 , t i ,   ξ ˜ 0 0 = ξ 0 ,   i = 1 , ,
where in x ˜ i 0 , i = 1 , the random component is additively included, and ξ ˜ ξ ˜ 0 , t satisfies the property of the group [47]:
ξ ˜ 0 = ξ ˜ ξ ˜ 0 , 0 ,   ξ ˜ ξ ˜ ξ ˜ 0 , τ , t = ξ ˜ ξ ˜ 0 , t + τ .
Hence, similarly as described above, the analytical approximation of the general solution of the equations system (13)–(17) given in the form (34)–(36) is correct (follows from [47]).
Thus, the solution of the transformed Equations (13)–(17) of the mathematical prototyping method of energy processes consists in such a selection of constant parameters γ , that residuals e y ¯ γ , ξ 0 , t 0 , t and e z ¯ γ , ξ 0 , t 0 , t , having the meaning of fluctuations [51,52], determined based on (13)–(17) [51,52]:
e y ¯ γ , ξ 0 , t 0 , t = y ¯ ξ 0 , t 0 , t , γ t B ˜ y y ¯ ξ 0 , t 0 , t , γ , z ¯ ξ 0 , t 0 , t , γ , U t δ Δ x t d t d y ¯ t d t ,
e z ¯ γ , ξ 0 , t 0 , t = z ¯ ξ 0 , t 0 , t , γ t B ˜ z y ¯ ξ 0 , t 0 , t , γ , z ¯ ξ 0 , t 0 , t , γ , U t δ Δ x t d t d z ¯ t d t .
should not exceed in modulus a certain fraction of the mean maximum of fluctuations [51,52]. The quantities included in (37) and (38) are determined in accordance with the expressions:
d y ¯ t d t = C y , x y ¯ ξ 0 , t 0 , t , γ , z ¯ ξ 0 , t 0 , t , γ , U t d x t d t + C y , U y ¯ ξ 0 , t 0 , t , γ , z ¯ ξ 0 , t 0 , t , γ , U t d U t d t
d z ¯ t d t = C z , x y ¯ ξ 0 , t 0 , t , γ , z ¯ ξ 0 , t 0 , t , γ , U t d x t d t + C z , U y ¯ ξ 0 , t 0 , t , γ , z ¯ ξ 0 , t 0 , t , γ , U t d U t d t
δ Δ x t d t = A y ¯ ξ 0 , t 0 , t , γ , z ¯ ξ 0 , t 0 , t , γ , U t Δ F y ¯ ξ 0 , t 0 , t , γ , z ¯ ξ 0 , t 0 , t , γ , U t
Because by selecting parameters, it is necessary to achieve values of residuals e y ¯ γ , ξ 0 , t 0 , t and e z ¯ γ , ξ 0 , t 0 , t , not exceeding in modulus a certain fraction of the average maximum of fluctuations [51,52], and also to simultaneously minimize the objective function Q , determined by (24) for y t , determined on the basis of (11) and (18), due to:
y ξ 0 , t 0 , t , γ = y ¯ T ξ 0 , t 0 , t , γ y ˜ T ξ 0 , t 0 , t , γ T ,   y ˜ ξ 0 , t 0 , t , γ = g ˜ y y ¯ ξ 0 , t 0 , t , γ , z ¯ ξ 0 , t 0 , t , γ , U t ,
then it is advisable in terms of parameters γ and constant coefficients included in (25) and (26) to minimize the objective function Q ˜ [10,51,52]:
Q ˜ = Q y t + 1 2 i = 1 n l = 1 N t , i e y ¯ T γ i , ξ 0 , i , t 0 , i , t l , i L y ¯ t l , i e y ¯ γ i , ξ 0 , i , t 0 , i , t l , i + e z ¯ T γ i , ξ 0 , i , t 0 , i , t l , i L z ¯ t l , i e z ¯ γ i , ξ 0 , i , t 0 , i , t l , i ,
where L y ¯ t , L z ¯ t are positively defined symmetric matrices, because in this case, the residuals e y ¯ γ , ξ 0 , t 0 , t and e z ¯ γ , ξ 0 , t 0 , t , and also y t y Э t determine both the proximity of the exact solution of the system (13)–(17) to the approximate one, and the proximity of the calculated values of the characteristics of the measured parameters y t to the corresponding measured values y Э t of these parameters y t [10,51,53]. According to (24), Equation (43) will take the form:
Q ˜ = 1 2 i = 1 n j = 1 N t , i y i t j y i E t j T L ˜ i y i t j y i E t j + + 1 2 i = 1 n l = 1 N t , i e y ¯ T γ i , ξ 0 , i , t 0 , i , t l , i L y ¯ t l , i e y ¯ γ i , ξ 0 , i , t 0 , i , t l , i + e z ¯ T γ i , ξ 0 , i , t 0 , i , t l , i L z ¯ t l , i e z ¯ γ i , ξ 0 , i , t 0 , i , t l , i .
Thus, the construction of a model of physical and chemical processes in RES is carried out by minimizing the objective function Q ˜ , determined by virtue of (16), (17), (25), (26), (34)–(42), (44) in constant parameters ξ 0 , i , γ i , i = 1 , n (at fixed initial times t 0 , i , i = 1 , n ) and constant parameters included in (25) and (26) [40]. The number i of modes can include both control and operational modes of operation. After optimizing the objective function, Q ˜ we determine the controlled parameters z in the form:
z ξ 0 , t 0 , t , γ = z ¯ T ξ 0 , t 0 , t , γ z ˜ T ξ 0 , t 0 , t , γ T ,   z ˜ ξ 0 , t 0 , t , γ = g ˜ z y ¯ ξ 0 , t 0 , t , γ , z ¯ ξ 0 , t 0 , t , γ , U t .
The target function Q ˜ includes complex expressions (25) and (26), so it is proposed to convert them to a piecewise simplified form (for example, piecewise constant). Similarly, to simplify the objective function Q ˜ it is advisable to simplify the analytical expressions (34) of the general solution of the system of Equations (13)–(17) [40,45,54].
Thus, a generalized method of mathematical prototyping of energy processes is proposed to use as a unified approach to designing models of physical and chemical systems, in the application of which the following methods are used:
  • Transformations of the coordinate system of the state space to simplify the expressions of the state functions of the properties of substances and processes;
  • Transformations of the equations of the mathematical prototyping method with respect to the measured and controlled parameters;
  • Splitting the space of the system state coordinate into regions in order to obtain piecewise simplified state functions for the properties of substances and processes;
  • Reduction of the procedure for obtaining analytical solutions to the equations of the mathematical prototyping method for state coordinates to the problem of finding a global minimum.

3. Algorithm of Obtaining the Model from Experimental Data

3.1. The Sequence of Obtaining Controlled Parameters

So, the determination of the controlled parameters of the system by its measured parameters consists of the following sequence of actions:
For the experimentally measured parameters y of the RES, we perform the minimization of the objective function Q ˜ , determined by virtue of (16), (17), (25), (26), (34)–(42), (44), with respect to the parameters ξ 0 , i , γ i , i = 1 , n , and the parameters included in (25) and (26).
For the optimized parameters ξ 0 , i , γ i , and the parameters included in (25) and (26), we compare the values of residuals e y ¯ γ i , ξ 0 , i , t 0 , i , t and e z ¯ γ i , ξ 0 , i , t 0 , i , t with the value of the average maximum of fluctuations; if they are higher than the average maximum of fluctuations we should conduct the correction of r y ¯ ξ , γ and r z ¯ ξ , γ then return to point 1; otherwise, go to the next item.
Substituting the optimal parameters ξ 0 , i and γ i into (42) and (45) we obtain the dynamics of the measured y t = y ξ 0 , t 0 , t , γ and control z t = z ξ 0 , t 0 , t , γ parameters of the system.
Equations (13)–(18), and hence the proposed method for determining the controlled parameters of the system from the measured ones, can be implemented based on model-based design [55] using a block diagram [37].

3.2. Implementation of the Algorithm for Obtaining Controlled Parameters

The practical implementation of the approach proposed in the article to building models of physical and chemical processes in RES (the method of mathematical prototyping of energy processes) involves the following sequence of actions:
Formation of a list of physical and chemical processes in the object of study.
Determination of a set of assumptions, including considered modes of operation and external disturbances.
Writing a complete system of equations for the dynamics of physical and chemical processes in accordance with the proposed method—the method of mathematical prototyping.
Determining the required set of experimental data based on analyzing the system state dependence of the properties of substances and processes.
Setting analytical expressions (25) and (26) for the properties of the substances and the processes up to constant coefficients.
Piecewise simplification of analytical expressions (25) and (26).
Identification of the coefficients of analytical expressions.
Checking the correctness of piecewise simplification procedures.
Validation of the obtained models according to the control experimental data.

4. Results

In the current chapter, two examples of the application described methods are presented. The first of them is the calculation of physical and chemical processes model parameters in a nickel–cadmium battery from a 20NKBN-25-U3 series of three cells. The second one is a mathematical model for the voltage and temperature of lithium-ion batteries of the US18650VTC6 series.

4.1. An Example of a Nickel–Cadmium Battery

4.1.1. Physical and Chemical Processes in a Nickel–Cadmium Accumulator

The electrochemical system of a nickel–cadmium accumulator is a positive nickel oxide electrode N i O O H and a negative cadmium electrode immersed in an electrolyte solution K O H . The substance K O H does not enter the reaction, it is only a carrier of hydroxide ions O H [56].
The main current-generating process taking place on the positive nickel oxide electrode [56]:
N i O O H + H 2 O + e N i O H 2 + O H
Hydroxide ions O H diffuse through the electrolyte and react on the negative cadmium electrode in accordance with the reaction [56]:
C d + 2 O H C d O H 2 + 2 e
When a nickel–cadmium accumulator is recharged, the process of oxygen evolution proceeds at the positive electrode [56]:
2 O H 1 2 O 2 + H 2 O + 2 e
Oxygen O 2 diffuses through the porous separator to the negative electrode and is reduced on it [56]:
1 2 O 2 + C d + H 2 O C d O H 2
Reaction (49) of oxygen reduction is exothermic, it leads to heating of the accumulator, which can cause its thermal runaway [56].

4.1.2. Assumptions Imposed on the Mathematical Model of Physicochemical Processes in the Nickel–Cadmium Accumulators

The modeling of physicochemical processes in the nickel–cadmium accumulators is carried out taking into account the following assumptions [57]:
  • Aging processes in the accumulators are not modeled (due to the fact that they proceed much more slowly than main processes);
  • Hydrogen release is not simulated (the nickel–cadmium accumulators are fairly new) [56];
  • Distribution of water in the volume of the separator is even;
  • The volume of the separator is divided into near-anode and near-cathode regions, and the state of each region is characterized by averaged values of the distributed quantities;
  • Physical and chemical processes between each pair of electrodes are identical, therefore, the accumulators are represented with one pair of electrodes, on which the above processes occur;
  • The temperature of the accumulators is uniform;
  • Cross-diffusion of hydroxide ions O H and oxygen molecules is absent;
  • The oxygen above the separator is in equilibrium with the oxygen in the separator;
  • The contact of the base of the positive and negative electrodes with their spattering is ideal;
  • The “memory effect” [56] of a nickel–cadmium accumulators is not taken into account;
  • Capacities of double layers of positive and negative electrodes are not taken into account.
Within the framework of these assumptions, the model of a nickel–cadmium battery is built using the mathematical prototyping methods of energy processes [57].

4.1.3. Mathematical Model of Physical and Chemical Processes in a Nickel–Cadmium Accumulator

The mathematical model of a nickel–cadmium accumulator includes the following quantities [57,58,59]:
  • Charge transferred through an external circuit, Δ q , A∙h;
  • Current in the external circuit, I = δ Δ q / d t , A;
  • Current through the membrane δ Δ q m / d t , A
  • Component of the current in the external circuit, due to the main current-generating processes, δ Δ q b a s e / d t , A;
  • Charge accumulated in the membrane q a c c , A∙h;
  • Component of the current in the external circuit, due to the release of oxygen, δ Δ q O 2 / d t , A;
  • Membrane resistance R m , Ohm;
  • Membrane capacity C m , F;
  • Internal resistance due to the main current-generating processes r i n . b a s e , Ohm;
  • Internal resistance due to the release of oxygen r i n . O 2 , Ohm;
  • EMF of the main current-generating processes E b a s e , V;
  • EMF due to the release of oxygen E O 2 , V;
  • Cross EMF for the main current, due to the release of oxygen E O 2 b a s e , V;
  • The coefficient of cross-EMF for the main current, due to the release of oxygen ε O 2 b a s e , V;
  • Cross EMF for the current associated with the release of oxygen, due to the main current, E b a s e O 2 , V;
  • Cross-EMF coefficient for the current associated with the release of oxygen, due to the main current, ε b a s e O 2 , V;
  • Electrical potentials of the positive ϕ + and negative ϕ electrodes, respectively, V;
  • The number of moles of accumulated oxygen in the electrode regions of the positive Δ ν O 2 + and negative Δ ν O 2 electrodes, respectively;
  • The number of moles of oxygen released at the positive electrode δ Δ ν O 2 + c and utilized at the negative electrode δ Δ ν O 2 c ;
  • The number of moles of oxygen diffused through the membrane δ Δ ν O 2 ;
  • Chemical potentials of oxygen in the anode region μ O 2 + and cathode region μ O 2 , J/mol;
  • Equilibrium chemical potential of oxygen μ O 2 , J/mol;
  • Thermal coefficient for the main current q ˜ b a s e ;
  • Thermal coefficient for the main current associated with the release of oxygen at the positive electrode q ˜ O 2 + ;
  • Thermal coefficients of oxygen diffusion through the electrolyte membrane q ˜ O 2 d and oxygen utilization at the negative electrode q ˜ O 2 u ;
  • Heat capacity of the accumulator, C p , J/K;
  • Heat release power in the accumulator, Q V , W;
  • Heat transfer coefficient of the accumulator, K , W/(m2∙K);
  • Heat transfer area, S , m2;
  • Accumulator temperature, T , ambient temperature, T 0 , K.
The mentioned mathematical model of the physical and chemical processes dynamics in a nickel–cadmium accumulator, obtained by mathematical prototyping, has the following form [57,58,59]:
  • Stoichiometric ratios:
    δ Δ ν N i O H 2 = F δ Δ q b a s e ,   δ Δ ν O 2 + c = 1 4 F δ Δ q O 2 ,   δ Δ ν C d O H 2 = 1 2 F δ Δ q b a s e + δ Δ q O 2 + 2 δ Δ ν O 2 c ,
  • Equations of the conservation law of oxygen mass:
    δ Δ ν O 2 + = δ Δ ν O 2 + c δ Δ ν O 2 ,   δ Δ ν O 2 = δ Δ ν O 2 δ Δ ν O 2 c ;
  • Equivalent circuit equations for electrochemical processes (Figure 1):
    r i n . b a s e δ Δ q b a s e d t + R m δ Δ q m d t + ϕ + ϕ = E b a s e E O 2 b a s e ,   r b a s e . O 2 δ Δ q O 2 d t + R m δ Δ q m d t + ϕ + ϕ = E O 2 E b a s e O 2 ,   E O 2 = μ O 2 + μ O 2 4 F ,   E O 2 b a s e = ε O 2 b a s e δ Δ q O 2 d t ,   E b a s e O 2 = ε b a s e O 2 δ Δ q b a s e d t ,   ε O 2 b a s e = ε b a s e O 2 ,   I = δ Δ q d t = δ Δ q m d t + δ Δ q a c c d t = δ Δ q b a s e d t + δ Δ q O 2 d t ,   R m δ Δ q m d t = q a c c C m ;
  • Oxygen diffusion equations and oxygen utilization equations:
    δ Δ ν O 2 d t = D O 2 μ O 2 + μ O 2 ,   d Δ ν O 2 c d t = R O 2 μ O 2 μ O 2 ;
  • Heat release power:
    Q V = q ˜ b a s e r i n . b a s e d Δ q b a s e d t 2 + q ˜ O 2 + r i n . O 2 d Δ q O 2 d t 2 + R m d Δ q m d t 2 + q ˜ O 2 d 1 D O 2 d Δ ν O 2 d t 2 + + q ˜ b a s e ε O 2 b a s e + q ˜ O 2 + ε b a s e O 2 d Δ q O 2 d t d Δ q b a s e d t + q ˜ O 2 u R O 2 d Δ ν O 2 c d t 2 ;
  • Equations of thermal coefficients:
    q ˜ b a s e = 1 T E b a s e T E b a s e q a c c C m ϕ + ϕ ,   q ˜ O 2 + = 1 T E O 2 T E O 2 q a c c C m ϕ + ϕ ,   q ˜ O 2 d = 1 T μ O 2 + μ O 2 T μ O 2 + μ O 2 ,   q ˜ O 2 u = 1 T μ O 2 μ O 2 T μ O 2 μ O 2 ;
  • Heat balance equation:
    Q V = C p d T d t + K S T T 0
The resulting system of equations is closed, and we will be able to predict the physical and chemical processes in a nickel–cadmium accumulator [57,58,59] by knowing the following parameters:
  • The parameters of the equivalent circuit (Figure 1);
  • The function of the oxygen chemical potential;
  • The oxygen chemical potential at equilibrium;
  • Coefficients of diffusion and utilization of oxygen.

4.1.4. Identification of the Substances and Processes Properties in a Nickel–Cadmium Battery

The identification of the normalized variables is completed based on the following assumptions (and simplifications) [57,59]:
  • The oxygen cycle is assumed to be stationary—how much oxygen was released on the nickel oxide electrode, the same amount diffused to the cadmium electrode, the same amount was utilized on the electrode;
  • The diffusion coefficient of oxygen through the membrane is constant;
  • The main current-forming reactions (46) and (47) in the forward direction proceed on the regions of the corresponding electrodes free from the hydroxide film, in the reverse direction, on the regions covered with the hydroxide film (Figure 2);
  • The reaction of release (48) and utilization (49) of oxygen on the corresponding electrodes proceeds only on the areas of these electrodes free from the hydroxide film (Figure 2);
  • The area of the hydroxide film covering the electrodes is directly proportional to the number of moles of the corresponding hydroxides (Figure 2);
  • The reactivity coefficients of the main current-forming reactions do not depend on the number of moles of oxygen in the near-electrode region;
  • The capacity and resistance of the membrane do not depend on the current δ Δ q m / d t and the redistributed charge q a c c ;
  • The dependences of the reactivity coefficients of electrochemical reactions on the redistribution of the electrolyte in the active sites are similar to each other;
  • The parameters of a nickel–cadmium accumulator at temperatures below the critical temperature at which the development of thermal acceleration begins do not depend on temperature;
  • The heat capacity and heat transfer coefficient of the battery are constant throughout the charge/discharge process.
The assumptions formulated above related to the coating of the electrodes with a hydroxide film are taken into account when obtaining analytical dependences of the resistances r i n . b a s e , r i n . O 2 and the coefficient of transfer EMF ε b a s e O 2 = ε O 2 b a s e on the battery discharge level [57]. The coefficients of the transfer EMF do not depend on the degree of coating of the electrodes, and the resistances are r i n . b a s e , r i n . O 2 inversely proportional to the active areas of the electrodes, which means they are inversely proportional to the numbers of moles of hydroxide films covering the electrodes [57]. Hence, the analytical functions of the properties of the main current-forming reactions (46) and (47), as well as the oxygen cycle reactions (48) and (49), it is advisable to set in the form (depending on the charge or discharge of a nickel–cadmium battery) [57]:
r i n . b a s e = r i n . b a s e d + Δ r i n . b a s e d 2 1 C 0 + + Δ q + 4 F Δ ν O 2 + c C + + r i n . b a s e d Δ r i n . b a s e d 2 1 C 0 + Δ q + 4 F Δ ν O 2 c C ,    d Δ q b a s e d t 0 d Δ q b a s e d t + d Δ q O 2 d t 0 r i n . b a s e d + Δ r i n . b a s e d 2 1 C 0 + + Δ q + 4 F Δ ν O 2 + c C + + r i n . b a s e c Δ r i n . b a s e c 2 C 0 + Δ q + 4 F Δ ν O 2 c C ,    d Δ q b a s e d t 0 d Δ q b a s e d t + d Δ q O 2 d t < 0 r i n . b a s e c + Δ r i n . b a s e c 2 C 0 + + Δ q + 4 F Δ ν O 2 + c C + + r i n . b a s e d Δ r i n . b a s e d 2 1 C 0 + Δ q + 4 F Δ ν O 2 c C ,    d Δ q b a s e d t < 0 d Δ q b a s e d t + d Δ q O 2 d t 0 r i n . b a s e c + Δ r i n . b a s e c 2 C 0 + + Δ q + 4 F Δ ν O 2 + c C + + r i n . b a s e c Δ r i n . b a s e c 2 C 0 + Δ q + 4 F Δ ν O 2 c C ,    d Δ q b a s e d t < 0 d Δ q b a s e d t + d Δ q O 2 d t < 0
r i n . O 2 = 2 r i n . O 2 d r i n . b a s e d + Δ r i n . b a s e d 2 1 C 0 + + Δ q + 4 F Δ ν O 2 + c C + + r i n . b a s e d Δ r i n . b a s e d 2 1 C 0 + Δ q + 4 F Δ ν O 2 c C ,    d Δ q O 2 d t 0 d Δ q b a s e d t + d Δ q O 2 d t 0 2 r i n . O 2 d r i n . b a s e d + Δ r i n . b a s e d 2 1 C 0 + + Δ q + 4 F Δ ν O 2 + c C + + r i n . b a s e c Δ r i n . b a s e c 2 C 0 + Δ q + 4 F Δ ν O 2 c C ,    d Δ q O 2 d t 0 d Δ q b a s e d t + d Δ q O 2 d t < 0 2 r i n . O 2 c r i n . b a s e c + Δ r i n . b a s e c 2 1 C 0 + + Δ q + 4 F Δ ν O 2 + c C + + r i n . b a s e d Δ r i n . b a s e d 2 1 C 0 + Δ q + 4 F Δ ν O 2 c C ,    d Δ q O 2 d t < 0 d Δ q b a s e d t + d Δ q O 2 d t 0 2 r i n . O 2 c r i n . b a s e c + Δ r i n . b a s e c 2 1 C 0 + + Δ q + 4 F Δ ν O 2 + c C + + r i n . b a s e c Δ r i n . b a s e c 2 C 0 + Δ q + 4 F Δ ν O 2 c C ,    d Δ q O 2 d t < 0 d Δ q b a s e d t + d Δ q O 2 d t < 0
ε O 2 b a s e = ε b a s e O 2 = r i n . b a s e d Δ r i n . b a s e d 2 1 C 0 + Δ q + 4 F Δ ν O 2 c C , d Δ q b a s e d t + d Δ q O 2 d t 0 r i n . b a s e c Δ r i n . b a s e c 2 C; 0 + Δ q + 4 F Δ ν O 2 c C , d Δ q b a s e d t + d Δ q O 2 d t < 0
R O 2 = R O 2 1 C 0 + Δ q + 4 F Δ ν O 2 c C
The parameters included in the above dependencies are taken only for pure nickel oxide and cadmium electrodes [57].
When identifying the parameters included in (50)–(53), the discharge curve (the space of state coordinates) of a nickel–cadmium accumulator is represented as three sections (Figure 3) [57]:
  • Section I (Figure 3b) corresponds to a decrease in the discharge voltage associated with the polarization (redistribution) of the electrolyte;
  • In sections II and III (Figure 3b), the battery is discharged solely by filling the electrodes with products of the corresponding electrochemical reactions (hydroxide films).
  • Section II differs from section III of the discharge curve in Figure 3b. In section II, an electrode with a smaller capacity is slightly coated with a hydroxide film, while in section III, the electrode is already significantly coated with a hydroxide film [57]. In section III, each fraction of the remaining small free area is already more significant than in section II, therefore, the voltage decreases much faster, and the rate of decrease increase [57].
Identification of the properties of substances and processes in a nickel–cadmium battery is carried out from the equations of physical and chemical processes by performing the following actions [57,59]:
  • The resistances of the free sections of the electrodes and the EMF of the electrode reactions are determined from the condition of coincidence of the calculated voltage curve with the experimental one for the considered discharge current and ambient temperature in sections II and III (Figure 3) corresponding to the steady-state concentrations of the electrolyte.
  • The capacity of the membrane is determined from the condition of coincidence of the calculated voltage curve with the experimental one for the discharge current and ambient temperature in section I (Figure 3).
  • The distribution of electrolyte concentrations over near-electrode regions, currents through the membrane during polarization, and the heat release power in a nickel–cadmium battery are determined for the specified discharge current of a nickel–cadmium battery and ambient temperature.
  • The heat capacity and heat transfer coefficient of a nickel–cadmium battery are determined from the condition of coincidence of the calculated temperature dynamics with the experimental one.
  • The following dependences are constructed:
    • The membrane capacity on the current through the membrane and the temperature of the battery;
    • The resistances of the active sections of the electrode on the charge or discharge currents;
    • The EMF of the double layers on the redistribution of electrolyte concentrations.
  • Analytical expressions of the properties of substances and processes in a nickel–cadmium battery are constructed by adding additional components in (50)–(53).
  • The additional coefficients are determined from the discharge voltage and temperature dynamics for different discharge modes of a nickel–cadmium battery.
Thus, the identification of the parameters of the nickel–cadmium battery model is carried out using a piecewise analytical solution of the equations of dynamics of physical and chemical processes in the battery [57].
As an example, let us consider the calculation of physical and chemical processes model parameters in a nickel–cadmium battery from a 20NKBN-25-U3 series 3 battery. The discharge curves (Figure 4) show the calculated curves for the identified values for sections II and III of internal resistances and EMF without taking into account polarization. In addition, the graphs (Figure 5) show the calculated curves, which additionally take into account the membrane capacity identified in section I. The coincidence of the calculated and experimental data, as can be seen in Figure 5, confirms the sufficiently high accuracy of the calculation of the proposed method of mathematical prototyping [57].
For different discharge currents, the dependence of the resistances of pure electrodes is shown in Figure 6 [57]. As it is easy to see from Figure 6, such resistances fall with increasing discharge current, which corresponds to the theoretical provisions of electrochemistry. Additionally, with an increase in the discharge current, the capacity of the membrane decreases (Figure 7) [57]. Such dependencies can be considered as dependencies on the corresponding currents since these currents in the steady state are equal to discharge currents.
Taking into account the dependencies shown in Figure 6 and Figure 7 in (50)–(53), as well as adding and identifying additional components in (50)–(53) (similarly to (25) and (26)), we will obtain a more complete model of physicochemical processes in a nickel–cadmium battery.
Such a model makes it possible, for example, to predict the temperature of a nickel–cadmium battery, thereby predicting and taking measures to prevent the thermal runaway of the battery [57,59].

4.2. Lithium-Ion Accumulator Simulation

The principle of operation of lithium-ion batteries is the intercalation/deintercalation of lithium ions into the electrodes (Figure 8) [22,60]. Similarly, to a nickel–cadmium battery, as the cells in the electrode are filled with lithium ions, the active area of the electrodes decreases. Based on this, in [61], a mathematical model for the voltage of a lithium-ion battery of the QL079KM series was obtained from the discharge voltage curves.
A mathematical model for the voltage and temperature of lithium-ion batteries of the US18650VTC6 series was obtained in [61,62] from test results similar to a nickel–cadmium battery. The largest relative error of this model did not exceed 12% [62]. Analytical expressions (25) and (26) in [62] were given by further modification of the Balter–Vollmer model [22,63].

5. Discussion

In this paper, a unified approach to designing models of physical and chemical processes is proposed—a generalized method of mathematical prototyping of energy processes.
The model obtained by mathematical prototyping incorporates the laws of thermodynamics, conservation laws, as well as some physical features of the processes in a RES. This guarantees the correctness of the desired model, i.e., its consistency with its general physical laws.
Obtaining a model of an arbitrary system is reduced to specifying classes of analytical expressions of properties of substances and processes satisfying the corresponding constraints with accuracy up to experimentally determined constant coefficients. The specified classes of analytical expressions cover the entire space of functions of the properties of substances and processes of the system, taking into account restrictions. The determination of the controlled system’s parameters by the measured parameters is reduced to the identification of constant coefficients of the specified analytical expressions from experimental data.
To simplify the integration of equations of the mathematical prototyping method in order to reduce computational costs, an analytical task for solving differential equations of the mathematical prototyping method with an accuracy of constant coefficients is proposed.
The proposed approach to piecewise simplification of analytical expressions also reduces computational costs. Additionally, the possibility of parallelization of calculations that appears as a result of piecewise analytical simplification significantly speeds up the execution of calculations.
The correctness of the models obtained by the proposed method of mathematical prototyping for specific instances of objects allows using this approach for:
  • Formation and refinement of real-time digital twins of objects and systems;
  • Synthesis of objects and systems governing laws;
  • Diagnostics and forecasting of the technical condition of systems, as well as medical diagnostics;
  • To form and optimize technological processes (in operation and maintenance, biochemistry and bioengineering, geoengineering and meteorology, aerospace technologies, etc.);
  • Designing new facilities and systems.
The physicality of the proposed method of mathematical prototyping allows using artificial intelligence methods to obtain models that do not contradict physics, unlike classical methods of constructing simulation models. The proposed piecewise analytical approach allows the processing of experimental data, and hence the training of models, in parts.

Author Contributions

Conceptualization, S.K.; Methodology, S.K. and I.S.; Software, I.S.; Formal analysis, I.A.; Investigation, I.S.; Data curation, I.A.; Writing—original draft, I.S.; Writing—review & editing, S.K. and I.A.; Supervision, S.K. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Norenkov, I.P. Fundamentals of Computer-Aided Design, 4th ed.; Publishing House of Bauman Moscow State Technical University: Moscow, Russia, 2009; pp. 12–40. [Google Scholar]
  2. Barzilovich, E.Y. Models of Maintenance of Complex Systems; Higher School: Moscow, Russia, 1982; pp. 5–74. [Google Scholar]
  3. Kolodezhny, L.P.; Chernodarov, A.V. Reliability and Technical Diagnostics; VVA Publishing House named after N.E. Zhukovsky and Yu.A. Gagarin: Moscow, Russia, 2010; pp. 7–178. [Google Scholar]
  4. Bessekersky, V.A.; Popov, E.P. Theory of Automatic Control Systems, 3rd ed.; Nauka: Moscow, Russia, 1975; pp. 102–113. [Google Scholar]
  5. Akimov, V.A. Scientific foundations of the general security theory. Civ. Secur. Technol. 2017, 14, 4–9. [Google Scholar]
  6. Zhilina, I.Y. Geoengineering as a way to combat climate change: Benefit or harm? Social and humanitarian sciences: Domestic and foreign literature. Econ. Abstr. J. 2020, 106–115. [Google Scholar]
  7. Khromov, S.P.; Petrosyants, M.A. Meteorology and Climatology, 7th ed.; Moscow University Press: Moscow, Russia, 2006; pp. 19–25. [Google Scholar]
  8. Gaydes, M.A. General Theory of Systems (Systems and System Analysis), 2nd ed.; Globus-Press: Moscow, Russia, 2005; pp. 24–168. [Google Scholar]
  9. Antonov, A.V. System Analysis; Publishing House “Higher School”: Moscow, Russia, 2004; pp. 37–124. [Google Scholar]
  10. Eykhoff, P. Systems Identification: Parametrs and State Estimation; University of Technology: Eindhoven, The Netherlands, 1975; pp. 17–49. [Google Scholar]
  11. Tsypkin, Y.Z. Information Theory of Identification; Science, Physical Education: Moscow, Russia, 1995; pp. 13–57. [Google Scholar]
  12. Vyugin, V.V. Mathematical Foundations of Machine Learning and Forecasting; ICNMO Publishing House: Moscow, Russia, 2014; pp. 16–66. [Google Scholar]
  13. Flach, P. Machine Learning. The Art and Science of Algorithms that Make Sense of Data; Cambridge University Press: Cambridge, UK, 2015; pp. 14–61. [Google Scholar]
  14. Haykin, S. Neural Networks. A Comprehensive Foundation, 2nd ed.; Prentice Hall: Upper Saddle River, NJ, USA, 2006; pp. 42–161. [Google Scholar]
  15. Strizhov, V.V. Methods of Inductive Generation of Regression Models; Computing Center of the Russian Academy of Sciences: Moscow, Russia, 2008; pp. 24–56. [Google Scholar]
  16. Goemans, M.X.; Williamson, D.P. A general approximation technique for constrained forest problems. SIAM J. Comput. 1995, 24, 296–317. [Google Scholar] [CrossRef] [Green Version]
  17. Hegde, C.; Indyk, P.; Schmidt, L. A nearly-linear time framework for graph-structured sparsity. In Proceedings of the International Conference on Machine Learning, Lille, France, 6–11 July 2015; pp. 928–937. [Google Scholar]
  18. Diveev, A.I. Properties of the superposition of functions for numerical methods of symbolic regression. Cloud Sci. 2016, 3, 290–301. [Google Scholar]
  19. Dang Thi Fook; Diveev, A.I.; Sofronova, E.A. Solving problems of identification of mathematical models of objects and processes by the method of symbolic regression. Cloud Sci. 2018, 5, 147–162. [Google Scholar]
  20. Diveev, A.I.; Lomakova, E.M. The method of binary genetic programming for the search for a mathematical expression. In Bulletin of the Peoples’ Friendship University of Russia; Engineering Research: Moscow, Russia, 2017; Volume 18, pp. 125–134. [Google Scholar] [CrossRef] [Green Version]
  21. Dzyadyk, V.K. Introduction to the Theory of Uniform Approximation of Functions by Polynomials; Nauka Publishing House, Main Editorial Office of Physical and Mathematical Literature: Moscow, Russia, 1977; pp. 9–202. [Google Scholar]
  22. Borisevich, A.V. Modeling of Lithium-Ion Batteries for Battery Management Systems: An Overview of the Current State. Electron. Sci. Pract. J. Mod. Eng. Technol. 2014. Available online: http://technology.snauka.ru/2014/05/3542 (accessed on 23 October 2022).
  23. Wang, Y.; Tian, J.; Sun, Z.; Wang, L.; Xu, R.; Li, M.; Chen, Z. A comprehensive review of battery modeling and state estimation approaches for advanced battery management systems. Renew. Sustain. Energy Rev. 2020, 131, 110015. [Google Scholar] [CrossRef]
  24. Etkin, V.A. Energodynamics (Synthesis of Energy Transfer and Conversion Theories); SPb.: Nauka, Russia, 2008; 409p, ISBN 978-5-02-025318-6. [Google Scholar]
  25. Etkin, V.A. Towards a unified thermodynamic theory of real transfer processes. Inf. Process. Syst. Technol. 2021, 9–18. [Google Scholar]
  26. Etkin, V.A. Ergodynamic theory of the evolution of biological systems. Inf. Process. Syst. Technol. 2022, 12–24. [Google Scholar] [CrossRef]
  27. On the Fourth Principle of Thermodynamics. Available online: http://samlib.ru/e/etkin_w_a/ochetvernomnachaletermodynamiki.shtml (accessed on 24 October 2022).
  28. Zarubin, V.S.; Kuvyrkin, G.N. Mathematical Models of Mechanics and Electrodynamics of a Continuous Medium; Publishing House of Bauman Moscow State Technical University: Moscow, Russia, 2008; pp. 17–196. [Google Scholar]
  29. Zarubin, V.S. Mathematical Modeling in Engineering, 2nd ed.; Publishing House of Bauman Moscow State Technical University: Moscow, Russia, 2003; pp. 87–286. [Google Scholar]
  30. Belenky, I.M. Introduction to Analytical Mechanics; Higher School: Moscow, Russia, 1964; pp. 101–165. [Google Scholar]
  31. Gorshkov, A.G.; Rabinsky, L.N.; Tarlakovsky, D.V. Fundamentals of Tensor Analysis and Continuous Medium Mechanics; Publishing House Nauka: Moscow, Russia, 2000; pp. 86–213. [Google Scholar]
  32. Bessonov, L.A. Theoretical Foundations of Electrical Engineering; Higher School: Moscow, Russia, 1996; pp. 2–308, 404–554. [Google Scholar]
  33. Demirel, Y.; Gerbaud, V. Nonequilibrium Thermodynamics. Transport and Rate Processes in Physical, Chemical and Biological Systems, 3rd ed.; Elsevier: Amsterdam, the Netherlands, 2014; pp. 75–563. [Google Scholar]
  34. Krutov, V.I.; Isaev, S.I.; Kozhinov, I.A. Technical Thermodynamics, 3rd ed.; Higher School: Moscow, Russia, 1991; pp. 16–84, 206–277. [Google Scholar]
  35. Groot, S.R. Thermodynamics of Irreversible Processes; State Publishing House of Technical and Theoretical Literature: Moscow, Russia, 1956; pp. 19–233. [Google Scholar]
  36. Starostin, I.E.; Bykov, V.I. Kinetic Theorem of Modern Non-Equilibrim Thermodynamics; Open Science Publishing: Raley, NC, USA, 2017; pp. 132–206. [Google Scholar]
  37. Starostin, I.E.; Stepankin, A.G. Software Implementation of Modern Nonequilibrium Thermodynamics Methods. The Simulation System of Physico-Chemical Processes SimulationNonEqProcSS v.0.1.0; Lambert Academic Publishing: Beau Bassin, Mauritius, 2019; pp. 9–46. [Google Scholar]
  38. Starostin, I.E.; Khalyutin, S.P.; Bykov, V.I. Setting the State functions for the properties of substances and processes in a differential form. Complex Syst. 2022, 4–16. [Google Scholar]
  39. Khalyutin, S.P.; Starostin, I.E.; Davidov, A.O.; Kharkov, V.P.; Zhmurov, B.V. Digital twins in the theory and practice of aviation electric power industry. Electricity 2022, 4–13. [Google Scholar] [CrossRef]
  40. Starostin, I.E. Setting correct analytical expressions for transformed potentially streaming models. In Proceedings of the XIX International Scientific and Practical Conference “Innovative, Information and Communication Technologies”, Moscow, Russia, 21–25 November 2022; pp. 263–269. [Google Scholar]
  41. Gurov, A.A.; Badaev, F.Z.; Ovcharenko, L.P.; Shapoval, V.N. Chemistry, 3rd ed.; Publishing House of Bauman Moscow State Technical University: Moscow, Russia, 2007; pp. 175–592. [Google Scholar]
  42. Cybenko, G.V. Approximation by Superpositions of a Sigmoidal function. Math. Control. Signals Syst. 1989, 2, 303–314. [Google Scholar] [CrossRef]
  43. Gorban, A.N. Generalized approximation theorem and computational capabilities of artificial neural networks. Sib. Gournal Numer. Math. 1998, 1, 11–24. [Google Scholar]
  44. Kalitkin, N.N. Approximation of functions. In Numerical Methods; Nauka: Moscow, Russia, 1978; pp. 27–69. [Google Scholar]
  45. Chernorutsky, I.G. Optimization Methods; St. Petersburg State Polytechnic University Publishing House: St. Petersburg, Russia, 2012; pp. 5–46. [Google Scholar]
  46. Kalitkin, N.N. Numerical Methods; Nauka: Moscow, Russia, 1978; pp. 238–257. [Google Scholar]
  47. Arnold, V.I. Ordinary Differential Equations, 3rd ed.; ICNMO: Moscow, Russia, 2012; pp. 13–158. [Google Scholar]
  48. Prigozhin, I.; Kondepudi, D. Modern Thermodynamics: From Heat Engines to Dissipative Structures; Mir: Moscow, Russia, 2002; pp. 319–425. [Google Scholar]
  49. Bykov, V.I.; Starostin, I.E.; Khalyutin, S.P. Analysis of the formation of dissipative structures in complex concentrated systems based on the potentially streaming method. Cybernetic approach. Complex Syst. 2012, 72–89. [Google Scholar]
  50. Lantsov, V.N. Methods of Lowering the Order of Models of Complex Systems; Publishing House of Vladimir State University Named after A.G. and N.G. Stoletovs: Vladimir, Russia, 2017; pp. 16–75. [Google Scholar]
  51. Starostin, I.E. An algorithm for numerical integration of potential-flow equations in concentrated parameters with control of the correctness of the approximate solution. Com. Res. Model. 2014, 6, 479–493. [Google Scholar] [CrossRef] [Green Version]
  52. Starostin, I.E.; Bykov, V.I.; Stepankin, A.G. Numerical modeling of physical and physico-chemical processes taking into account fluctuations and checking the correctness of the approximate solution. Proc. Int. Symp. Reliab. Qual. 2018, 1, 79–83. [Google Scholar]
  53. Starostin, I.E.; Bykov, V.I. Estimation of the error tolerance of measuring the properties of substances and processes for modeling of nonequilibrium processes by the potential-flow method. Proc. Int. Symp. Reliab. Qual. 2015, 2, 203–206. [Google Scholar]
  54. Starostin, I.E.; Khalyutin, S.P.; Altukhov, A.V.; Davidov, A.O. Parallelization Applied to the Synthesis Methodology and Operation of Complex Systems Based on the Analysis and Modelling of their Physical and Chemical Processes. In Proceedings of the 2020 1st International Conference Problems of Informatics, Electronics, and Radio Engineering (PIERE), Novosibirsk, Russia, 10–11 December 2020; pp. 287–294. [Google Scholar]
  55. Vasiliev, S.N.; Yusupov, R.M. Simulation modeling. In Proceedings of the Theory and Practice: The Seventh All-Russian Sci-entific and Practical Conference, Moscow, Russia, 21–23 October 2015; V.A. Trapeznikov Institute of Management Problems of the Russian Academy of Sciences: Moscow, Russia, 2015; Volume 2, pp. 280–284. [Google Scholar]
  56. Taganova, A.A.; Bubnov, Y.I.; Orlov, S.B. Hermetic Chemical Current Sources: Cells and Accumulators. Equipment for Testing and Operation: Handbook; Chemical Publishing House: St-Petersburg, Russia, 2005; pp. 67–86. [Google Scholar]
  57. Khalyutin, S.P.; Starostin, I.E. Determination of the parameters of the substitution scheme of a potential-flow model of a nickel-cadmium battery by the method of hydroxide films. Reliab. Qual. Pr. International. Simp. (Penza, 2011) 2011, 318–324. [Google Scholar]
  58. Starostin, I.E.; Khalyutin, S.P. Modeling of physico-chemical processes in nickel-cadmium batteries by the potential-flow method. Innov. Based Inf. Commun. Technol. 2011, 517–520. [Google Scholar]
  59. Starostin, I.E. Determination of parameters of physico-chemical processes of nickel-cadmium batteries by the method of hydroxide films. Innov. Based Inf. Commun. Technol. 2011, 520–523. [Google Scholar]
  60. Kedrinsky, I.A.; Yakovlev, V.G. Lithium-ion Batteries; Platinum: Krasnoyarsk, Russia, 2002; pp. 19–22. [Google Scholar]
  61. Khalyutin, S.P.; Zhmurov, B.V.; Starostin, I.E. Mathematical modeling of electrochemical processes in lithium-ion batteries by the potential-flow method. Civ. Aviat. High Technol. 2014, 65–73. [Google Scholar]
  62. Starostin, I.E.; Khalyutin, S.P.; Davidov, A.O.; Punt, E.A.; Pavlova, V.I. Obtaining a Model for the Voltage and Temperature of the US18650VTC6 Series Lithium-ion Battery in Constant Current Discharge Mode from the Analysis of Physical and Chemical Processes in the Accumulator. In Proceedings of the 2021 XVIII Technical Scientific Conference on Aviation Dedicated to the Memory of N.E. Zhukovsky (TSCZh), Moscow, Russian, 29–30 October 2021; pp. 109–117. [Google Scholar] [CrossRef]
  63. Smith, K.; Wang, C.Y. Solid-state diffusion limitations on pulse operation of a lithium-ion cell for hybrid electric vehicles. J. Power Sources 2006, 161, 628–639. [Google Scholar] [CrossRef]
Figure 1. Equivalent circuit for electrochemical processes in a nickel–cadmium accumulator.
Figure 1. Equivalent circuit for electrochemical processes in a nickel–cadmium accumulator.
Energies 16 01933 g001
Figure 2. Active sites of nickel oxide (a) and cadmium (b) electrodes.
Figure 2. Active sites of nickel oxide (a) and cadmium (b) electrodes.
Energies 16 01933 g002
Figure 3. Characteristic view of discharge curves of nickel–cadmium batteries (a) and characteristic sections of the discharge curve of a nickel–cadmium battery (b).
Figure 3. Characteristic view of discharge curves of nickel–cadmium batteries (a) and characteristic sections of the discharge curve of a nickel–cadmium battery (b).
Energies 16 01933 g003
Figure 4. Calculation of internal resistances and EMF of a nickel–cadmium battery by comparing the calculated data with the experiment.
Figure 4. Calculation of internal resistances and EMF of a nickel–cadmium battery by comparing the calculated data with the experiment.
Energies 16 01933 g004
Figure 5. Calculation of the membrane capacities of a nickel–cadmium battery by comparing the calculated data with the experiment.
Figure 5. Calculation of the membrane capacities of a nickel–cadmium battery by comparing the calculated data with the experiment.
Energies 16 01933 g005
Figure 6. Dependence of internal resistance parameters on discharge current.
Figure 6. Dependence of internal resistance parameters on discharge current.
Energies 16 01933 g006
Figure 7. Dependence of the membrane capacity on the discharge current.
Figure 7. Dependence of the membrane capacity on the discharge current.
Energies 16 01933 g007
Figure 8. The principle of operation of a lithium-ion battery.
Figure 8. The principle of operation of a lithium-ion battery.
Energies 16 01933 g008
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Khalyutin, S.; Starostin, I.; Agafonkina, I. Generalized Method of Mathematical Prototyping of Energy Processes for Digital Twins Development. Energies 2023, 16, 1933. https://doi.org/10.3390/en16041933

AMA Style

Khalyutin S, Starostin I, Agafonkina I. Generalized Method of Mathematical Prototyping of Energy Processes for Digital Twins Development. Energies. 2023; 16(4):1933. https://doi.org/10.3390/en16041933

Chicago/Turabian Style

Khalyutin, Sergey, Igor Starostin, and Irina Agafonkina. 2023. "Generalized Method of Mathematical Prototyping of Energy Processes for Digital Twins Development" Energies 16, no. 4: 1933. https://doi.org/10.3390/en16041933

APA Style

Khalyutin, S., Starostin, I., & Agafonkina, I. (2023). Generalized Method of Mathematical Prototyping of Energy Processes for Digital Twins Development. Energies, 16(4), 1933. https://doi.org/10.3390/en16041933

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop