Synchronous machine electromagnetic transient modeling method, system and equipment based on saturation effect

Abstract

The embodiment of the present invention relates to a synchronous machine electromagnetic transient modeling method, system and equipment based on saturation effect. The method predicts the first rotor angular velocity, the first rotor angle, the first current q component, the first current d component, the first flux turning point d component, the first flux turning point q component and the flux turning data of the synchronous machine for analysis and processing to obtain the corresponding second current q component, the second current d component, the second rotor angular velocity, the second rotor angle, the stator flux d component and the stator flux q component. This not only avoids the occurrence of historical quantities and current quantities of the synchronous machine's rotating potential, but also improves the accuracy of the simulation calculation results. The calculated results can maintain the calculation efficiency of the dq0 model on the basis of achieving the accuracy of the phase domain model, thereby making the method have high simulation accuracy and fast calculation efficiency. The synchronous machine electromagnetic transient modeling method based on saturation effect can be suitable for the development of electromagnetic transient simulation software for power systems used in actual engineering calculations.

Description

Technical Field

[0001] The present invention relates to the field of electromagnetic transient technology, and in particular to a synchronous machine electromagnetic transient modeling method, system and equipment based on saturation effect. Background Art

[0002] With the rapid promotion and application of new energy, DC transmission, especially flexible DC transmission, large-scale grid electromagnetic transient simulation has become a new trend. How to significantly improve the simulation efficiency of electromagnetic transient models and algorithms while ensuring simulation accuracy is a topic that experts and scholars have been committed to studying.

[0003] As an important electrical component in electromagnetic transient simulation, the saturation characteristics of rotating motors have an impact on power flow, static and transient stability, and electromagnetic transients. The saturation characteristics of rotating motors are nonlinear, and it is very difficult to accurately simulate the saturation of motors. How to achieve efficient modeling and simulation of rotating motors considering saturation effects is crucial to the accuracy and efficiency of electromagnetic transient simulation of the entire system, especially the large-scale access system of new energy. Summary of the invention

[0004] The embodiments of the present invention provide a synchronous machine electromagnetic transient modeling method, system and device based on saturation effect, which are used to solve the technical problems of low simulation accuracy and low efficiency in implementing rotating motor modeling and simulation of saturation effect in existing electromagnetic transient simulation software.

[0005] In order to achieve the above purpose, the embodiment of the present invention provides the following technical solutions:

[0006] A synchronous machine electromagnetic transient modeling method based on saturation effect includes the following steps:

[0007] S1. Use linear extrapolation method to predict the first rotor angular velocity, the first rotor angle, the first current q component of the armature current, the first current d component of the armature current, the first flux turning point d component of the flux turning point and the first flux turning point q component of the flux turning point of the synchronous machine at a certain moment, and obtain the flux turning point data of the synchronous machine at the turning point between the unsaturated state and the saturated state;

[0008] S2. Processing the first current q component and the first current d component to obtain a first Norton circuit simulating a synchronous machine; converting the first Norton circuit from a dq0 quantity to a second Norton circuit in abc phasor form by coordinate transformation;

[0009] S3. Inverting the equivalent resistance matrix in the second Norton circuit to obtain an equivalent conductance matrix, and inputting the equivalent conductance matrix into the network conductance matrix for solution to obtain the three-phase voltage of the synchronous machine port;

[0010] S4. Processing the three-phase voltage to obtain the second current q component, the second current d component and the rotor current of the synchronous machine armature current; and processing the second current q component, the second current d component, the first flux turning point d component and the first flux turning point q component to obtain the stator flux d component, the stator flux q component and the air gap flux data of the synchronous machine;

[0011] S5. If the air gap flux data is not greater than the saturated air gap flux threshold, the second current q component, the second current d component, the stator flux d component and the stator flux q component are solved in the mechanical system equation to obtain the second rotor angular velocity and the second rotor angle of the synchronous machine; the stator flux d component, the stator flux q component and the flux turning data are also processed to obtain the second flux turning point d component and the second flux turning point q component of the synchronous machine flux turning point;

[0012] S6. Calculate the second current q component, the second current d component, the second rotor angular velocity, the second rotor angle, the second flux turning point d component and the second flux turning point q component respectively with the corresponding first current q component, the first current d component, the first rotor angular velocity, the first rotor angle, the first flux turning point d component and the first flux turning point q component to obtain the corresponding absolute values ​​of the errors; if all the absolute values ​​of the errors are less than the allowable error value, return to step S1.

[0013] Preferably, the synchronous machine electromagnetic transient modeling method based on saturation effect includes: if any one of the absolute values ​​of the error is not less than the allowable error value, returning to step S4.

[0014] Preferably, the first flux turning point d component, the first flux turning point q component and the flux turning data are processed to obtain a second flux turning point d component and a second flux turning point q component of the synchronous machine flux turning point, comprising:

[0015] Based on the parameter data of the synchronous machine, the first flux turning point d component, the first flux turning point q component, the second current q component and the second current d component, an air gap flux d component and an air gap flux q component are calculated, and an inverse trigonometric function is used to calculate the ratio of the air gap flux q component to the air gap flux d component to obtain a flux rotation angle at the flux turning point of the synchronous machine;

[0016] The flux turning point data and the flux turning angle are calculated using trigonometric functions to obtain a second flux turning point d component and a second flux turning point q component of the synchronous machine flux turning point.

[0017] Preferably, the synchronous machine electromagnetic transient modeling method based on saturation effect includes: if the air gap flux data is greater than the saturation air gap flux threshold, using a saturation correction parameter to correct the air gap flux d component, the air gap flux q component and the air gap flux data to obtain a modified stator flux d component and a stator flux q component.

[0018] Preferably, the mechanical system equation is:

[0019]

[0020]

[0021]

[0022] Where p is the number of pole pairs of the synchronous machine, λ q is the stator flux q component, λ d is the stator flux d component, is the second current d component, is the second current q component, J is the moment of inertia of the synchronous machine, D is the viscosity and air friction damping coefficient of the synchronous machine, T is the constant of the synchronous machine, ω is the second rotor angular velocity, θ is the second rotor angle, and t is the simulation time.

[0023] Preferably, the first current q component and the first current d component are processed to obtain a first Norton circuit simulating a synchronous machine; and the first Norton circuit is converted from a dq0 quantity to a second Norton circuit in abc phase form by coordinate transformation, comprising:

[0024] Obtaining the stator and rotor voltage equations of the synchronous machine, and performing discretization processing using an implicit trapezoidal integration method according to the stator and rotor voltage equations to obtain a first transformation equation;

[0025] Performing Park transformation on the first transformation equation, eliminating rotor variables, and using average resistance processing on the dq axis to obtain the Thevenin equation on the stator side;

[0026] The Thévenin equation on the stator side is converted into a first Norton circuit simulating a synchronous machine by mathematical transformation;

[0027] The first Norton circuit is converted from a dq0 quantity to a second Norton circuit of abc phase using a phase coordinate transformation formula;

[0028] Wherein, the first Norton circuit is:

[0029]

[0030] R ave =(R d +R q ) / 2;

[0031] The phase coordinate transformation formula is:

[0032]

[0033] In the formula, is the first current d component, is the first current q component, R d , R q R0 and R1 are the resistance parameters of the resistance matrix in the Thevenin equation on the stator side, e d 、e q , e0 are the voltage parameters of the voltage source matrix in the Thevenin equation on the stator side, i d，source is the first current d value of the first Norton circuit, i q，source is the second current q value of the first Norton circuit, i 0，source is the third current 0 value of the first Norton circuit, θ1 is the first rotor angle, i a，source is the first current of the a-phase current source of the second Norton circuit, i b，source is the second current of the second Norton circuit b-phase current source, i c，source is the third current of the second Norton circuit c-phase current source.

[0034] Preferably, the three-phase voltage is processed to obtain a second current q component, a second current d component and a rotor current of an armature current of a synchronous machine; and the second current q component, the second current d component, the first flux turning point d component and the first flux turning point q component are processed to obtain a stator flux d component, a stator flux q component and air gap flux data of the synchronous machine, including:

[0035] The three-phase voltage is converted by Park transformation to obtain a dq0-axis voltage component corresponding to the three-phase voltage;

[0036] According to the matrix parameters of the Thevenin equation on the stator side and the dq0-axis voltage component, the second current q component and the second current d component of the armature current of the synchronous machine are calculated by the armature current calculation formula;

[0037] Based on the parameter data of the synchronous machine, the dq0-axis voltage component, the second current q component and the second current d component, the rotor current of the synchronous machine is calculated by a rotor current calculation formula;

[0038] Based on the parameter data of the synchronous machine and the second current q component, the second current d component, the first flux turning point d component and the first flux turning point q component, the stator flux d component, the stator flux q component and the air gap flux data of the synchronous machine are calculated by a stator flux dq component calculation formula;

[0039] Wherein, the Parker transform is:

[0040]

[0041] The armature current calculation formula is:

[0042]

[0043] R ave =(R d +R q ) / 2;

[0044] The rotor current calculation formula is:

[0045]

[0046] i r =[i f i D i g i Q ] T

[0047]

[0048]

[0049]

[0050]

[0051]

[0052] The calculation formula of the stator flux dq component is:

[0053] λ d =λ ld +λ md ; q =λ lq +λ mq ;

[0054]

[0055]

[0056] In the formula, is the second current d component, is the second current q component, R d , R qR0 and R1 are the resistance parameters of the resistance matrix in the Thevenin equation on the stator side, e d 、e q , e0 are the voltage parameters of the voltage source matrix in the Thevenin equation on the stator side, θ1 is the first rotor angle, v a is the phase a voltage of the three-phase voltage, v b is the b-phase voltage of the three-phase voltage, v c is the c-phase voltage of the three-phase voltage, v d is the first voltage d component of the dq0 axis voltage component, v q is the second voltage q component of the dq0 axis voltage component, v0 is the third voltage 0 component of the dq0 axis voltage component, λ d is the stator flux d component, λ q is the stator flux q component, λ m is the air gap flux data, λ md is the air gap flux d component, λ mq is the air gap flux q component, λ lq is the leakage flux q component, λ ld is the leakage flux d component, L md，u is the unsaturated mutual inductance of the d-axis, L ld is the d-axis leakage inductance, Lmq, and u is the q-axis unsaturated mutual inductance, L lq is the q-axis leakage inductance, λ knee1，d is the d component of the first flux turning point, λ knee1，q is the q component of the first flux turning point; the parameter data of the synchronous machine include the saturation parameter b of the synchronous machine, the excitation current i f , direct axis damping winding D current i D , quadrature-axis damping winding g current i g and the quadrature-axis damping winding Q current i Q ,i r is the rotor current matrix, is the stator self-inductance dq0 matrix of the synchronous machine, R s is the stator resistance matrix of the synchronous machine, k is 2 / Δt, is the stator-rotor mutual inductance dq0 matrix of the synchronous machine, and are the stator current, stator voltage and stator flux phase domain matrices of the previous time step respectively.

[0057] Preferably, inverting the equivalent resistance matrix in the second Norton circuit to obtain an equivalent conductance matrix, and inputting the equivalent conductance matrix into the network conductance matrix for solving to obtain the three-phase voltage of the synchronous machine port comprises: inverting the equivalent resistance in the second Norton circuit to obtain the equivalent conductance matrix, and inputting the equivalent conductance matrix into the network conductance matrix before the time step cycle and solving it through the network solution equation to obtain the three-phase voltage of the synchronous machine port; wherein the network solution equation is YV=I, Y is the network conductance matrix, I is the current matrix composed of currents in the second Norton circuit, and V is the voltage matrix composed of the solved three-phase voltages of the synchronous machine port.

[0058] The present application also provides a synchronous machine electromagnetic transient modeling system based on saturation effect, comprising: a prediction data module, a first processing module, a first calculation and solution module, a second processing module, a second calculation and solution module and a judgment module;

[0059] The prediction data module is used to predict the first rotor angular velocity, the first rotor angle, the first current q component of the armature current, the first current d component of the armature current, the first flux turning point d component of the flux turning point and the first flux turning point q component of the flux turning point of the synchronous machine at a certain moment by using a linear extrapolation method, and to obtain the flux turning point data of the synchronous machine at a turning point between an unsaturated state and a saturated state;

[0060] The first processing module is used to process the first current q component and the first current d component to obtain a first Norton circuit simulating a synchronous machine; convert the first Norton circuit from a dq0 quantity to a second Norton circuit in abc phase form through coordinate transformation;

[0061] The first calculation and solution module is used to invert the equivalent resistance matrix in the second Norton circuit to obtain an equivalent conductance matrix, and input the equivalent conductance matrix into the network conductance matrix for solution to obtain the three-phase voltage of the synchronous machine port;

[0062] The second processing module processes the three-phase voltage to obtain a second current q component, a second current d component and a rotor current of the synchronous machine armature current; and processes the second current q component, the second current d component, the first flux turning point d component and the first flux turning point q component to obtain a stator flux d component, a stator flux q component and air gap flux data of the synchronous machine;

[0063] The second calculation and solution module is used to obtain the second rotor angular velocity and the second rotor angle of the synchronous machine by solving the second current q component, the second current d component, the stator flux d component and the stator flux q component in the mechanical system equation according to the air gap flux data being not greater than the saturated air gap flux threshold; and process the first flux turning point d component, the first flux turning point q component and the flux turning data to obtain the second flux turning point d component and the second flux turning point q component of the synchronous machine flux turning point;

[0064] The judgment module is used to calculate the second current q component, the second current d component, the second rotor angular velocity, the second rotor angle, the second flux turning point d component and the second flux turning point q component respectively with the corresponding first current q component, the first current d component, the first rotor angular velocity, the first rotor angle, the first flux turning point d component and the first flux turning point q component to obtain the corresponding absolute values ​​of the errors; if all the absolute values ​​of the errors are less than the allowable error value, the second rotor angular velocity and the second rotor angle of the synchronous machine are output.

[0065] The present application also provides a synchronous machine electromagnetic transient modeling device based on saturation effect, including a processor and a memory;

[0066] The memory is used to store program codes and transmit the program codes to the processor;

[0067] The processor is used to execute the synchronous machine electromagnetic transient modeling method based on saturation effect according to the instructions in the program code.

[0068] It can be seen from the above technical solutions that the embodiments of the present invention have the following advantages: The electromagnetic transient modeling method, system and device of the synchronous machine based on the saturation effect provided in the embodiments of the present application, the method comprises the following steps: S1. Using the linear extrapolation method to predict the first rotor angular velocity, the first rotor angle, the first current q component of the armature current, the first current d component of the armature current, the first flux turning point d component of the flux turning point and the first flux turning point q component of the flux turning point of the synchronous machine at a certain moment, and obtaining the flux turning point data of the synchronous machine at the turning point between the unsaturated state and the saturated state; S2. S1. Process the first current q component and the first current d component to obtain the first current d value, the second current q value and the third current 0 value of the equivalent resistance parallel current source; perform phase conversion on the first current d value, the second current value q and the third current 0 value to obtain the corresponding first current, second current and third current; S2. Input the first current, the second current and the third current into the network conductivity matrix for solution to obtain the first phase voltage, the second phase voltage and the third phase voltage corresponding to the first current, the second current and the third current; S3. Process the first phase voltage, the second phase voltage and the third phase voltage to obtain the synchronous The second current q component, the second current d component and the rotor current of the armature current of the synchronous machine are processed; and the second current q component, the second current d component, the first flux turning point d component and the first flux turning point q component are processed to obtain the stator flux d component, the stator flux q component and the air gap flux data of the synchronous machine; S5. If the air gap flux data is not greater than the saturated air gap flux threshold, the second current q component, the second current d component, the stator flux d component and the stator flux q component are solved in the mechanical system equation to obtain the second rotor angular velocity and the second rotor angle of the synchronous machine; the stator flux d component, the stator flux The flux q component and the flux turning point data are processed to obtain the second flux turning point d component and the second flux turning point q component of the synchronous machine flux turning point; S6. The second current q component, the second current d component, the second rotor angular velocity, the second rotor angle, the second flux turning point d component, and the second flux turning point q component are respectively subtracted from the corresponding first current q component, the first current d component, the first rotor angular velocity, the first rotor angle, the first flux turning point d component, and the first flux turning point q component to obtain the corresponding absolute value of the error; if all the absolute values ​​of the errors are less than the allowable error value, return to step S1.The synchronous machine electromagnetic transient modeling method based on saturation effect predicts the first rotor angular velocity, the first rotor angle, the first current q component, the first current d component, the first flux turning point d component, the first flux turning point q component and the flux turning data of the synchronous machine for analysis and processing, and obtains the corresponding second current q component, the second current d component, the second rotor angular velocity, the second rotor angle, the stator flux d component and the stator flux q component, which not only avoids the occurrence of historical quantities and current quantities of the synchronous machine's rotating potential, but also improves the accuracy of the simulation calculation results. The calculated results can maintain the calculation efficiency of the dq0 model on the basis of achieving the accuracy of the phase domain model, so that the synchronous machine electromagnetic transient modeling method based on saturation effect has high simulation accuracy and fast calculation efficiency. The synchronous machine electromagnetic transient modeling method based on saturation effect can be suitable for the development of electromagnetic transient simulation software for power system for actual engineering calculations, and solves the technical problems of low simulation accuracy and low efficiency in the existing electromagnetic transient simulation software for rotating motor modeling and simulation of saturation effect. BRIEF DESCRIPTION OF THE DRAWINGS

[0069] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative labor.

[0070] Figure 1 This is a flowchart of the steps of the synchronous machine electromagnetic transient modeling method based on the saturation effect according to an embodiment of the present application;

[0071] Figure 2 A magnetic flux line diagram of the saturated and unsaturated characteristic states of a synchronous machine according to a synchronous machine electromagnetic transient modeling method based on saturation effect described in an embodiment of the present application;

[0072] Figure 3 This is a framework diagram of a synchronous machine electromagnetic transient modeling system based on saturation effect according to an embodiment of the present application. DETAILED DESCRIPTION

[0073] In order to make the purpose, features and advantages of the present invention more obvious and easy to understand, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the embodiments described below are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0074] The embodiments of the present application provide a synchronous machine electromagnetic transient modeling method, system and device based on saturation effect, which are used to solve the technical problems of low simulation accuracy and low efficiency in the existing electromagnetic transient simulation software for rotating motor modeling and simulation with saturation effect.

[0075] Embodiment 1:

[0076] Figure 1 This is a flowchart of the steps of the synchronous machine electromagnetic transient modeling method based on saturation effect described in an embodiment of the present application. In the embodiment of the present application, a synchronous machine such as an engine is used as an example for explanation.

[0077] like Figure 1 As shown, the embodiment of the present application provides a synchronous machine electromagnetic transient modeling method based on saturation effect, comprising the following steps:

[0078] S1. Use linear extrapolation method to predict the first rotor angular velocity, the first rotor angle, the first current q component of the armature current, the first current d component of the armature current, the first flux turning point d component of the flux turning point and the first flux turning point q component of the flux turning point of the synchronous machine at a certain moment, and obtain the flux turning point data of the turning point between the unsaturated state and the saturated state of the synchronous machine.

[0079] In the embodiment of the present application, the linear extrapolation method is used to predict the first rotor angular velocity of the synchronous machine at a certain moment: ω1(t)=2ω1(t-Δt)-ω1(t-2Δt), t is the simulation time of the synchronous machine at a certain moment, and Δt is the simulation step size. Then, the first rotor angular velocity is processed by the trapezoidal integration method to obtain the first rotor angle.

[0080] It should be noted that the expression of the trapezoidal integration method is:

[0081]

[0082] In the embodiment of the present application, a linear extrapolation method is used to predict the first current q component and the first current d component of the armature current of the synchronous machine at a certain moment.

[0083] It should be noted that the expression for predicting the first current d component of the synchronous machine armature current at a certain moment by using the linear extrapolation method is: The expression for predicting the first current q component of the synchronous machine armature current at a certain moment by linear extrapolation method is:

[0084] In the embodiment of the present application, the linear extrapolation method is used to predict the first flux turning point d component and the first flux turning point q component of the synchronous machine flux turning point at a certain moment, and to obtain the flux turning point data λ of the synchronous machine at the turning point between the unsaturated state and the saturated state. knee .

[0085] It should be noted that the linear extrapolation method is used to predict the first flux turning point d component λ of the synchronous machine flux turning point at a certain moment. knee1，d and the first flux turning point q component λ knee1，q The expression is:

[0086] λ knee1，d (t) = 2λ knee1，d (t-Δt)-λ knee1，d (t-2Δt);

[0087] λ knee1，q (t) = 2λ knee1，q (t-Δt)-λ knee1，q (t-2Δt).

[0088] S2. Process the first current q component and the first current d component to obtain the first Norton circuit of the simulated synchronous machine; transform the first Norton circuit from dq0 quantity to the second Norton circuit in the form of abc phasor by coordinate transformation. It can be understood that the first current q component and the first current d component are processed to obtain the first current d value, the second current q value and the third current 0 value in the first Norton circuit; coordinate transformation is performed on the first current d value, the second current value q and the third current 0 value to obtain the first current, the second current and the third current of the abc phasor in the second Norton circuit.

[0089] It should be noted that the first current q component and the first current d component predicted in step S1 may be processed to construct an equivalent first Norton circuit of the synchronous machine, thereby obtaining a first current d value, a second current q value, and a third current 0 value of the first Norton circuit when the synchronous machine is in the case of an equivalent resistance parallel current source. Then, the first current d value, the second current q value, and the third current 0 value are converted from dq0 quantities into the first current, the second current, and the third current of the abc vector in the second Norton circuit.

[0090] S3. Invert the equivalent resistance matrix in the second Norton circuit to obtain an equivalent conductance matrix, and input the equivalent conductance matrix into the network conductance matrix for solution to obtain the three-phase voltage of the synchronous machine port. It can be understood that the first current, the second current and the third current are input into the network conductance matrix for solution to obtain the three-phase voltage of the synchronous machine port corresponding to the first current, the second current and the third current, and the three-phase voltages are a phase voltage, b phase voltage and c phase voltage respectively.

[0091] It should be noted that the process of inputting the first current, the second current and the third current into the network conductivity matrix for solving includes: before the time step cycle, inverting the equivalent resistance matrix in the second Norton circuit to obtain the equivalent conductance matrix, and then inputting it into the network conductivity matrix and solving it through the network solution equation to obtain the three-phase voltage of the synchronous machine port; wherein the network solution equation is YV=I, Y is the network conductivity matrix; I is the historical current source of the entire network, which includes the current matrix composed of the current in the second Norton circuit; V is the three-phase voltage value of the node to be solved in the entire network, which includes the voltage matrix composed of the three-phase voltage of the synchronous machine port to be solved.

[0092] S4. Process the three-phase voltage to obtain the second current q component, the second current d component and the rotor current of the synchronous machine's armature current, and process the second current q component, the second current d component, the first flux turning point d component and the first flux turning point q component to obtain the stator flux d component, the stator flux q component and the air gap flux data of the synchronous machine.

[0093] It should be noted that the second current q component, the second current d component and the rotor current of the synchronous machine's armature current are calculated based on the three data of phase a voltage, phase b voltage and phase c voltage of the three-phase voltage after Park transformation, and the stator flux d component, stator flux q component and air gap flux data of the synchronous machine are calculated based on the second current q component, the second current d component, the first flux turning point d component and the first flux turning point q component and the parameters in the inductance parameter matrix of the synchronous machine.

[0094] S5. If the air gap flux data is not greater than the saturated air gap flux threshold, the second current q component, the second current d component, the stator flux d component and the stator flux q component are solved in the mechanical system equation to obtain the second rotor angular velocity and the second rotor angle of the synchronous machine; the first flux turning point d component, the first flux turning point q component and the flux turning data are also processed to obtain the second flux turning point d component and the second flux turning point q component of the synchronous machine flux turning point.

[0095] It should be noted that the second current q component, the second current d component, the stator flux d component, the stator flux q component and the rotor current data obtained in step S4 are mainly input into the mechanical system equation to calculate the second rotor angular velocity and the second rotor angle of the synchronous machine. The first flux turning point d component, the first flux turning point q component and the flux turning data are processed to obtain the second flux turning point d component and the second flux turning point q component of the synchronous machine flux turning point to provide data for step S6 analysis. In this embodiment, the saturated air gap flux threshold can be set according to demand and is not limited in detail here.

[0096] In the embodiment of the present application, the mechanical system equation is:

[0097]

[0098]

[0099]

[0100] Where p is the number of pole pairs of the synchronous machine, λ q is the stator flux q component, λ d is the stator flux d component, is the second current d component, is the second current q component, J is the moment of inertia of the synchronous machine, D is the viscosity and air friction damping coefficient of the synchronous machine, T is the constant of the synchronous machine, ω is the second rotor angular velocity, θ is the second rotor angle, and t is the simulation time. Among them, T is the power P0 of the synchronous machine and the initial angular velocity ω of the synchronous machine s The ratio of, that is, T = P0 / ω s .

[0101] S6. Calculate the second current q component, the second current d component, the second rotor angular velocity, the second rotor angle, the second magnetic flux turning point d component, and the second magnetic flux turning point q component respectively with the corresponding first current q component, the first current d component, the first rotor angular velocity, the first rotor angle, the first magnetic flux turning point d component, and the first magnetic flux turning point q component to obtain the corresponding absolute values ​​of the errors; if all the absolute values ​​of the errors are less than the allowable error value, return to step S1.

[0102] The main step is to perform difference processing on the second current q component, the second current d component, the second rotor angular velocity, the second rotor angle, the second flux turning point d component, and the second flux turning point q component calculated in steps S4 and S5 respectively with the first current q component, the first current d component, the first rotor angular velocity, the first rotor angle, the first flux turning point d component, and the first flux turning point q component predicted in step S1 to obtain the corresponding absolute value of the error, and then determine whether all the absolute values ​​of the errors are less than the allowable error value. If so, return to step S1 to perform the next synchronous machine electromagnetic transient modeling. If not, the absolute value of any one of the errors is not less than the allowable error value, return to step S4, and recalculate the second current q component, the second current d component, the second rotor angular velocity, the second rotor angle, the second flux turning point d component, and the second flux turning point q component.

[0103] The present application provides a synchronous machine electromagnetic transient modeling method based on saturation effect, comprising the following steps: S1. using a linear extrapolation method to predict the first rotor angular velocity, the first rotor angle, the first current q component of the armature current, the first current d component of the armature current, the first flux turning point d component of the flux turning point and the first flux turning point q component of the flux turning point of the synchronous machine at a certain moment, and obtaining the flux turning point data of the synchronous machine at the turning point between the unsaturated state and the saturated state; S2. processing the first current q component and the first current d component to obtain a first Norton circuit simulating the synchronous machine; converting the first Norton circuit from a dq0 quantity to a second Norton circuit in the form of abc phasor through coordinate transformation; S3. inverting the equivalent resistance matrix in the second Norton circuit to obtain an equivalent conductance matrix, and inputting the equivalent conductance matrix into the network conductance matrix for solution to obtain the three-phase voltage of the synchronous machine port; S4. processing the three-phase voltage to obtain the second current q component, the second current d component and the rotor current of the armature current of the synchronous machine; and The d component, the first flux turning point d component and the first flux turning point q component are processed to obtain the stator flux d component, the stator flux q component and the air gap flux data of the synchronous machine; S5. If the air gap flux data is not greater than the saturated air gap flux threshold, the second current q component, the second current d component, the stator flux d component and the stator flux q component are solved in the mechanical system equation to obtain the second rotor angular velocity and the second rotor angle of the synchronous machine; the stator flux d component, the stator flux q component and the flux turning point data are also processed to obtain the synchronous machine. The second flux turning point d component and the second flux turning point q component of the machine flux turning point; S6. Calculate the second current q component, the second current d component, the second rotor angular velocity, the second rotor angle, the second flux turning point d component, and the second flux turning point q component with the corresponding first current q component, the first current d component, the first rotor angular velocity, the first rotor angle, the first flux turning point d component, and the first flux turning point q component to obtain the corresponding absolute values ​​of the errors; if all the absolute values ​​of the errors are less than the allowable error value, return to step S1.The synchronous machine electromagnetic transient modeling method based on saturation effect predicts the first rotor angular velocity, the first rotor angle, the first current q component, the first current d component, the first flux turning point d component, the first flux turning point q component and the flux turning data of the synchronous machine for analysis and processing, and obtains the corresponding second current q component, the second current d component, the second rotor angular velocity, the second rotor angle, the stator flux d component and the stator flux q component, which not only avoids the occurrence of historical quantities and current quantities of the synchronous machine's rotating potential, but also improves the accuracy of the simulation calculation results. The calculated results can maintain the calculation efficiency of the dq0 model on the basis of achieving the accuracy of the phase domain model, so that the synchronous machine electromagnetic transient modeling method based on saturation effect has high simulation accuracy and fast calculation efficiency. The synchronous machine electromagnetic transient modeling method based on saturation effect can be suitable for the development of electromagnetic transient simulation software for power system for actual engineering calculations, and solves the technical problems of low simulation accuracy and low efficiency in the existing electromagnetic transient simulation software for rotating motor modeling and simulation of saturation effect.

[0104] In one embodiment of the present application, the first current q component and the first current d component are processed to obtain a first Norton circuit simulating a synchronous machine; and the first Norton circuit is converted from a dq0 quantity to a second Norton circuit in abc phasor form by coordinate transformation, including:

[0105] Obtaining the stator and rotor voltage equations of the synchronous machine, and performing discretization processing using an implicit trapezoidal integration method according to the stator and rotor voltage equations to obtain a first transformation equation;

[0106] The first transformation equation is subjected to Park transformation, rotor variables are eliminated, and the dq axis is processed using average resistance to obtain the Thevenin equation on the stator side;

[0107] The Thevenin equation on the stator side is converted into the first Norton circuit simulating the synchronous machine through mathematical transformation;

[0108] The first Norton circuit is converted from the dq0 quantity to the second Norton circuit of the abc phase using the phase coordinate transformation formula;

[0109] Among them, the first Norton circuit equation is:

[0110]

[0111] R ave =(R d +R q ) / 2;

[0112] The phase conversion formula is:

[0113]

[0114] In the formula, is the first current d component, is the first current q component, R d , R q R0 and R1 are the resistance parameters of the resistance matrix in the Thevenin equation on the stator side, e d 、e q , e0 are the voltage parameters of the voltage source matrix in the Thevenin equation on the stator side, i d，source is the first current d value, i q，source is the second current q value, i 0，source is the third current 0 value, θ1 is the first rotor angle, i a，source is the first current, i b，source is the second current, i c，source is the third current.

[0115] In an embodiment of the present application, the stator and rotor voltage equations and the flux equations of the synchronous machine are obtained, and the stator and rotor voltage equations are discretized using the implicit trapezoidal integration method to obtain a first transformation equation; the first transformation equation is subjected to a Park transformation and the rotor variables are eliminated to obtain the Thevenin equation on the stator side.

[0116] It should be noted that the stator and rotor voltage equations are:

[0117]

[0118] The magnetic flux equation is:

[0119]

[0120]

[0121]

[0122] The first transformation equation is:

[0123]

[0124]

[0125]

[0126]

[0127] The Thevenin equation on the stator side is:

[0128]

[0129]

[0130]

[0131] In the formula, are the stator voltage, current and flux parameters of the synchronous machine phase domain matrix, v r 、i r , r are the stator voltage, current, and flux parameters of the synchronous machine flux matrix, R s , R r are the stator resistance and rotor resistance of the synchronous machine, L(θ1) is the inductance related to the first rotor angle in the synchronous machine, and L ss , L rr are the stator inductance and rotor inductance in the self-inductance matrix of the synchronous machine, L sr , L rs are the stator inductance and rotor inductance in the mutual inductance matrix of the synchronous machine, k is 2 / Δt, and the variables with ^ are the values ​​of the variables in the previous time step, that is, the historical quantities, R dq0 、e dq0 They are the resistance matrix and the series voltage source matrix in the Thevenin equation on the stator side, are the phase domain value of the synchronous machine stator saturation flux turning point and its corresponding dq0 component, λ rk is the phase domain value of the synchronous machine rotor saturation flux turning point. middle Directly obtained from the historical quantities of the network solution, and The historical variables of current and flux dq0 are obtained by Park transformation, v r Use the value from the previous moment.

[0132] In an embodiment of the present application, a Park transformation is performed on the first transformation equation to obtain a second transformation equation, and then the rotor variable in the second transformation equation is eliminated to obtain the Thevenin equation on the stator side.

[0133] It should be noted that the second transformation equation is:

[0134]

[0135]

[0136] In the formula, are the stator voltage dq0 matrix and stator current dq0 matrix of the synchronous machine respectively, is the stator self-inductance dq0 matrix of the synchronous machine, They are the stator-rotor mutual inductance dq0 matrix and the rotor mutual inductance dq0 matrix of the synchronous machine respectively.

[0137] Figure 2 This is a magnetic flux line diagram of the saturated and unsaturated characteristic states of a synchronous machine according to the electromagnetic transient modeling method of a synchronous machine based on the saturation effect described in an embodiment of the present application.

[0138] In one embodiment of the present application, the inductance matrix L of the synchronous machine dq0 The expression is:

[0139]

[0140] Where b is the saturation parameter of the synchronous machine, L md，u is the unsaturated mutual inductance of the d-axis, L ld is the d-axis leakage inductance, L mq，u is the unsaturated mutual inductance of the q-axis, L lq is the q-axis leakage inductance, L lQ is the leakage inductance of the synchronous machine quadrature-axis damping winding Q, L lg is the leakage inductance of the synchronous machine quadrature-axis damping winding g, L lf is the leakage inductance of the synchronous machine excitation winding.

[0141] It should be noted that if Figure 2 As shown, the inductance matrix L of the synchronous machine dq0 In the expression of, b is M when the synchronous machine is saturated. slope / M du , b is 1 when unsaturated, and the no-load saturation curve of the synchronous machine is approximated by a broken line with two slopes. The slopes of unsaturated and saturated are M du and M slope Here, the subscript m represents the mutual inductance value, u represents the unsaturated value, and d, q, and 0 represent the dq0 axis components, respectively.

[0142] In the embodiment of the present application, the above-mentioned Thevenin equation on the stator side is obtained under a dq0 model. In order to avoid generating a time-varying asymmetric 3×3 resistive matrix and improve the calculation accuracy, the average resistance is used in the dq axis to obtain the modified Thevenin equation converted to the stator side. The converted second Norton circuit is:

[0143]

[0144] In one embodiment of the present application, the three-phase voltage is processed to obtain the second current q component, the second current d component and the rotor current of the armature current of the synchronous machine; and the second current q component, the second current d component, the first flux turning point d component and the first flux turning point q component are processed to obtain the stator flux d component, the stator flux q component and the air gap flux data of the synchronous machine, including:

[0145] The three-phase voltage is converted by Park transformation to obtain the dq0 axis voltage component corresponding to the three-phase voltage;

[0146] According to the matrix parameters of the Thevenin equation on the stator side and the dq0 axis voltage component, the second current q component and the second current d component of the armature current of the synchronous machine are calculated by the armature current calculation formula;

[0147] Based on the parameter data of the synchronous machine, the dq0 axis voltage component, the second current q component and the second current d component, the rotor current of the synchronous machine is calculated by a rotor current calculation formula;

[0148] Based on the parameter data of the synchronous machine and the second current q component, the second current d component, the first flux turning point d component and the first flux turning point q component, the stator flux d component, the stator flux q component and the air gap flux data of the synchronous machine are calculated by the stator flux dq component calculation formula;

[0149] The Parker transform is:

[0150]

[0151] The armature current calculation formula is:

[0152]

[0153] R ave =(R d +R q ) / 2;

[0154] The rotor current calculation formula is:

[0155]

[0156]

[0157]

[0158]

[0159]

[0160]

[0161]

[0162] The calculation formula for the stator flux dq component is:

[0163] λ d =λ ld +λ md ; q =λ lq +λ mq ;

[0164]

[0165]

[0166] In the formula, is the second current d component, is the second current q component, R d , R q R0 and R1 are the resistance parameters of the resistance matrix in the Thevenin equation on the stator side, e d 、e q , e0 are the voltage parameters of the voltage source matrix in the Thevenin equation on the stator side, θ1 is the first rotor angle, v a is the phase a voltage of the three-phase voltage, v b is the b-phase voltage of the three-phase voltage, v c is the c-phase voltage of the three-phase voltage, v d is the first voltage d component of the dq0 axis voltage component, v q is the second voltage q component of the dq0 axis voltage component, v0 is the third voltage 0 component of the dq0 axis voltage component, λ d is the stator flux d component, λ q is the stator flux q component, λ m is the air gap flux data, λ md is the air gap flux d component, λ mq is the air gap flux q component, λ lq is the leakage flux q component, λ ld is the leakage flux d component, L md，u is the unsaturated mutual inductance of the d-axis, L ld is the d-axis leakage inductance, L mq，u is the unsaturated mutual inductance of the q-axis, L lq is the q-axis leakage inductance, λ knee1，d is the d component of the first flux turning point, λ knee1，q is the q component of the first flux turning point; the parameter data of the synchronous machine include the saturation parameter b of the synchronous machine, the excitation current i f , direct-axis damping winding D current i D , quadrature-axis damping winding g current i g and the quadrature-axis damping winding Q current i Q ,i r is the rotor current matrix, is the stator self-inductance dq0 matrix of the synchronous machine, R s is the stator resistance matrix of the synchronous machine, k is 2 / Δt, is the stator-rotor mutual inductance dq0 matrix of the synchronous machine, and are the stator current, stator voltage and stator flux phase domain matrices of the previous time step respectively.

[0167] In the embodiment of the present application, in step S3, an equivalent resistance matrix R is constructed according to the Thevenin equation on the stator side. equiv The synchronous machine electromagnetic transient modeling method based on saturation effect obtains the equivalent conductance matrix by inverting the equivalent resistance matrix, and inputs the first current, the second current and the third current and the equivalent conductance matrix into the network equation at one time before the simulation time step cycle and solves them through the network solution equation to obtain the a-phase voltage, the b-phase voltage and the c-phase voltage.

[0168] It should be noted that the equivalent resistance matrix is:

[0169]

[0170] In one embodiment of the present application, the first flux turning point d component, the first flux turning point q component and the flux turning data are processed to obtain the second flux turning point d component and the second flux turning point q component of the synchronous machine flux turning point, including:

[0171] Based on the parameter data of the synchronous machine, the first flux turning point d component, the first flux turning point q component, the second current q component and the second current d component, the air gap flux d component and the air gap flux q component are calculated, and the ratio of the air gap flux q component to the air gap flux d component is calculated by using an inverse trigonometric function to obtain the flux rotation angle of the synchronous machine flux turning point;

[0172] The flux turning point data and the flux turning angle are calculated using trigonometric functions to obtain the second flux turning point d component and the second flux turning point q component of the synchronous machine flux turning point.

[0173] It should be noted that the calculation of the air gap flux d component λ md and the air gap flux q component λ mq The expressions are:

[0174]

[0175] The inverse trigonometric function expression for calculating the magnetic flux rotation angle β is: β = tan -1 (λ mq / λ md ) Calculate the second flux turning point d component λ of the synchronous machine flux turning point knee2，d and the second flux turning point q component λ knee2，q The expressions are:

[0176] λ knee2，d =λ knee cosβ,λ knee2，q =λ kneesinβ,λ knee It is the magnetic flux turning data.

[0177] in, λ rk =[λ knee，d λ knee，d λ knee，q λ knee，q ] T .

[0178] In an embodiment of the present application, the synchronous machine electromagnetic transient modeling method based on the saturation effect includes: if the air gap flux data is greater than the saturation air gap flux threshold, the air gap flux d component, the air gap flux q component and the air gap flux data are corrected using a saturation correction parameter to obtain a modified stator flux d component and a stator flux q component.

[0179] It should be noted that when the air gap flux data is greater than the saturation air gap flux threshold, it means that the synchronous machine has entered the saturation state. At this time, the saturation parameter b of the synchronous machine is the value when the synchronous machine is saturated, that is, Figure 2 b=M slope / M du The inductance matrix L of the synchronous machine is modified according to the saturation parameter b dq0 , and obtain the corrected L md，u and L mq，u The values ​​of these two parameters are corrected using the air gap flux d component λ md and the air gap flux q component λ mq The corrected air gap flux d component λ is calculated from the expression md and the air gap flux q component λ mq , thereby obtaining the corrected air gap flux data and realizing the correction of the air gap flux data.

[0180] In the embodiment of the present application, the synchronous machine electromagnetic transient modeling method based on the saturation effect only needs to predict the armature current and the dq component of the flux turning point to establish the synchronous machine equivalent Norton circuit through step S2, and does not need to predict the rotating potential, avoiding the cumulative error caused by the predicted value, thereby improving the simulation accuracy. Through steps S2 and S3, the phase domain discrete model of the synchronous machine can be transformed by Parker, ensuring that the equivalent resistance matrix in the synchronous machine equivalent Norton circuit is a constant, and only when the synchronous machine enters and exits the saturation state, the equivalent resistance matrix is ​​changed and the triangular decomposition is re-performed to obtain the air gap flux data, thereby maintaining the computational efficiency of the synchronous machine dq0 model. Through step S5, the next simulation is performed only when the air gap flux data is not greater than the saturated air gap flux threshold. Compared with the synchronous machine electromagnetic transient modeling without considering saturation, the synchronous machine electromagnetic transient modeling method based on the saturation effect considers the situation of the saturated synchronous machine, and improves the electromagnetic steady-state and transient modeling accuracy of the synchronous machine.

[0181] Embodiment 2:

[0182] Figure 3 This is a framework diagram of a synchronous machine electromagnetic transient modeling system based on saturation effect according to an embodiment of the present application.

[0183] like Figure 3 As shown, the present application also provides a synchronous machine electromagnetic transient modeling system based on saturation effect, including a prediction data module 10, a first processing module 20, a first calculation and solution module 30, a second processing module 40, a second calculation and solution module 50 and a judgment module 60;

[0184] The prediction data module 10 is used to predict the first rotor angular velocity, the first rotor angle, the first current q component of the armature current, the first current d component of the armature current, the first flux turning point d component of the flux turning point and the first flux turning point q component of the flux turning point of the synchronous machine at a certain moment by using a linear extrapolation method, and to obtain the flux turning point data of the synchronous machine at the turning point between the unsaturated state and the saturated state;

[0185] The first processing module 20 is used to process the first current q component and the first current d component to obtain a first Norton circuit simulating a synchronous machine; and to convert the first Norton circuit from a dq0 quantity into a second Norton circuit in the form of an abc phase quantity through coordinate transformation;

[0186] The first calculation and solution module 30 is used to invert the equivalent resistance matrix in the second Norton circuit to obtain an equivalent conductance matrix, and input the equivalent conductance matrix into the network conductance matrix for solution to obtain the three-phase voltage of the synchronous machine port;

[0187] The second processing module 40 is used to process the three-phase voltage to obtain the second current q component, the second current d component and the rotor current of the synchronous machine armature current; and to process the second current q component, the second current d component, the first flux turning point d component and the first flux turning point q component to obtain the stator flux d component, the stator flux q component and the air gap flux data of the synchronous machine;

[0188] The second calculation and solution module 50 is used to obtain the second rotor angular velocity and the second rotor angle of the synchronous machine by solving the second current q component, the second current d component, the stator flux d component and the stator flux q component in the mechanical system equation according to the air gap flux data being not greater than the saturated air gap flux threshold; and to process the first flux turning point d component, the first flux turning point q component and the flux turning data to obtain the second flux turning point d component and the second flux turning point q component of the synchronous machine flux turning point;

[0189] The judgment module 60 is used to calculate the second current q component, the second current d component, the second rotor angular velocity, the second rotor angle, the second flux turning point d component and the second flux turning point q component with the corresponding first current q component, the first current d component, the first rotor angular velocity, the first rotor angle, the first flux turning point d component and the first flux turning point q component, respectively, to obtain the corresponding absolute values ​​of the errors; if all the absolute values ​​of the errors are less than the allowable error value, the second rotor angular velocity and the second rotor angle of the synchronous machine are output.

[0190] In the embodiment of the present application, the mechanical system equation is:

[0191]

[0192]

[0193]

[0194] Where p is the number of pole pairs of the synchronous machine, λ q is the stator flux q component, λ d is the stator flux d component, is the second current d component, is the second current q component, J is the moment of inertia of the synchronous machine, D is the viscosity and air friction damping coefficient of the synchronous machine, T is the constant of the synchronous machine, ω is the second rotor angular velocity, θ is the second rotor angle, and t is the simulation time.

[0195] It should be noted that the contents of the modules in the system of the second embodiment have been described in detail in the contents of the steps in the method of the first embodiment, and the contents of the modules in the system of the second embodiment will not be described in detail here.

[0196] Embodiment three:

[0197] The present application also provides a synchronous machine electromagnetic transient modeling device based on saturation effect, including a processor and a memory;

[0198] A memory, used for storing program codes and transmitting the program codes to a processor;

[0199] The processor is used to execute the above-mentioned synchronous machine electromagnetic transient modeling method based on saturation effect according to the instructions in the program code.

[0200] It should be noted that the synchronous machine electromagnetic transient modeling method based on saturation effect has been described in detail in the first embodiment and will not be elaborated here. The processor is used to execute the steps in the embodiment of the synchronous machine electromagnetic transient modeling method based on saturation effect according to the instructions in the program code. Alternatively, the processor implements the functions of each module / unit in the above-mentioned system / device embodiments when executing the computer program.

[0201] Exemplarily, the computer program may be divided into one or more modules / units, one or more modules / units are stored in a memory and executed by a processor to complete the present application. One or more modules / units may be a series of computer program instruction segments capable of completing a specific function, and the instruction segments are used to describe the execution process of the computer program in a terminal device.

[0202] The terminal device may be a computing device such as a desktop computer, a notebook, a PDA, or a cloud server. The terminal device may include, but is not limited to, a processor and a memory. Those skilled in the art will appreciate that this does not constitute a limitation on the terminal device, and may include more or fewer components than shown in the figure, or a combination of certain components, or different components. For example, the terminal device may also include input and output devices, network access devices, buses, etc.

[0203] The processor may be a central processing unit (CPU), other general-purpose processors, digital signal processors (DSP), application-specific integrated circuits (ASIC), field-programmable gate arrays (FPGA) or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. A general-purpose processor may be a microprocessor or any conventional processor, etc.

[0204] The memory may be an internal storage unit of the terminal device, such as a hard disk or memory of the terminal device. The memory may also be an external storage device of the terminal device, such as a plug-in hard disk, a smart memory card (SmartMedia Card, SMC), a secure digital (Secure Digital, SD) card, a flash card (Flash Card), etc. equipped on the terminal device. Furthermore, the memory may also include both an internal storage unit of the terminal device and an external storage device. The memory is used to store computer programs and other programs and data required by the terminal device. The memory may also be used to temporarily store data that has been output or is to be output.

[0205] Those skilled in the art can clearly understand that, for the convenience and brevity of description, the specific working processes of the systems, devices and units described above can refer to the corresponding processes in the aforementioned method embodiments and will not be repeated here.

[0206] In the several embodiments provided in the present application, it should be understood that the disclosed systems, devices and methods can be implemented in other ways. For example, the device embodiments described above are only schematic. For example, the division of the units is only a logical function division. There may be other division methods in actual implementation, such as multiple units or components can be combined or integrated into another system, or some features can be ignored or not executed. Another point is that the mutual coupling or direct coupling or communication connection shown or discussed can be an indirect coupling or communication connection through some interfaces, devices or units, which can be electrical, mechanical or other forms.

[0207] The units described as separate components may or may not be physically separated, and the components shown as units may or may not be physical units, that is, they may be located in one place or distributed on multiple network units. Some or all of the units may be selected according to actual needs to achieve the purpose of the solution of this embodiment.

[0208] In addition, each functional unit in each embodiment of the present invention may be integrated into one processing unit, or each unit may exist physically separately, or two or more units may be integrated into one unit. The above-mentioned integrated unit may be implemented in the form of hardware or in the form of software functional units.

[0209] If the integrated unit is implemented in the form of a software functional unit and sold or used as an independent product, it can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present invention is essentially or the part that contributes to the prior art or all or part of the technical solution can be embodied in the form of a software product. The computer software product is stored in a storage medium, including several instructions for a computer device (which can be a personal computer, a server, or a network device, etc.) to perform all or part of the steps of the method described in each embodiment of the present invention. The aforementioned storage medium includes: U disk, mobile hard disk, read-only memory (ROM, Read-Only Memory), random access memory (RAM, Random Access Memory), disk or optical disk and other media that can store program codes.

[0210] As described above, the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit the same. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that the technical solutions described in the aforementioned embodiments may still be modified, or some of the technical features thereof may be replaced by equivalents. However, these modifications or replacements do not deviate the essence of the corresponding technical solutions from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A synchronous machine electromagnetic transient modeling method based on saturation effect, characterized in that: The following steps are involved: S1. Use linear extrapolation method to predict the first rotor angular velocity, the first rotor angle, the first current q component of the armature current, the first current d component of the armature current, the first flux turning point d component of the flux turning point and the first flux turning point q component of the flux turning point of the synchronous machine at a certain moment, and obtain the flux turning point data of the synchronous machine at the turning point between the unsaturated state and the saturated state; S2. Processing the first current q component and the first current d component to obtain a first Norton circuit simulating a synchronous machine; converting the first Norton circuit from a dq0 quantity to a second Norton circuit in abc phasor form by coordinate transformation; S3. Inverting the equivalent resistance matrix in the second Norton circuit to obtain an equivalent conductance matrix, and inputting the equivalent conductance matrix into the network conductance matrix for solution to obtain the three-phase voltage of the synchronous machine port; S4. Processing the three-phase voltage to obtain the second current q component, the second current d component and the rotor current of the synchronous machine armature current; and processing the second current q component, the second current d component, the first flux turning point d component and the first flux turning point q component to obtain the stator flux d component, the stator flux q component and the air gap flux data of the synchronous machine; S5. If the air gap flux data is not greater than the saturated air gap flux threshold, the second current q component, the second current d component, the stator flux d component and the stator flux q component are solved in the mechanical system equation to obtain the second rotor angular velocity and the second rotor angle of the synchronous machine; the stator flux d component, the stator flux q component and the flux turning data are also processed to obtain the second flux turning point d component and the second flux turning point q component of the synchronous machine flux turning point; S6. Calculate the second current q component, the second current d component, the second rotor angular velocity, the second rotor angle, the second flux turning point d component and the second flux turning point q component respectively with the corresponding first current q component, the first current d component, the first rotor angular velocity, the first rotor angle, the first flux turning point d component and the first flux turning point q component to obtain corresponding absolute values ​​of errors; if all absolute values ​​of errors are less than the allowable error value, return to step S1; Processing the three-phase voltage to obtain the second current q component, the second current d component and the rotor current of the synchronous machine armature current includes: The three-phase voltage is converted by Park transformation to obtain a dq0-axis voltage component corresponding to the three-phase voltage; According to the matrix parameters of the Thevenin equation on the stator side and the dq0-axis voltage component, the second current q component and the second current d component of the armature current of the synchronous machine are calculated by the armature current calculation formula; Based on the parameter data of the synchronous machine, the dq0-axis voltage component, the second current q component and the second current d component, the rotor current of the synchronous machine is calculated by a rotor current calculation formula; Wherein, the Parker transform is: ； The armature current calculation formula is: ； The rotor current calculation formula is: ； ； ； ； In the formula, is the first current d component, is the first current q component, is the second current d component, is the second current q component, R d , R q R0 and R1 are the resistance parameters of the resistance matrix in the Thevenin equation on the stator side, e d 、e q , e0 are the voltage parameters of the voltage source matrix in the Thevenin equation on the stator side, is the first rotor angle, v a is the phase a voltage of the three-phase voltage, v b is the b-phase voltage of the three-phase voltage, v c is the c-phase voltage of the three-phase voltage, v d is the first voltage d component of the dq0 axis voltage component, v q is the second voltage q component of the dq0 axis voltage component, v0 is the third voltage 0 component of the dq0 axis voltage component, and the parameter data of the synchronous machine include the saturation parameter b of the synchronous machine, the excitation current i f , direct axis damping winding D current i D , quadrature-axis damping winding g current i g and the quadrature-axis damping winding Q current i Q ,i r is the rotor current matrix, is the stator self-inductance dq0 matrix of the synchronous machine, R s is the stator resistance matrix of the synchronous machine, k is 2 / , is the stator-rotor mutual inductance dq0 matrix of the synchronous machine, , and are the stator current, stator voltage and stator flux phase domain matrices of the previous time step respectively.

2. The synchronous machine electromagnetic transient modeling method based on saturation effect according to claim 1 is characterized in that: include: If the absolute value of any one of the errors is not less than the allowable error value, return to step S4.

3. The synchronous machine electromagnetic transient modeling method based on saturation effect according to claim 1 is characterized in that: Processing the first flux turning point d component, the first flux turning point q component and the flux turning data to obtain a second flux turning point d component and a second flux turning point q component of the synchronous machine flux turning point includes: Based on the parameter data of the synchronous machine, the first flux turning point d component, the first flux turning point q component, the second current q component and the second current d component, an air gap flux d component and an air gap flux q component are calculated, and an inverse trigonometric function is used to calculate the ratio of the air gap flux q component to the air gap flux d component to obtain a flux rotation angle at the flux turning point of the synchronous machine; The flux turning point data and the flux turning angle are calculated using trigonometric functions to obtain a second flux turning point d component and a second flux turning point q component of the synchronous machine flux turning point.

4. The synchronous machine electromagnetic transient modeling method based on saturation effect according to claim 3 is characterized in that: include: If the air gap flux data is greater than a saturated air gap flux threshold, a saturation correction parameter is used to correct the air gap flux d component, the air gap flux q component and the air gap flux data to obtain modified stator flux d component and stator flux q component.

5. The synchronous machine electromagnetic transient modeling method based on saturation effect according to claim 1 is characterized in that: The mechanical system equation is: ； ； Where p is the number of pole pairs of the synchronous machine, is the stator flux q component, is the stator flux d component, J is the moment of inertia of the synchronous machine, D is the viscosity and air friction damping coefficient of the synchronous machine, T is the constant of the synchronous machine, is the second rotor angular velocity, is the second rotor angle, and t is the simulation time.

6. The synchronous machine electromagnetic transient modeling method based on saturation effect according to claim 1 is characterized in that: Processing the first current q component and the first current d component to obtain a first Norton circuit simulating a synchronous machine; The second Norton circuit that converts the first Norton circuit from the dq0 quantity into the abc phase quantity form by coordinate transformation includes: Obtaining the stator and rotor voltage equations of the synchronous machine, and performing discretization processing using an implicit trapezoidal integration method according to the stator and rotor voltage equations to obtain a first transformation equation; Performing Park transformation on the first transformation equation, eliminating rotor variables, and using average resistance processing on the dq axis to obtain the Thevenin equation on the stator side; The Thévenin equation on the stator side is converted into a first Norton circuit simulating a synchronous machine by mathematical transformation; The first Norton circuit is converted from a dq0 quantity to a second Norton circuit of abc phase using a phase coordinate transformation formula; Wherein, the first Norton circuit is: ； The phase coordinate transformation formula is: ； In the formula, is the first current d value of the first Norton circuit, is the second current q value of the first Norton circuit, is the third current 0 value of the first Norton circuit, is the first rotor angle, is the first current of the a-phase current source of the second Norton circuit, is the second current of the second Norton circuit b-phase current source, is the third current of the second Norton circuit c-phase current source.

7. The synchronous machine electromagnetic transient modeling method based on saturation effect according to claim 1 is characterized in that: Processing the second current q component, the second current d component, the first flux turning point d component, and the first flux turning point q component to obtain the stator flux d component, the stator flux q component, and the air gap flux data of the synchronous machine includes: Based on the parameter data of the synchronous machine and the second current q component, the second current d component, the first flux turning point d component and the first flux turning point q component, the stator flux d component, the stator flux q component and the air gap flux data of the synchronous machine are calculated by a stator flux dq component calculation formula; The calculation formula of the stator flux dq component is: ； In the formula, is the stator flux d component, is the stator flux q component, is the air gap flux data, is the air gap flux d component, is the air gap flux q component, is the leakage flux q component, is the leakage flux d component, L md，u is the unsaturated mutual inductance of the d-axis, L ld is the d-axis leakage inductance, L mq，u is the unsaturated mutual inductance of the q-axis, L lq is the q-axis leakage inductance, is the d component of the first flux turning point, is the q component of the first flux turning point.

8. The synchronous machine electromagnetic transient modeling method based on saturation effect according to claim 1 is characterized in that: Inverting the equivalent resistance matrix in the second Norton circuit to obtain an equivalent conductance matrix, and inputting the equivalent conductance matrix into the network conductance matrix for solving to obtain the three-phase voltage at the synchronous machine port includes: inverting the equivalent resistance in the second Norton circuit to obtain the equivalent conductance matrix, and inputting the equivalent conductance matrix into the network conductance matrix before the time step cycle and solving it through the network solution equation to obtain the three-phase voltage at the synchronous machine port; wherein the network solution equation is YV=I, Y is the network conductance matrix, I is the current matrix composed of currents in the second Norton circuit, and V is the voltage matrix composed of the solved three-phase voltages at the synchronous machine port.

9. A synchronous machine electromagnetic transient modeling system based on saturation effect, characterized in that: include: A prediction data module, a first processing module, a first calculation and solution module, a second processing module, a second calculation and solution module and a judgment module; The prediction data module is used to predict the first rotor angular velocity, the first rotor angle, the first current q component of the armature current, the first current d component of the armature current, the first flux turning point d component of the flux turning point and the first flux turning point q component of the flux turning point of the synchronous machine at a certain moment by using a linear extrapolation method, and to obtain the flux turning point data of the synchronous machine at a turning point between an unsaturated state and a saturated state; The first processing module is used to process the first current q component and the first current d component to obtain a first Norton circuit simulating a synchronous machine; Converting the first Norton circuit from a dq0 quantity to a second Norton circuit in abc phase form by coordinate transformation; The first calculation and solution module is used to invert the equivalent resistance matrix in the second Norton circuit to obtain an equivalent conductance matrix, and input the equivalent conductance matrix into the network conductance matrix for solution to obtain the three-phase voltage of the synchronous machine port; The second processing module is used to process the three-phase voltage to obtain the second current q component, the second current d component and the rotor current of the synchronous machine armature current; and process the second current q component, the second current d component, the first flux turning point d component and the first flux turning point q component to obtain the stator flux d component, the stator flux q component and the air gap flux data of the synchronous machine; The second calculation and solution module is used to obtain the second rotor angular velocity and the second rotor angle of the synchronous machine by solving the second current q component, the second current d component, the stator flux d component and the stator flux q component in the mechanical system equation according to the air gap flux data being not greater than the saturated air gap flux threshold; and process the first flux turning point d component, the first flux turning point q component and the flux turning data to obtain the second flux turning point d component and the second flux turning point q component of the synchronous machine flux turning point; The judgment module is used to calculate the second current q component, the second current d component, the second rotor angular velocity, the second rotor angle, the second flux turning point d component and the second flux turning point q component respectively with the corresponding first current q component, the first current d component, the first rotor angular velocity, the first rotor angle, the first flux turning point d component and the first flux turning point q component to obtain corresponding absolute values ​​of errors; if all the absolute values ​​of the errors are less than the allowable error value, output the second rotor angular velocity and the second rotor angle of the synchronous machine; Processing the three-phase voltage to obtain the second current q component, the second current d component and the rotor current of the synchronous machine armature current includes: The three-phase voltage is converted by Park transformation to obtain a dq0-axis voltage component corresponding to the three-phase voltage; According to the matrix parameters of the Thevenin equation on the stator side and the dq0-axis voltage component, the second current q component and the second current d component of the armature current of the synchronous machine are calculated by the armature current calculation formula; Based on the parameter data of the synchronous machine, the dq0-axis voltage component, the second current q component and the second current d component, the rotor current of the synchronous machine is calculated by a rotor current calculation formula; Wherein, the Parker transform is: ； The armature current calculation formula is: ； The rotor current calculation formula is: ； ； ； ； In the formula, is the first current d component, is the first current q component, is the second current d component, is the second current q component, R d , R q R0 and R1 are the resistance parameters of the resistance matrix in the Thevenin equation on the stator side, e d 、e q , e0 are the voltage parameters of the voltage source matrix in the Thevenin equation on the stator side, is the first rotor angle, v a is the phase a voltage of the three-phase voltage, v b is the b-phase voltage of the three-phase voltage, v c is the c-phase voltage of the three-phase voltage, v d is the first voltage d component of the dq0 axis voltage component, v q is the second voltage q component of the dq0 axis voltage component, v0 is the third voltage 0 component of the dq0 axis voltage component, and the parameter data of the synchronous machine include the saturation parameter b of the synchronous machine, the excitation current i f , direct axis damping winding D current i D , quadrature-axis damping winding g current i g and the quadrature-axis damping winding Q current i Q ,i r is the rotor current matrix, is the stator self-inductance dq0 matrix of the synchronous machine, R s is the stator resistance matrix of the synchronous machine, k is 2 / , is the stator-rotor mutual inductance dq0 matrix of the synchronous machine, , and are the stator current, stator voltage and stator flux phase domain matrices of the previous time step respectively.

10. A synchronous machine electromagnetic transient modeling device based on saturation effect, characterized in that: including a processor and a memory; The memory is used to store program codes and transmit the program codes to the processor; The processor is used to execute the synchronous machine electromagnetic transient modeling method based on saturation effect as described in any one of claims 1-8 according to the instructions in the program code.

Images