Reduced-order parallel acceleration calculation method for finite element model of electromagnetic field of relay

Abstract

The invention discloses a relay electromagnetic field finite element model reduced-order parallel acceleration calculation method, and belongs to the technical field of relays. Constructing electromagnetic dynamic characteristics and an electromagnetic field equation of the relay; establishing a relay electromagnetic finite element model; a transmission line iteration method is introduced, a high-order nonlinear equation of the electromagnetic finite element model is decoupled into a high-order linear equation and a plurality of low-order nonlinear equations, and parallel acceleration calculation of the low-order nonlinear equations of the electromagnetic finite element model is achieved; a relay improved eigenorthogonal decomposition model order reduction parallel finite element method is provided, order reduction is carried out on a high-order linear equation after the electromagnetic finite element model is decoupled, and rapid calculation of the electromagnetic finite element model of the relay is achieved. According to the method, the calculation and analysis efficiency of the electromagnetic field characteristics of the relay is improved, the design and research and development cycle of the relay is shortened, a new method is provided for rapid calculation of the electromagnetic characteristics of the relay, and intelligent optimization and high-end intelligent manufacturing of the relay can be effectively assisted.

Description

Technical Field

[0001] The invention relates to a parallel accelerated calculation method for reducing the order of a relay electromagnetic field finite element model, and belongs to the technical field of relays. Background Art

[0002] Relay is an electrical control device with functions such as amplification, isolation, and protection. It is mainly used to realize automatic control, remote control, and conversion control of circuits. In view of the complex and changeable working environment of relays, it is particularly important to optimize the static and dynamic characteristics, temperature performance, mechanical resistance, life, and quality consistency of relays according to actual conditions before production.

[0003] In order to achieve the optimal design of relays, the modern power system design and optimization process is highly dependent on the modeling and simulation of relay performance based on advanced intelligent algorithms. This link requires the ability to accurately simulate the dynamic response of the relay under various working conditions to ensure the efficiency and practicality of the model, thereby providing solid data support for subsequent optimization work. In this context, the finite element method has become the mainstream technology in the field of relay performance simulation with its powerful ability to handle complex electromagnetic field distribution and accurately calculate the electromagnetic force inside the relay.

[0004] The finite element method discretizes the relay structure into a series of small units, applies electromagnetic field theories such as Maxwell's equations, and conducts a detailed analysis of each unit, and finally integrates the performance parameters of the entire system. This method has significant advantages in ensuring simulation accuracy, can accurately reflect the behavioral characteristics of the relay in a real working environment, and provides a reliable theoretical basis for the optimal design of the relay.

[0005] However, despite the excellent performance of the finite element method in terms of solution accuracy, its computational efficiency is relatively low, which is mainly attributed to the huge computing resources required for detailed analysis of a large number of units. In the actual process of relay optimization design, the intelligent optimization algorithm needs to solve the target performance multiple times in order to find the global optimal solution. This requirement is in sharp conflict with the high computational cost of the finite element method, resulting in the overall optimization process taking too long and being unable to meet the urgent needs of modern industry for fast design and real-time response.

[0006] Therefore, how to improve the calculation efficiency while maintaining the high accuracy of the finite element method, or explore other simulation methods that can balance accuracy and efficiency, has become an important issue that needs to be solved in the current relay simulation calculation field. The solution to this problem is not only related to the performance improvement and cost reduction of relay products, but also the key to promoting the intelligent and efficient development of power systems. Summary of the invention

[0007] In order to solve the problems existing in the background technology, the present invention provides a parallel accelerated calculation method for reducing the order of a relay electromagnetic field finite element model.

[0008] To achieve the above object, the present invention adopts the following technical solution: a method for parallel acceleration calculation of a relay electromagnetic field finite element model with reduced order, the method comprising the following steps:

[0009] S1: Construct the electromagnetic dynamic characteristics and electromagnetic field equations of the relay;

[0010] S101: According to the working principle of the relay, the finite element method is used to solve the vector magnetic potential A. The electromagnetic field equation of the relay is as follows:

[0011]

[0012] n·B=0 (2)

[0013] n×H=0 (3)

[0014] In formula (1):

[0015] represents the sign of partial derivative;

[0016] v represents magnetic permeability;

[0017] A represents the vector magnetic potential;

[0018] σ s represents electrical conductivity;

[0019] Js represents the current density applied by the outside world;

[0020] ▽ represents the vector differential operator, which is used to calculate the divergence of the magnetic induction intensity B;

[0021] B r It indicates the remanence of permanent magnet;

[0022] x is a unit vector, indicating the direction of the x-axis;

[0023] y is a unit vector, indicating the direction of the y-axis;

[0024] z is a unit vector, indicating the direction of the z-axis;

[0025] t represents time;

[0026] Equations (2) and (3) are the boundary conditions of the electromagnetic field of the relay, where:

[0027] n is a unit vector, indicating the normal direction of the interface, surface or boundary;

[0028] B represents the magnetic induction intensity;

[0029] H represents the magnetic field strength;

[0030] S102: According to the vector magnetic potential A, the magnetic flux Ψ is calculated as follows:

[0031]

[0032] In formula (4):

[0033] S represents the area enclosed by the coil;

[0034] l represents the boundary of the area enclosed by the coil;

[0035] S103: Solve the electromagnetic force F1 on the armature in the magnetic field region according to the virtual work method:

[0036]

[0037] In formula (5):

[0038] F1 represents the electromagnetic force on the armature;

[0039] W represents the magnetic field energy;

[0040] V represents the magnetic field region;

[0041] u represents the displacement direction;

[0042] S104: Construct the electromagnetic dynamic characteristics and electromagnetic field equations of the relay, and solve the dynamic characteristics of the armature according to the electromagnetic dynamic characteristics equations of the relay:

[0043] in:

[0044] The electromagnetic dynamic characteristic equations of direct-acting relays are as follows:

[0045]

[0046] In formula (6):

[0047] Ψ represents magnetic flux;

[0048] U represents the coil excitation voltage;

[0049] i represents the coil current;

[0050] R represents the coil resistance;

[0051] F l Indicates the electromagnetic force on the armature;

[0052] F f Indicates the reaction force of the contact spring system on the armature;

[0053] m represents the armature mass;

[0054] x1 represents the armature displacement;

[0055] The electromagnetic dynamic characteristic equations of the rotary relay are as follows:

[0056]

[0057] In formula (7):

[0058] T l It represents the electromagnetic torque on the armature;

[0059] T f Indicates the reaction torque of the contact spring system on the armature;

[0060] I x It represents the moment of inertia of the armature;

[0061] θ represents the armature rotation angle.

[0062] S2: Establish the electromagnetic finite element model of the relay;

[0063] S201: Each solution area of ​​the relay is arbitrarily divided into multiple independent triangular areas. For any triangular area Ω e The three vertices are numbered K, M and N respectively;

[0064] S202: Representing the triangular area Ω by linear interpolation e The vector magnetic potential function A at any point inside e as follows:

[0065]

[0066] In formula (8):

[0067] x t and t All represent vertex coordinates;

[0068] N i represents the shape function;

[0069] A i Represents the vector magnetic potential of the vertex, i = K, M, N;

[0070] Δ e Represents the area of ​​a triangle;

[0071] p i ,q i and r i Both represent shape function N i The coefficient of

[0072]

[0073] In formula (9):

[0074] x K and K Respectively represent the x-axis and y-axis coordinates of point K;

[0075] x M and M Respectively represent the x-axis and y-axis coordinates of point M;

[0076] x N and N Respectively represent the x-axis and y-axis coordinates of point N.

[0077] S203: Using the natural boundary conditions through the Galerkin method, set the residual R e With the weighting function W e The relationship is as follows:

[0078]

[0079] In formula (10):

[0080] R e Represents the residual number;

[0081] W e represents the weighting function;

[0082] v e represents the magnetic permeability of the triangular area;

[0083] σ se represents the conductivity of the triangular area;

[0084] Js e represents the externally applied current density in the triangular area;

[0085] S204: Set the residual number R e = 0, the weighting function W e Set it as the shape function and integrate it, we can get:

[0086]

[0087] In formula (11):

[0088] m n Indicates the number of divided triangular areas;

[0089] S205: Arrange equation (11) into a matrix form according to the vertex numbers as follows:

[0090]

[0091] In formula (12):

[0092] F c represents the current density matrix;

[0093] K c It is expressed by formula (11) The matrix is ​​accumulated and calculated according to the vertex number to obtain the coefficient matrix;

[0094] M c It is expressed by formula (11) The matrix is ​​accumulated and calculated according to the vertex number to obtain the coefficient matrix;

[0095] represents the column matrix of the vector magnetic potential of each vertex, ni represents the number of vertices to be found;

[0096] S206: After solving the vector magnetic potential according to equation (12), the magnetic induction intensity B in the triangular area can be solved. e :

[0097]

[0098] S207: Solve the flux linkage according to equations (4) and (8):

[0099]

[0100] In formula (14):

[0101] m coil Indicates the number of coil units;

[0102] N c Indicates the number of coil turns;

[0103] S coil Indicates the cross-sectional area of ​​the coil;

[0104] S208: Solve the electromagnetic force on the relay armature according to equation (15):

[0105]

[0106] m ch Indicates the number of vertices of the contour of the moving part;

[0107] S208: Simulating the electromagnetic performance of the relay according to the magnetic induction intensity, magnetic flux and electromagnetic force on the armature;

[0108] S209: Since the vector magnetic potential and current density of the two-dimensional electromagnetic field only contain the component in the z direction, equation (1) can be transformed into:

[0109]

[0110] In formula (16):

[0111] Js zrepresents the current density in the z direction;

[0112] S2010: Using the axisymmetric model, according to the cylindrical coordinate system transformation, the vector magnetic potential and current density are perpendicular to the zr plane, and considering the boundary conditions, we can get:

[0113]

[0114] In formula (17):

[0115] A'=rA represents the transformed vector magnetic potential;

[0116] v'=v / r represents the transformed magnetic permeability;

[0117] z is a unit vector, indicating the direction of the z-axis;

[0118] r is a unit vector, indicating the direction of the r-axis;

[0119] Γ1 represents the boundary condition 1 of equation (17);

[0120] Γ2 represents the second boundary condition of equation (17);

[0121] S2011: Due to the axisymmetry, after replacing the z-axis and r-axis of the cylindrical coordinate system with the x-axis and y-axis, equation (17) is the same as equation (16), so finite element modeling and solution can be performed according to equation (16).

[0122] S3: Introducing the transmission line iteration method, the high-order nonlinear equations of the electromagnetic finite element model of S2 are decoupled into high-order linear equations and several low-order nonlinear equations, realizing the parallel accelerated calculation of the low-order nonlinear equations of the electromagnetic finite element model;

[0123] S301: Set the conductivity and capacitance parameters as follows. The vector magnetic potential difference of each vertex is the voltage value:

[0124]

[0125] In formula (18):

[0126] G KM represents the distributed conductance between vertex K and vertex M;

[0127] G NK represents the distributed conductance between vertex N and vertex K;

[0128] G MN represents the distributed conductance between vertex M and vertex N;

[0129] C KM represents the distributed capacitance between vertex K and vertex M;

[0130] C NKrepresents the distributed capacitance between vertex N and vertex K;

[0131] C MN represents the distributed capacitance between vertex M and vertex N;

[0132] C K represents the distributed capacitance of vertex K;

[0133] C M represents the distributed capacitance of vertex M;

[0134] C N represents the distributed capacitance of vertex C;

[0135] σ e Indicates the electrical conductivity of the component material;

[0136] S302: adding transmission lines at both ends of the nonlinear element;

[0137] S303: Perform Norton equivalence on the circuit containing the transmission line, and use the resistor equivalent capacitor with a time step of Δt. The capacitor equivalent admittance parameters are as follows:

[0138]

[0139] In formula (19):

[0140] G CKM represents the distributed conductance between vertex K and vertex M in the equivalent Norton circuit;

[0141] G CNK represents the distributed conductance between vertex N and vertex K in the equivalent Norton circuit;

[0142] G CMN represents the distributed conductance between vertex M and vertex N in the equivalent Norton circuit;

[0143] G CK represents the distributed conductance of vertex K in the equivalent Norton circuit;

[0144] G CM represents the distributed conductance of vertex M in the equivalent Norton circuit;

[0145] G CN represents the distributed conductance of vertex N in the equivalent Norton circuit;

[0146] S304: Obtain the current source through capacitor equivalent, and complete the calculation of any triangular area Ω in the finite element model e The circuit equivalent of:

[0147]

[0148] In formula (20):

[0149] I CKM represents the current source current between vertex K and vertex M in the equivalent Norton circuit;

[0150] I CMN represents the current source current between vertex M and vertex N in the equivalent Norton circuit;

[0151] I CNK represents the current source current between vertex N and vertex K in the equivalent Norton circuit;

[0152] I CK represents the current source current at vertex K in the equivalent Norton circuit;

[0153] I CM represents the current source current at vertex M in the equivalent Norton circuit;

[0154] I CN represents the current source current at vertex N in the equivalent Norton circuit;

[0155] Indicates the voltage value of vertex K at the previous moment;

[0156] represents the voltage value of vertex M at the previous moment;

[0157] Indicates the voltage value of vertex N at the previous moment;

[0158] S305: Based on the Euler method and equation (12), the characteristic admittance parameters of the transmission line are set as follows through the initial value of the equivalent capacitance of the current source:

[0159]

[0160] In formula (21):

[0161] Y GKM represents the admittance between vertex K and vertex M;

[0162] Y GNK represents the admittance between vertex N and vertex K;

[0163] Y GMN represents the admittance between vertex M and vertex N;

[0164] v eg Represents the guessed magnetic permeability. The closer the guessed magnetic permeability is to the real magnetic permeability, the faster the transmission line iteration converges.

[0165] S306: Assume that the nth t The voltage between nodes at the time of the transmission line iteration is U X (n t ):

[0166] U X (n t )=U i (n t )+U r (n t ) (twenty two)

[0167] U r (n t ) indicates the nth t The reflected voltage at the first transmission line iteration is the input value;

[0168] U i (n t ) indicates the nth t The incident voltage at the transmission line iteration;

[0169] S307: Calculate the nth t +1 transmission line iteration with incident voltage U i (n t +1),

[0170] Since all triangular areas share the same magnetic permeability v e Therefore, when calculating the incident voltage, the nonlinear resistance in all triangular areas must satisfy equation (23) at the same time, and equation (24) must be solved at the same time:

[0171]

[0172] In formula (23):

[0173] G represents the characteristic conductance of the transmission line;

[0174] Z represents the characteristic impedance of the transmission line;

[0175] In formula (24):

[0176] V x , V y , V z is the incident wave in the x, y, and z directions;

[0177] V x0 , V y0 , V z0 It is the reflected wave in the x, y, and z directions;

[0178] S308: Solve equation (24) by Newton's method to calculate the derivative of the required conductance with respect to the incident wave:

[0179]

[0180] S309: Calculate the current source current according to the voltage difference between the two ends of the node:

[0181]

[0182] In formula (26):

[0183] I GKM (n t +1) indicates the nth t +1 transmission line iteration Current of the current source between nodes K and M;

[0184] I GMN (n t +1) indicates the nth t +1 transmission line iteration current source between nodes M and N;

[0185] I GNK (n t +1) indicates the nth t +1 transmission line iteration Current sourced between nodes N and K.

[0186] U KM (n t +1) indicates the nth t +1 transmission line iteration voltage difference between node K and node M;

[0187] U MN (n t +1) indicates the nth t +1 transmission line iteration voltage difference between node M and node N;

[0188] U NK (n t +1) indicates the nth t +1 transmission line iteration voltage difference between node N and node K;

[0189] S3010: Integrating the parameters of equations (18) to (26) can achieve the purpose of solving equation (12). At this time, equation (12) becomes a linear equation (27):

[0190] K ctlm A+M ctlm A=F c +F tlm +F Alast (27)

[0191] In formula (27):

[0192] K ctlm Represents the original coefficient matrix K c The linear part and the characteristic admittance in equation (21) are calculated according to the corresponding vertices and;

[0193] M ctlmIn expression (19), all capacitive admittances are expressed by corresponding vertices and;

[0194] F tlm The current sources in the Norton equivalent circuit representing the transmission line are arranged in corresponding vertices and;

[0195] F Alast The current sources equivalently generated by the capacitors in expression (20) are summed up at each vertex.

[0196] S4: A parallel finite element method for reducing the order of the improved intrinsic orthogonal decomposition model of the relay is proposed to reduce the order of the high-order linear equations after decoupling the S3 electromagnetic finite element model, thereby realizing the rapid calculation of the relay electromagnetic finite element model.

[0197] S401: Based on the eigenorthogonal decomposition theory, equation (27) is reduced to equation (28):

[0198]

[0199] In formula (28):

[0200] K ctlmr It represents the characteristic admittance after reduction according to the corresponding vertices and;

[0201] A r represents the reduced-order vector magnetic potential;

[0202] M ctlmr It indicates that the capacitance admittance after reduction is in accordance with the corresponding vertices and;

[0203] F tlmr Represent the current sources in the reduced-order Norton equivalent circuit according to the corresponding vertices and;

[0204] F Alastr The current sources equivalent to the capacitors after the reduction are summed up according to the vertices;

[0205] Ψ P represents an orthogonal matrix;

[0206] Denotes the orthogonal matrix Ψ P The transposed matrix of

[0207] F r Represents the reduced orthogonal matrix Ψ P The transposed matrix of With the current density matrix F c The product of

[0208] S402: In each transmission line iteration process, according to formula (29), the reduced vector magnetic potential A r Calculate the unreduced vector magnetic potential A,

[0209] xp =Ψ p x pl (29)

[0210] In formula (29):

[0211] x p Represents a high-dimensional space vector, x p ∈R n ;

[0212] x pl represents a low-dimensional space vector, x pl ∈R k ;

[0213] Current density matrix F c The current source generated by the equivalent capacitor is summed up at each vertex F Alast will change at each time step, so it needs to be solved again at each time step;

[0214] The current sources in the Norton equivalent circuit are arranged according to the corresponding vertices and F tlm It will change after each transmission line iteration, so it needs to be recalculated after each transmission line iteration;

[0215] S402: Further calculate formula (22):

[0216] When the characteristic impedance of the transmission line does not change, that is, when the magnetic permeability does not change, the characteristic admittance is calculated according to the corresponding vertex and K ctlm and the capacitance admittance according to the corresponding vertex and M ctlm All parameters are unchanged, and only one reduction calculation is required.

[0217] Compared with the prior art, the present invention has the following beneficial effects:

[0218] The present invention provides a specific simulation method for the electromagnetic force of a relay based on mathematical equations, and clarifies that the solution of high-order nonlinear equations is the main factor affecting the calculation efficiency of the finite element method; introduces a transmission line iteration method, and decouples the electromagnetic high-order nonlinear equations into high-order linear equations and several low-order nonlinear equations, thereby realizing parallel calculation of the low-order nonlinear equations and improving the calculation efficiency of the nonlinear equations of the relay electromagnetic finite element model; determines that the transmission line iteration method has a bottleneck of low efficiency in solving high-order linear equations, and proposes an improved intrinsic orthogonal decomposition model reduced-order parallel finite element method, which accelerates the solution efficiency of the high-order linear equations of the relay electromagnetic characteristics and solves the nonlinear equations in parallel, thereby improving the calculation and analysis efficiency of the relay electromagnetic field characteristics, shortening the relay design and development cycle, and providing a new method for the rapid solution of the relay electromagnetic characteristics, which can effectively assist the intelligent optimization and high-end intelligent manufacturing of relays. BRIEF DESCRIPTION OF THE DRAWINGS

[0219] Figure 1: is a schematic diagram of the equivalent circuit of the finite element unit of the present invention, wherein: (a) is a schematic diagram of a triangular area; (b) is a schematic diagram of the equivalent circuit of the triangular area; (c) is a schematic diagram of the equivalent Norton circuit;

[0220] Figure 2 It is a flow chart of the transmission line iterative finite element method;

[0221] Figure 3 It is a flowchart of the reduced-order parallel finite element method;

[0222] Figure 4 This is a schematic diagram of the structure of a commercial relay, where: 1 is a static contact, 2 is a moving contact, 3 is an overtravel spring, 4 is a connecting rod, 5 is a reaction spring, 6 is an armature, 7 is an iron core, 8 is a coil, and 9 is a yoke;

[0223] Figure 5 The magnetic vector potential cloud diagram comparison diagram when the armature is in the release position, where: (a) is the Newton iteration comparison diagram, (b) is the transmission line iteration comparison diagram, and (c) is the reduced-order parallel comparison diagram;

[0224] Figure 6 The magnetic vector potential cloud diagrams when the armature is in the attracted position are compared, where: (a) is the Newton iteration comparison diagram, (b) is the transmission line iteration comparison diagram, and (c) is the reduced-order parallel comparison diagram. DETAILED DESCRIPTION

[0225] The technical solution 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 described embodiments are only part of the embodiments of the invention, rather than all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of the present invention.

[0226] A method for parallel accelerated calculation of a relay electromagnetic field finite element model with reduced order, comprising the following steps:

[0227] S1: Construct the electromagnetic dynamic characteristics and electromagnetic field equations of the relay;

[0228] When the relay is working, current flows through the coil to generate a magnetic field. The change in current causes the change in the magnetic field, and the electromagnetic force on the armature between the iron core and the armature also changes accordingly. Under the action of the changing electromagnetic force on the armature, the armature will overcome the reaction force provided by the contact spring system and move (translation or rotation). The movement of the armature will drive the moving contact to move. When the moving contact moves to contact or separate from the static contact, the output circuit jump is realized.

[0229] The electromagnetic dynamic characteristics of the relay are important indicators for describing the working process of the relay, which are mainly determined by two parts:

[0230] One part is the electromagnetic force F exerted by the control circuit and electromagnetic structure on the relay armature. l And the influence of magnetic flux Ψ;

[0231] Part of it is the electromagnetic force F on the armature l The reaction force F of the spring system on the armature f The above two parts are coupled with each other and can be described by the electromagnetic field equation and electromagnetic dynamic characteristic equation group of the relay. The electromagnetic relay has a complex structure and is mainly divided into two types according to its action form: direct-acting relay and rotary relay.

[0232] S101: According to the working principle of the relay, the finite element method is used to solve the vector magnetic potential A. The electromagnetic field equation of the relay is as follows:

[0233]

[0234] n·B=0 (2)

[0235] n×H=0 (3)

[0236] In formula (1):

[0237] represents the sign of partial derivative;

[0238] v represents magnetic permeability;

[0239] A represents the vector magnetic potential;

[0240] σ s represents electrical conductivity;

[0241] Js represents the current density applied by the outside world;

[0242] ▽ represents the vector differential operator, which is used to calculate the divergence of the magnetic induction intensity B;

[0243] B r It indicates the remanence of permanent magnet;

[0244] x is a unit vector, indicating the direction of the x-axis;

[0245] y is a unit vector, indicating the direction of the y-axis;

[0246] z is a unit vector, indicating the direction of the z-axis;

[0247] t represents time;

[0248] Equations (2) and (3) are the boundary conditions of the electromagnetic field of the relay, where:

[0249] n is a unit vector, indicating the normal direction of the interface, surface or boundary;

[0250] B represents the magnetic induction intensity;

[0251] H represents the magnetic field strength;

[0252] S102: According to the vector magnetic potential A, the magnetic flux Ψ is calculated as follows:

[0253]

[0254] In formula (4):

[0255] S represents the area enclosed by the coil;

[0256] l represents the boundary of the area enclosed by the coil;

[0257] S103: Solve the electromagnetic force F1 on the armature in the magnetic field region according to the virtual work method:

[0258]

[0259] In formula (5):

[0260] F1 represents the electromagnetic force on the armature;

[0261] W represents the magnetic field energy;

[0262] V represents the magnetic field region;

[0263] u represents the displacement direction;

[0264] S104: Construct the electromagnetic dynamic characteristics and electromagnetic field equations of the relay, and solve the dynamic characteristics of the armature according to the electromagnetic dynamic characteristics equations of the relay:

[0265] in:

[0266] The electromagnetic dynamic characteristic equations of direct-acting relays are as follows:

[0267]

[0268] In formula (6):

[0269] Ψ represents magnetic flux;

[0270] U represents the coil excitation voltage;

[0271] i represents the coil current;

[0272] R represents the coil resistance;

[0273] F l Indicates the electromagnetic force on the armature;

[0274] F f Indicates the reaction force of the contact spring system on the armature;

[0275] m represents the armature mass;

[0276] x1 represents the armature displacement;

[0277] The electromagnetic dynamic characteristic equations of the rotary relay are as follows:

[0278]

[0279] In formula (7):

[0280] T l It represents the electromagnetic torque on the armature;

[0281] T f Indicates the reaction torque of the contact spring system on the armature;

[0282] I x It represents the moment of inertia of the armature;

[0283] θ represents the armature rotation angle.

[0284] S2: Establish the electromagnetic finite element model of the relay;

[0285] S201: Each solution region of the relay is arbitrarily divided into multiple independent, non-overlapping triangular regions. For any triangular region Ω e The three vertices are numbered K, M and N respectively;

[0286] The solution areas related to the relay and the electromagnetic field mainly include the core area, yoke area, air area, coil area and armature area. The conductivity of the core area, yoke area and armature area is a nonlinear area, and the conductivity of the remaining areas is a linear area.

[0287] S202: Representing the triangular area Ω by linear interpolation e The vector magnetic potential function A at any point inside e as follows:

[0288]

[0289] In formula (8):

[0290] x t and t All represent vertex coordinates;

[0291] N i Represents the shape function, which is only related to the vertex coordinates;

[0292] A i Represents the vector magnetic potential of the vertex, i = K, M, N;

[0293] Δ e Represents the area of ​​a triangle;

[0294] pi ,q i and r i Both represent shape function N i The coefficient of is related to the vertex coordinates;

[0295]

[0296] In formula (9):

[0297] x K and K Respectively represent the x-axis and y-axis coordinates of point K;

[0298] x M and M Respectively represent the x-axis and y-axis coordinates of point M;

[0299] x N and N Respectively represent the x-axis and y-axis coordinates of point N.

[0300] S203: Using the natural boundary conditions through the Galerkin method, set the residual R e With the weighting function W e The relationship is as follows:

[0301]

[0302] In formula (10):

[0303] R e Represents the residual number;

[0304] W e represents the weighting function;

[0305] v e represents the magnetic permeability of the triangular area;

[0306] σ se represents the conductivity of the triangular area;

[0307] Js e represents the externally applied current density in the triangular area;

[0308] S204: Set the residual number R e = 0, the weighting function W e Set it as the shape function and integrate it, we can get:

[0309]

[0310] In formula (11):

[0311] m n Indicates the number of divided triangular areas;

[0312] S205: Arrange equation (11) into a matrix form according to the vertex numbers as follows:

[0313]

[0314] In formula (12):

[0315] F c represents the current density matrix;

[0316] K c It is expressed by formula (11) The matrix is ​​accumulated and calculated according to the vertex number to obtain the coefficient matrix;

[0317] M c It is expressed by formula (11) The matrix is ​​accumulated and calculated according to the vertex number to obtain the coefficient matrix;

[0318] represents the column matrix of the vector magnetic potential of each vertex, ni represents the number of vertices to be found;

[0319] S206: After solving the vector magnetic potential according to equation (12), the magnetic induction intensity B in the triangular area can be solved. e :

[0320]

[0321] S207: Solve the flux linkage according to equations (4) and (8):

[0322]

[0323] In formula (14):

[0324] m coil Indicates the number of coil units;

[0325] N c Indicates the number of coil turns;

[0326] S coil Indicates the cross-sectional area of ​​the coil;

[0327] S208: Solve the electromagnetic force on the relay armature according to equation (15):

[0328]

[0329] m ch Indicates the number of vertices of the contour of the moving part;

[0330] S208: Simulating the electromagnetic performance of the relay according to the magnetic induction intensity, magnetic flux and electromagnetic force on the armature;

[0331] S209: Since the vector magnetic potential and current density of the two-dimensional electromagnetic field only contain the component in the z direction, equation (1) can be transformed into:

[0332]

[0333] In formula (16):

[0334] Js z represents the current density in the z direction;

[0335] S2010: Using the axisymmetric model according to the cylindrical coordinate system transformation, the vector magnetic potential and current density are perpendicular to the zr plane, and considering the boundary conditions, we can get:

[0336]

[0337] In formula (17):

[0338] A'=rA represents the transformed vector magnetic potential;

[0339] v'=v / r represents the transformed magnetic permeability;

[0340] z is a unit vector, indicating the direction of the z-axis;

[0341] r is a unit vector, indicating the direction of the r-axis;

[0342] Γ1 represents the boundary condition 1 of equation (17);

[0343] Γ2 represents the second boundary condition of equation (17);

[0344] S2011: Due to the axisymmetry, after replacing the z-axis and r-axis of the cylindrical coordinate system with the x-axis and y-axis, equation (17) is the same as equation (16), so finite element modeling and solution can be performed according to equation (16).

[0345] S3: Introducing the transmission line iteration method, the high-order nonlinear equations of the electromagnetic finite element model of S2 are decoupled into high-order linear equations and several low-order nonlinear equations, realizing the parallel accelerated calculation of the low-order nonlinear equations of the electromagnetic finite element model;

[0346] Solving S2 in equation (12) is the core of solving the finite element model. Due to the nonlinearity of magnetic permeability, the coefficient matrix K in equation (12) is c It will change nonlinearly with the magnetic induction intensity, so equation (12) is an N-order nonlinear equation. Common commercial finite element software uses the Newton-Raphson method to solve this equation, but due to the high order of the nonlinear equation, this algorithm is very inefficient. Here, the transmission line iteration method is introduced to improve the calculation efficiency.

[0347] S301: Set the conductivity and capacitance parameters as follows: Figure 1The vector magnetic potential difference of each vertex is the voltage value:

[0348]

[0349] In formula (18):

[0350] G KM represents the distributed conductance between vertex K and vertex M;

[0351] G NK represents the distributed conductance between vertex N and vertex K;

[0352] G MN represents the distributed conductance between vertex M and vertex N;

[0353] C KM represents the distributed capacitance between vertex K and vertex M;

[0354] C NK represents the distributed capacitance between vertex N and vertex K;

[0355] C MN represents the distributed capacitance between vertex M and vertex N;

[0356] C K represents the distributed capacitance of vertex K;

[0357] C M represents the distributed capacitance of vertex M;

[0358] C N represents the distributed capacitance of vertex C;

[0359] σ e Indicates the electrical conductivity of the component material;

[0360] S302: The magnetic permeability v of the triangular region in equation (18) e Transmission lines are added at both ends of the nonlinear element. Figure 1 (b) is shown by the thick black solid line;

[0361] S303: Perform Norton equivalence on the circuit containing the transmission line, and use the resistor equivalent capacitor with a time step of Δt. The capacitor equivalent admittance parameters are as follows:

[0362]

[0363] In formula (19):

[0364] G CKM represents the distributed conductance between vertex K and vertex M in the equivalent Norton circuit;

[0365] G CNKrepresents the distributed conductance between vertex N and vertex K in the equivalent Norton circuit;

[0366] G CMN represents the distributed conductance between vertex M and vertex N in the equivalent Norton circuit;

[0367] G CK represents the distributed conductance of vertex K in the equivalent Norton circuit;

[0368] G CM represents the distributed conductance of vertex M in the equivalent Norton circuit;

[0369] G CN represents the distributed conductance of vertex N in the equivalent Norton circuit;

[0370] S304: Obtain the current source through capacitor equivalent, and complete the calculation of any triangular area Ω in the finite element model e like Figure 1 As shown in (a), Figure 1 The circuit equivalent of (b) is:

[0371]

[0372] In formula (20):

[0373] I CKM represents the current source current between vertex K and vertex M in the equivalent Norton circuit;

[0374] I CMN represents the current source current between vertex M and vertex N in the equivalent Norton circuit;

[0375] I CNK represents the current source current between vertex N and vertex K in the equivalent Norton circuit;

[0376] I CK represents the current source current at vertex K in the equivalent Norton circuit;

[0377] I CM represents the current source current at vertex M in the equivalent Norton circuit;

[0378] I CN represents the current source current at vertex N in the equivalent Norton circuit;

[0379] Indicates the voltage value of vertex K at the previous moment;

[0380] represents the voltage value of vertex M at the previous moment;

[0381] Indicates the voltage value of vertex N at the previous moment;

[0382] S305: Based on the Euler method and equation (12), the initial value of the equivalent capacitance of the current source is calculated. Figure 1 (b) is converted into Figure 1 (c). Figure 1 The characteristic admittance parameters of the transmission line are set as follows:

[0383]

[0384] In formula (21):

[0385] Y GKM represents the admittance between vertex K and vertex M;

[0386] Y GNK represents the admittance between vertex N and vertex K;

[0387] Y GMN represents the admittance between vertex M and vertex N;

[0388] v eg Represents the guessed magnetic permeability. The closer the guessed magnetic permeability is to the real magnetic permeability, the faster the transmission line iteration converges.

[0389] S306: Assume that the nth t The voltage between nodes at the time of the transmission line iteration is U X (n t ), which is determined by the incident voltage U i (n t ) and reflected voltage U r (n t ) jointly determine, thus we can get the nth t Secondary reflected voltage U r (n t ) as formula (22):

[0390] U X (n t )=U i (n t )+U r (n t ) (twenty two)

[0391] U r (n t ) indicates the nth t The reflected voltage at the first transmission line iteration is the input value;

[0392] U i (n t ) indicates the nth t The incident voltage at the transmission line iteration;

[0393] S307: Calculate the nth t+1 transmission line iteration with incident voltage U i (n t +1),

[0394] Since all triangular areas share the same magnetic permeability v e Therefore, when calculating the incident voltage, the nonlinear resistance in all triangular areas must satisfy equation (23) at the same time, and equation (24) must be solved at the same time, but each triangular area can be calculated in parallel.

[0395]

[0396] In formula (23):

[0397] G represents the characteristic conductance of the transmission line;

[0398] Z represents the characteristic impedance of the transmission line;

[0399] In formula (24):

[0400] V x , V y , V z is the incident wave in the x, y, and z directions;

[0401] V x0 , V y0 , V z0 It is the reflected wave in the x, y, and z directions;

[0402] S308: Solve equation (24) by Newton's method to calculate the derivative of the required conductance with respect to the incident wave:

[0403]

[0404] S309: Figure 1 In (c), the current source current is calculated based on the voltage difference across the node:

[0405]

[0406] In formula (26):

[0407] I GKM (n t +1) indicates the nth t +1 transmission line iteration Current of the current source between nodes K and M;

[0408] I GMN (n t +1) indicates the nth t +1 transmission line iteration current source between nodes M and N;

[0409] I GNK (nt +1) indicates the nth t +1 transmission line iteration Current sourced between nodes N and K.

[0410] U KM (n t +1) indicates the nth t +1 transmission line iteration voltage difference between node K and node M;

[0411] U MN (n t +1) indicates the nth t +1 transmission line iteration voltage difference between node M and node N;

[0412] U NK (n t +1) indicates the nth t +1 transmission line iteration voltage difference between node N and node K;

[0413] S3010: Integrating the parameters of equations (18) to (26) can achieve the purpose of solving equation (12). At this time, equation (12) becomes a linear equation (27):

[0414] K ctlm A+M ctlm A=F c +F tlm +F Alast (27)

[0415] In formula (27):

[0416] K ctlm Represents the original coefficient matrix K c The linear part and the characteristic admittance in equation (21) are calculated according to the corresponding vertices and;

[0417] M ctlm In expression (19), all capacitive admittances are expressed by corresponding vertices and;

[0418] F tlm The current sources in the Norton equivalent circuit representing the transmission line are arranged in corresponding vertices and;

[0419] F Alast The current sources equivalently generated by the capacitors in expression (20) are summed up at each vertex.

[0420] The establishment and solution process of formula (27) is as follows Figure 2 As shown, here t max To maximize the solution time, the nonlinear units in the figure can be solved in parallel, which greatly improves the solution efficiency of the finite element method.

[0421] S4: A parallel finite element method for reducing the order of the improved intrinsic orthogonal decomposition model of the relay is proposed to reduce the order of the high-order linear equations after decoupling the S3 electromagnetic finite element model, thereby realizing the rapid calculation of the relay electromagnetic finite element model.

[0422] The transmission line iteration method in S3 is very effective in guessing the magnetic permeability v eg The value of the real magnetic permeability v e When the difference is large, the convergence speed will be very slow, and more iterations are required to obtain the final result. That is, the transmission line iteration method needs to solve many high-order linear equations (27). Although the solution efficiency of linear equations is higher than that of nonlinear equations, it still consumes a lot of time due to the increase in the number of iterations. This problem leads to the fact that the transmission line iteration method cannot fully satisfy users in improving the efficiency of finite elements. This problem can be improved by the following methods:

[0423] (1) To speed up the convergence of transmission line iteration, the estimated magnetic permeability v can be introduced eg Update mechanism, update the guessed permeability v after a certain number of iterations eg , which makes it closer to the true magnetic permeability in subsequent iterations, thereby reducing the number of solutions to the high-order linear equations.

[0424] (2) Speed ​​up the solution of linear equations, because the magnetic permeability v is guessed eg When the characteristic admittance does not change, it is calculated by the corresponding vertex and K ctlm and the capacitance admittance according to the corresponding vertex and M ctlm There is no change. We can process these two matrices in advance through LU decomposition and then solve the nonlinear equations, thereby improving the speed of solving the linear equations in a single time. However, the time of LU decomposition also needs to be considered. At the same time, the high-order characteristic admittance is calculated according to the corresponding vertex and K ctlm and the capacitance admittance according to the corresponding vertex and M ctlm It will greatly occupy parallel device resources.

[0425] The above two methods are contradictory, and the updated guess of magnetic permeability v eg The more times, the better the convergence. It seems that the linear equation solution time can be reduced by reducing the number of linear equation solution times. However, at this time, the characteristic admittance is calculated according to the corresponding vertex and K ctlm and the capacitance admittance according to the corresponding vertex and M ctlm The more times of LU decomposition are required, the longer the linear equation solving time may be. Therefore, this patent further improves the efficiency of the finite element method from this perspective, and proposes a parallel finite element method based on the improved intrinsic orthogonal decomposition model of the relay to solve the problem of time-consuming electromagnetic finite element model in the linear equation solving process. The specific steps are as follows:

[0426] S401: Based on the eigenorthogonal decomposition theory, equation (27) is reduced to equation (28),

[0427]

[0428] In formula (28):

[0429] K ctlmr It represents the characteristic admittance after reduction according to the corresponding vertices and;

[0430] A r represents the reduced-order vector magnetic potential;

[0431] M ctlmr It indicates that the capacitance admittance after reduction is in accordance with the corresponding vertices and;

[0432] F tlmr Represent the current sources in the reduced-order Norton equivalent circuit according to the corresponding vertices and;

[0433] F Alastr The current sources equivalent to the capacitors after the reduction are summed up according to the vertices;

[0434] Ψ P represents an orthogonal matrix;

[0435] Denotes the orthogonal matrix Ψ P The transposed matrix of

[0436] F r Represents the reduced orthogonal matrix Ψ P The transposed matrix of With the current density matrix F c The product of

[0437] S402: Improve the overall process of the parallel finite element method for the reduction of the eigen-orthogonal decomposition model. Figure 3 shown.

[0438] Depend on Figure 3 It can be seen that in each transmission line iteration process, according to formula (29), the reduced vector magnetic potential A r Calculate the unreduced vector magnetic potential A,

[0439] x p =Ψ p x pl (29)

[0440] In formula (29):

[0441] x p Represents a high-dimensional space vector, x p ∈R n ;

[0442] xpl represents a low-dimensional space vector, x pl ∈R k ;

[0443] Current density matrix F c The current source generated by the equivalent capacitor is summed up at each vertex F Alast will change at each time step, so it needs to be solved again at each time step;

[0444] The current sources in the Norton equivalent circuit are arranged according to the corresponding vertices and F tlm It will change after each transmission line iteration, so it needs to be recalculated after each transmission line iteration;

[0445] S402: Further calculate formula (22):

[0446] When the characteristic impedance of the transmission line does not change, that is, when the magnetic permeability does not change, the characteristic admittance is calculated according to the corresponding vertex and K ctlm and the capacitance admittance according to the corresponding vertex and M ctlm All parameters are unchanged, and only one order reduction calculation is required, which has a low impact on the overall computational complexity.

[0447] Embodiment 1:

[0448] Select a commercial relay (such as Figure 4 An example analysis is performed as shown).

[0449] The electromagnetic field of the relay is solved by the present invention, the transmission line iterative finite element method and the Newton iterative finite element method, and the calculation accuracy and calculation time of the vector magnetic potential and the electromagnetic force are compared.

[0450] When comparing the accuracy of the vector magnetic potential, the results obtained using the commercial software COMSOL were used as a benchmark. This software uses Newton iteration to solve the nonlinear finite element equations.

[0451] When comparing the electromagnetic force accuracy, a relay static characteristic electromagnetic force test device (existing technology, no longer described) is used to measure the electromagnetic force of the relay as a benchmark, and the electromagnetic force at the position where the armature is attracted and released is selected for actual measurement to compare and verify with the results of the above three algorithms.

[0452] (1) The electromagnetic force of the relay armature in the energized position is simulated. At this time, the number of meshes in the finite element model is 36095, the number of nodes is 18124, the sample data used in the model of the present invention is 32, and the order selected for simulation is 16. The simulation results are shown in Figure 5 As shown, the vector magnetic potential clouds obtained by the three methods are Figure 1 Based on the COMSOL calculation results, the average deviation of each node of the three algorithms is within 0.60%.

[0453] (2) The electromagnetic force of the relay armature in the released position is simulated. At this time, the number of meshes in the finite element model is 37010, the number of nodes is 18582, the sample data used in the model of the present invention is 32, and the order selected for simulation is 16. The simulation results are shown in Figure 6 As shown, the vector magnetic potential clouds obtained by the three methods are Figure 1 Based on the COMSOL calculation results, the average deviation of each node of the three algorithms is within 0.90%.

[0454] (3) The comparison between the electromagnetic force results at rated voltage calculated by the three algorithms and the measured results is shown in Table 1. The error between the algorithm results and the measured results is shown in brackets. The errors of the three algorithms in the table are all small.

[0455] Table 1 Performance comparison of transmission line iterative finite element method, Newton iterative finite element method and the present invention

[0456]

[0457] The present invention reduces the order of the linear matrix and combines the advantages of the transmission line iterative finite element method and the Newton iterative finite element method. By comparison, the calculation efficiency of the present invention for calculating the attraction electromagnetic characteristics is 8 times that of the Newton iterative finite element method and 2 times that of the transmission line iterative finite element method, and the calculation efficiency of the release electromagnetic characteristics is 7 times that of the Newton iterative finite element method and 2 times that of the transmission line iterative finite element method, which greatly improves the calculation efficiency while ensuring the calculation accuracy.

[0458] It will be apparent to those skilled in the art that the invention is not limited to the details of the exemplary embodiments described above and that the invention can be implemented in other forms without departing from the spirit or essential features of the invention. Therefore, the embodiments should be considered in all respects as exemplary and non-limiting, and the scope of the invention is defined by the appended claims rather than the foregoing description, and it is intended that all variations within the meaning and range of equivalents of the claims be included in the invention. Any reference numeral in a claim should not be considered as limiting the claim to which it relates.

[0459] In addition, it should be understood that although the present specification is described according to implementation modes, not every implementation mode contains only one independent technical solution. This description of the specification is only for the sake of clarity. Those skilled in the art should regard the specification as a whole. The technical solutions in each embodiment may also be appropriately combined to form other implementation modes that can be understood by those skilled in the art.

Claims

1. A parallel accelerated calculation method for reducing the order of a relay electromagnetic field finite element model, characterized in that: The method comprises the following steps: S1: Construct the electromagnetic dynamic characteristics and electromagnetic field equations of the relay; S2: Establish the electromagnetic finite element model of the relay; S3: Introducing the transmission line iteration method, the high-order nonlinear equations of the electromagnetic finite element model of S2 are decoupled into high-order linear equations and several low-order nonlinear equations, realizing the parallel accelerated calculation of the low-order nonlinear equations of the electromagnetic finite element model; S4: A parallel finite element method for reducing the order of the improved intrinsic orthogonal decomposition model of the relay is proposed to reduce the order of the high-order linear equations after decoupling the S3 electromagnetic finite element model, thereby realizing the rapid calculation of the relay electromagnetic finite element model.

2. The method for parallel accelerated calculation of a relay electromagnetic field finite element model with reduced order according to claim 1, characterized in that: The S1 comprises the following steps: S101: According to the working principle of the relay, the finite element method is used to solve the vector magnetic potential A. The electromagnetic field equation of the relay is as follows: n·B=0 (2) n×H=0 (3) In formula (1): represents the sign of partial derivative; v represents magnetic permeability; A represents the vector magnetic potential; σ s represents electrical conductivity; Js represents the current density applied by the outside world; represents the vector differential operator, which is used to calculate the divergence of the magnetic induction intensity B; B r It indicates the remanence of permanent magnet; x is a unit vector, indicating the direction of the x-axis; y is a unit vector, indicating the direction of the y-axis; z is a unit vector, indicating the direction of the z-axis; t represents time; Equations (2) and (3) are the boundary conditions of the electromagnetic field of the relay, where: n is a unit vector, indicating the normal direction of the interface, surface or boundary; B represents the magnetic induction intensity; H represents the magnetic field strength; S102: According to the vector magnetic potential A, the magnetic flux Ψ is calculated as follows: In formula (4): S represents the area enclosed by the coil; l represents the boundary of the area enclosed by the coil; S103: Solve the electromagnetic force F1 on the armature in the magnetic field region according to the virtual work method: In formula (5): F1 represents the electromagnetic force on the armature; W represents the magnetic field energy; V represents the magnetic field region; u represents the displacement direction; S104: Construct the electromagnetic dynamic characteristics and electromagnetic field equations of the relay, and solve the dynamic characteristics of the armature according to the electromagnetic dynamic characteristics equations of the relay: in: The electromagnetic dynamic characteristic equations of direct-acting relays are as follows: In formula (6): Ψ represents magnetic flux; U represents the coil excitation voltage; i represents the coil current; R represents the coil resistance; F l Indicates the electromagnetic force on the armature; F f Indicates the reaction force of the contact spring system on the armature; m represents the armature mass; x1 represents the armature displacement; The electromagnetic dynamic characteristic equations of the rotary relay are as follows: In formula (7): T l It represents the electromagnetic torque on the armature; T f Indicates the reaction torque of the contact spring system on the armature; I x It represents the moment of inertia of the armature; θ represents the armature rotation angle.

3. The method for parallel accelerated calculation of a relay electromagnetic field finite element model with reduced order according to claim 2, characterized in that: The S2 comprises the following steps: S201: Each solution area of ​​the relay is arbitrarily divided into multiple independent triangular areas. For any triangular area Ω e The three vertices are numbered K, M and N respectively; S202: Representing the triangular area Ω by linear interpolation e The vector magnetic potential function A at any point inside e as follows: In formula (8): x t and t All represent vertex coordinates; N i represents the shape function; A i Represents the vector magnetic potential of the vertex, i = K, M, N; Δ e Represents the area of ​​a triangle; p i ,q i and r i Both represent shape function N i The coefficient of In formula (9): x K and K Respectively represent the x-axis and y-axis coordinates of point K; x M and M Respectively represent the x-axis and y-axis coordinates of point M; x N and N Respectively represent the x-axis and y-axis coordinates of point N. S203: Using the natural boundary conditions through the Galerkin method, set the residual R e With the weighting function W e The relationship is as follows: In formula (10): R e Display remaining number; W e represents the weighting function; v e represents the magnetic permeability of the triangular area; σ se represents the conductivity of the triangular area; Js e represents the externally applied current density in the triangular area; S204: Set the residual number R e = 0, the weighting function W e Assuming it as shape function and integrating it, we can get: In formula (11): m n Indicates the number of divided triangular areas; S205: Arrange equation (11) into a matrix form according to the vertex numbers as follows: In formula (12): F c represents the current density matrix; K c It is expressed by formula (11) The matrix is ​​accumulated and calculated according to the vertex number to obtain the coefficient matrix; M c It is expressed by formula (11) The matrix is ​​accumulated and calculated according to the vertex number to obtain the coefficient matrix; represents the column matrix of the vector magnetic potential of each vertex, ni represents the number of vertices to be found; S206: After solving the vector magnetic potential according to equation (12), the magnetic induction intensity B in the triangular area can be solved. e : S207: Solve the flux linkage according to equations (4) and (8): In formula (14): m coil Indicates the number of coil units; N c Indicates the number of coil turns; S coil Indicates the cross-sectional area of ​​the coil; S208: Solve the electromagnetic force on the relay armature according to equation (15): m ch Indicates the number of vertices of the contour of the moving part; S208: Simulating the electromagnetic performance of the relay according to the magnetic induction intensity, magnetic flux and electromagnetic force on the armature; S209: Since the vector magnetic potential and current density of the two-dimensional electromagnetic field only contain the component in the z direction, equation (1) can be transformed into: In formula (16): Js z represents the current density in the z direction; S2010: Using the axisymmetric model, according to the cylindrical coordinate system transformation, the vector magnetic potential and current density are perpendicular to the zr plane, and considering the boundary conditions, we can get: In formula (17): A'=rA represents the transformed vector magnetic potential; v'=v / r represents the transformed magnetic permeability; z is a unit vector, indicating the direction of the z-axis; r is a unit vector, indicating the direction of the r-axis; Γ1 represents the boundary condition 1 of equation (17); Γ2 represents the second boundary condition of equation (17); S2011: Due to the axisymmetry, after replacing the z-axis and r-axis of the cylindrical coordinate system with the x-axis and y-axis, equation (17) is the same as equation (16), so finite element modeling and solution can be performed according to equation (16).

4. The method for parallel accelerated calculation of a relay electromagnetic field finite element model with reduced order according to claim 3, characterized in that: The S3 comprises the following steps: S301: Set the conductivity and capacitance parameters as follows. The vector magnetic potential difference of each vertex is the voltage value: In formula (18): G KM represents the distributed conductance between vertex K and vertex M; G NK represents the distributed conductance between vertex N and vertex K; G MN represents the distributed conductance between vertex M and vertex N; C KM represents the distributed capacitance between vertex K and vertex M; C NK represents the distributed capacitance between vertex N and vertex K; C MN represents the distributed capacitance between vertex M and vertex N; C K represents the distributed capacitance of vertex K; C M represents the distributed capacitance of vertex M; C N represents the distributed capacitance of vertex C; σ e Indicates the electrical conductivity of the component material; S302: adding transmission lines at both ends of the nonlinear element; S303: Perform Norton equivalence on the circuit containing the transmission line, and use the resistor equivalent capacitor with a time step of Δt. The capacitor equivalent admittance parameters are as follows: In formula (19): G CKM represents the distributed conductance between vertex K and vertex M in the equivalent Norton circuit; G CNK represents the distributed conductance between vertex N and vertex K in the equivalent Norton circuit; G CMN represents the distributed conductance between vertex M and vertex N in the equivalent Norton circuit; G CK represents the distributed conductance of vertex K in the equivalent Norton circuit; G CM represents the distributed conductance of vertex M in the equivalent Norton circuit; G CN represents the distributed conductance of vertex N in the equivalent Norton circuit; S304: Obtain the current source through capacitor equivalent, and complete the calculation of any triangular area Ω in the finite element model e The circuit equivalent of: In formula (20): I CKM represents the current source current between vertex K and vertex M in the equivalent Norton circuit; I CMN represents the current source current between vertex M and vertex N in the equivalent Norton circuit; I CNK represents the current source current between vertex N and vertex K in the equivalent Norton circuit; I CK represents the current source current at vertex K in the equivalent Norton circuit; I CM represents the current source current at vertex M in the equivalent Norton circuit; I CN represents the current source current at vertex N in the equivalent Norton circuit; Indicates the voltage value of vertex K at the previous moment; represents the voltage value of vertex M at the previous moment; Represents the voltage value of vertex N at the previous moment; S305: Based on the Euler method and equation (12), the characteristic admittance parameters of the transmission line are set as follows through the initial value of the equivalent capacitance of the current source: In formula (21): Y GKM represents the admittance between vertex K and vertex M; Y GNK represents the admittance between vertex N and vertex K; Y GMN represents the admittance between vertex M and vertex N; v eg Represents the guessed magnetic permeability. The closer the guessed magnetic permeability is to the real magnetic permeability, the faster the transmission line iteration converges. S306: Assume that the nth t The voltage between nodes at the time of the transmission line iteration is U X (n t ): U X (n t )=U i (n t )+U r (n t ) (22) U r (n t ) indicates the nth t The reflected voltage at the first transmission line iteration is the input value; U i (n t ) indicates the nth t The incident voltage at the transmission line iteration; S307: Calculate the nth t +1 transmission line iteration with incident voltage U i (n t +1), Since all triangular areas share the same magnetic permeability v e Therefore, when calculating the incident voltage, the nonlinear resistance in all triangular areas must satisfy equation (23) at the same time, and equation (24) must be solved at the same time: In formula (23): G represents the characteristic conductance of the transmission line; Z represents the characteristic impedance of the transmission line; In formula (24): V x , V y , V z It is the incident wave in the x, y, and z directions; V x0 , V y0 , V z0 It is the reflected wave in the x, y, and z directions; S308: Solve equation (24) by Newton's method to calculate the derivative of the required conductance with respect to the incident wave: S309: Calculate the current source current according to the voltage difference between the two ends of the node: In formula (26): I GKM (n t +1) indicates the nth t +1 transmission line iteration Current of the current source between nodes K and M; I GMN (n t +1) indicates the nth t +1 transmission line iteration current source between nodes M and N; I GNK (n t +1) indicates the nth t +1 transmission line iteration Current sourced between nodes N and K. U KM (n t +1) indicates the nth t +1 transmission line iteration voltage difference between node K and node M; U MN (n t +1) indicates the nth t +1 transmission line iteration voltage difference between node M and node N; U NK (n t +1) indicates the nth t +1 transmission line iteration voltage difference between node N and node K; S3010: Integrating the parameters of equations (18) to (26) can achieve the purpose of solving equation (12). At this time, equation (12) becomes a linear equation (27): K ctlm A+M ctlm A=F c +F tlm +F Alast (27) In formula (27): K ctlm Represents the original coefficient matrix K c The linear part and the characteristic admittance in equation (21) are calculated according to the corresponding vertices and; M ctlm In expression (19), all capacitive admittances are expressed by corresponding vertices and; F tlm The current sources in the Norton equivalent circuit representing the transmission line are arranged in corresponding vertices and; F Alast The current sources equivalently generated by the capacitors in expression (20) are summed up at each vertex.

5. The method for parallel accelerated calculation of a relay electromagnetic field finite element model with reduced order according to claim 4, characterized in that: The S4 comprises the following steps: S401: Based on the eigenorthogonal decomposition theory, equation (27) is reduced to equation (28), In formula (28): K ctlmr It represents the characteristic admittance after reduction according to the corresponding vertices and; A r represents the reduced-order vector magnetic potential; M ctlmr It indicates that the capacitance admittance after reduction is in terms of the corresponding vertices and; F tlmr Represent the current sources in the reduced-order Norton equivalent circuit according to the corresponding vertices and; F Alastr The current sources equivalent to the capacitors after the reduction are summed up according to the vertices; Ψ P represents an orthogonal matrix; Denotes the orthogonal matrix Ψ P The transposed matrix of F r Represents the reduced orthogonal matrix Ψ P The transposed matrix of With the current density matrix F c The product of S402: In each transmission line iteration process, according to formula (29), the reduced vector magnetic potential A r Calculate the unreduced vector magnetic potential A, x p =Ψ p x pl (29) In formula (29): x p Represents a high-dimensional space vector, x p ∈R n ; x pl represents a low-dimensional space vector, x pl ∈R k ; Current density matrix F c The current source generated by the equivalent capacitor is summed up at each vertex F Alast will change at each time step, so it needs to be solved again at each time step; The current sources in the Norton equivalent circuit are arranged according to the corresponding vertices and F tlm It will change after each transmission line iteration, so it needs to be recalculated after each transmission line iteration; S402: Further calculate formula (22): When the characteristic impedance of the transmission line does not change, that is, when the magnetic permeability does not change, the characteristic admittance is calculated according to the corresponding vertex and K ctlm and the capacitance admittance according to the corresponding vertex and M ctlm All parameters are unchanged, and only one reduction calculation is required.

Images