A Method for Reducing Order and Parallel Accelerated Calculation of the Finite Element Model of Relay Electromagnetic Field

Abstract

A method for reducing the order and parallel accelerating the calculation of the finite element model of the electromagnetic field of a relay belongs to the technical field of relays. The electromagnetic dynamic characteristics and electromagnetic field equations of the relay are constructed; an electromagnetic finite element model of the relay is established; the transmission line iteration method is introduced to decouple the high-order nonlinear equations of the electromagnetic finite element model into high-order linear equations and several low-order nonlinear equations, realizing the parallel acceleration calculation of the low-order nonlinear equations of the electromagnetic finite element model; a method for reducing the order and parallel finite element method of the improved proper orthogonal decomposition model of the relay is proposed to reduce the order of the high-order linear equations after decoupling the electromagnetic finite element model, realizing the rapid calculation of the electromagnetic finite element model of the relay. The present invention improves the calculation and analysis efficiency of the electromagnetic field characteristics of the relay, shortens the design and research and development cycle of the relay, provides a new method for the rapid solution of the electromagnetic characteristics of the relay, and can effectively assist the intelligent optimization and high-end intelligent manufacturing of the relay.

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] To solve the problems existing in the background art, the present invention provides a method for reducing the order and parallel accelerating the calculation of the relay electromagnetic field finite element model.

[0008] To achieve the above object, the present invention adopts the following technical solutions: A method for reducing the order and parallel accelerating the calculation of the relay electromagnetic field finite element model, the method comprising the following steps:

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

[0010] S101: According to the working principle of the relay, use the finite element method to solve the vector magnetic potential A, and 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 partial derivative symbol;

[0016] v represents the magnetic permeability;

[0017] A represents the vector magnetic potential;

[0018] σ s Represents the conductivity;

[0019] Js represents the externally applied current density;

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

[0021] B r Represents the remanence of the permanent magnet;

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

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

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

[0025] t represents time;

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

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

[0028] B represents the magnetic induction intensity;

[0029] H represents the magnetic field intensity;

[0030] S102: Calculate the magnetic flux linkage Ψ according to the vector magnetic potential A as follows:

[0031]

[0032] In Equation (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 for the electromagnetic force F1 acting on the armature in the magnetic field region according to the virtual work method:

[0036]

[0037] In Equation (5):

[0038] F1 represents the electromagnetic force acting 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 for the dynamic characteristics of the armature according to the electromagnetic dynamic characteristic equations of the relay:

[0043] Where:

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

[0045]

[0046] In Equation (6):

[0047] Ψ represents the magnetic flux linkage;

[0048] U represents the coil excitation voltage;

[0049] i represents the coil current;

[0050] R represents the coil resistance;

[0051] F l represents the electromagnetic force acting on the armature;

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

[0053] m represents the mass of the armature;

[0054] x1 represents the displacement of the armature;

[0055] The electromagnetic dynamic characteristic equations of a rotating relay are as follows:

[0056]

[0057] In Equation (7):

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

[0059] T f represents the counter torque of the contact spring system on the armature;

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

[0061] θ represents the rotation angle of the armature.

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

[0063] S201: Arbitrarily divide each solution region of the relay into multiple independent triangular regions. For any triangular region Ω e the three vertices are numbered K, M, and N respectively;

[0064] S202: Represent the vector magnetic potential function A e at any point inside the triangular region Ω e as follows:

[0065]

[0066] In Equation (8):

[0067] x t and y t both represent the 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 the triangular region;

[0071] p i 、q i and r i all represent the coefficients of the shape function N i ;

[0072]

[0073] In Equation (9):

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

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

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

[0077] S203: Using the Galerkin method with natural boundary conditions, set the residue R e and the weighting function W e to have the following relationship:

[0078]

[0079] In equation (10):

[0080] R e represents the residue;

[0081] W e represents the weighting function;

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

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

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

[0085] S204: Set the residue R e = 0, and set the weighting function W e as the shape function for integration, and we can obtain:

[0086]

[0087] In equation (11):

[0088] m n represents the number of divided triangular regions;

[0089] S205: According to the vertex numbers, organize equation (11) into the matrix form as follows:

[0090]

[0091] In equation (12):

[0092] F c represents the current density matrix;

[0093] K c represents the coefficient matrix obtained by cumulative calculation of the matrix in Equation (11) according to the vertex numbers;

[0094] M c represents the coefficient matrix obtained by cumulative calculation of the matrix in Equation (11) according to the vertex numbers;

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

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

[0097]

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

[0099]

[0100] In Equation (14):

[0101] m coil represents the number of coil elements;

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

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

[0104] S208: Solve the electromagnetic force received by the relay armature according to Equation (15):

[0105]

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

[0107] S208: Simulate the electromagnetic performance of the relay according to the magnetic induction intensity, magnetic flux linkage and electromagnetic force received by the armature;

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

[0109]

[0110] In Equation (16):

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

[0112] S2010: By using the axisymmetric model and transforming according to the cylindrical coordinate system, the vector magnetic potential and the current density are perpendicular to the z-r plane. Considering the boundary conditions, we can obtain:

[0113]

[0114] In Equation (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 representing the direction of the z-axis;

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

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

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

[0121] S2011: Due to 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). Therefore, the finite element modeling and solution can be carried out according to Equation (16).

[0122] S3: Introduce the transmission line iteration method to decouple the high-order nonlinear equations of the electromagnetic finite element model in S2 into high-order linear equations and several low-order nonlinear equations, and realize the parallel acceleration calculation of the low-order nonlinear equations of the electromagnetic finite element model;

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

[0124]

[0125] In Equation (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 Represents the conductivity of the component material;

[0136] S302: Add a transmission line at both ends of the nonlinear component;

[0137] S303: Perform Norton equivalent on the circuit containing the transmission line. Through the resistance equivalent capacitance with time step Δt, the capacitance equivalent admittance parameters are as follows:

[0138]

[0139] In Equation (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 capacitance equivalence to complete the circuit equivalence of any triangular region Ω e in the finite element model:

[0147]

[0148] In Equation (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 of vertex K in the equivalent Norton circuit;

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

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

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

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

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

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

[0159]

[0160] In Equation (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 permeability. The closer the guessed permeability is to the true permeability, the faster the transmission line iteration convergence speed;

[0165] S306: Assume that the voltage between nodes during the n t th transmission line iteration is U X (n t ):

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

[0167] U r (n t ) represents the reflected voltage during the nth t transmission line iteration, which is the input value;

[0168] U i (n t ) represents the incident voltage during the nth t transmission line iteration;

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

[0170] Since all triangular regions share the same magnetic permeability v of a triangular region e , when calculating the incident voltage, the nonlinear resistors in all triangular regions must simultaneously satisfy Equation (23), and at the same time, Equation (24) needs to be solved:

[0171]

[0172] In Equation (23):

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

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

[0175] In Equation (24):

[0176] V x , V y , V z are the incident waves in the x, y, and z directions;

[0177] V x0 , V y0 , V z0 are the reflected waves in the x, y, and z directions;

[0178] S308: Solve Equation (24) by the Newton 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 across the node:

[0181]

[0182] In formula (26):

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

[0184] I GMN (n t +1) represents the current of the current source between the n t +1-th transmission line iteration nodes M and N;

[0185] I GNK (n t +1) represents the current of the current source between the n t +1-th transmission line iteration nodes N and K.

[0186] U KM (n t +1) represents the voltage difference between the n t +1-th transmission line iteration nodes K and M;

[0187] U MN (n t +1) represents the voltage difference between the n t +1-th transmission line iteration nodes M and N;

[0188] U NK (n t +1) represents the voltage difference between the n t +1-th transmission line iteration nodes N and K;

[0189] S3010: Integrating the parameters of formulas (18) to (26) can achieve the purpose of solving formula (12). At this time, formula (12) becomes the 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 linear part and the sum of the characteristic admittances at the corresponding vertices in formula (21);

[0193] M ctlmThe sum of all capacitance admittances in Expression (19) is taken according to the corresponding vertices;

[0194] F tlm The sum of the current sources in the Norton equivalent circuit of the transmission line is taken according to the corresponding vertices;

[0195] F Alast The current sources generated by the equivalent capacitance in Expression (20) are accumulated and summed according to the vertices.

[0196] S4: Propose a relay-improved proper orthogonal decomposition model reduction parallel finite element method to reduce the order of the high-order linear equations decoupled from the electromagnetic finite element model in S3, and realize the rapid calculation of the relay electromagnetic finite element model.

[0197] S401: Based on the proper orthogonal decomposition theory, Expression (27) is reduced to Expression (28):

[0198]

[0199] In Expression (28):

[0200] K ctlmr Represents the sum of the reduced characteristic admittances according to the corresponding vertices;

[0201] A r Represents the reduced vector magnetic potential;

[0202] M ctlmr Represents the sum of the reduced capacitance admittances according to the corresponding vertices;

[0203] F tlmr Represents the sum of the reduced current sources in the Norton equivalent circuit according to the corresponding vertices;

[0204] F Alastr Represents the sum of the reduced current sources generated by the equivalent capacitance by vertex accumulation;

[0205] Ψ P Represents the orthogonal matrix;

[0206] Represents the orthogonal matrix Ψ P of the transpose matrix;

[0207] F r Represents the transpose matrix of the reduced orthogonal matrix Ψ P of and the current density matrix F c of the product;

[0208] S402: In each iteration process of the transmission line, according to Expression (29), the unreduced vector magnetic potential A is calculated from the reduced vector magnetic potential A r ,

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

[0210] In Equation (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] The current density matrix F c and the current source generated by the equivalent capacitance are summed up according to the vertices, and F Alast will change at each time step, so it needs to be re-solved at each time step;

[0214] The current source in the Norton equivalent circuit is summed up according to the corresponding vertices and F tlm will change after each transmission line iteration. Therefore, it needs to be recalculated after each transmission line iteration;

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

[0216] When the characteristic impedance of the transmission line does not change, that is, when it is guessed that the magnetic permeability does not change, the characteristic admittance is summed up according to the corresponding vertices and K ctlm and the capacitance admittance is summed up according to the corresponding vertices and M ctlm are all invariant parameters, and a single order reduction calculation can be performed.

[0217] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0218] The present invention provides a specific simulation method for the electromagnetic force of a relay according to mathematical physics equations, and clarifies that the solution of high-order nonlinear equations is the main factor affecting the calculation efficiency of the finite element method; the transmission line iteration method is introduced, and by decoupling the electromagnetic high-order nonlinear equations into high-order linear equations and several low-order nonlinear equations, the parallel calculation of low-order nonlinear equations is realized, and the calculation efficiency of the nonlinear equations of the relay electromagnetic finite element model is improved; it is determined that there is a bottleneck in the low efficiency of solving high-order linear equations in the transmission line iteration method, and an improved proper orthogonal decomposition model reduction parallel finite element method is proposed, which accelerates the solution efficiency of high-order linear equations of relay electromagnetic characteristics while parallelly solving nonlinear equations, improves the calculation and analysis efficiency of relay electromagnetic field characteristics, shortens the relay design and R & D cycle, provides a new method for rapid calculation of relay electromagnetic characteristics, and can effectively assist relay intelligent optimization and high-end intelligent manufacturing. BRIEF DESCRIPTION OF THE DRAWINGS

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

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

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

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

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

[0224] Figure 6 It is the comparison of the magnetic vector potential nephogram when the armature is in the attracted position, where: (a) is the comparison diagram of Newton iteration, (b) is the comparison diagram of transmission line iteration, and (c) is the comparison diagram of reduced-order parallel. Specific implementation mode

[0225] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the invention, rather than all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

[0226] A method for reducing the order and accelerating the parallel calculation of the finite element model of the electromagnetic field of a relay, including the following steps:

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

[0228] When the relay works, the current flows through the coil to generate a magnetic field. The change of the current causes the change of 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 (translate or rotate). The movement of the armature will drive the movement of the moving contact. When the moving contact moves to contact or separate from the static contact, the jump of the output circuit is realized.

[0229] The electromagnetic dynamic characteristics of the relay are an important index to describe the working process of the relay, which is mainly determined by two parts:

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

[0231] One part is the influence of the electromagnetic force F l acting on the armature and the reaction force F f of the contact spring system acting on the armature on the motion state of the armature. The above two parts are mutually coupled and can be described by the electromagnetic field equations and electromagnetic dynamic characteristic equations of the relay. The electromagnetic relay has a complex structure. According to its action form, there are mainly two types: direct-acting relay and rotating relay.

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

[0233]

[0234] n·B = 0 (2)

[0235] n×H = 0 (3)

[0236] In Equation (1):

[0237] represents the partial derivative symbol;

[0238] v represents the magnetic permeability;

[0239] A represents the vector magnetic potential;

[0240] σ s represents the conductivity;

[0241] Js represents the externally applied current density;

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

[0243] B r represents the remanence of the permanent magnet;

[0244] x is the unit vector representing the x-axis direction;

[0245] y is the unit vector representing the y-axis direction;

[0246] z is the unit vector representing the z-axis direction;

[0247] t represents time;

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

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

[0250] B represents magnetic induction intensity;

[0251] H represents magnetic field strength;

[0252] S102: Calculate the magnetic flux linkage Ψ according to the vector magnetic potential A 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 acting 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 acting on the armature;

[0261] W represents 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 characteristic equations of the relay:

[0265] Where:

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

[0267]

[0268] In formula (6):

[0269] Ψ represents magnetic flux linkage;

[0270] U represents the coil excitation voltage;

[0271] i represents the coil current;

[0272] R represents the coil resistance;

[0273] F l represents the electromagnetic force acting on the armature;

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

[0275] m represents the mass of the armature;

[0276] x1 represents the armature displacement;

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

[0278]

[0279] In Equation (7):

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

[0281] T f represents the counter torque of the contact spring system on the armature;

[0282] I x 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: Arbitrarily divide each solution region of the relay into multiple independent and non - overlapping triangular regions. For any triangular region Ω e the three vertices are numbered K, M, and N respectively;

[0286] The solution regions of the relay related to the electromagnetic field mainly include the iron core region, yoke region, air region, coil region, and armature region. Among them, the conductivities of the iron core region, yoke region, and armature region are non - linear regions, and the conductivities of the remaining regions are linear regions.

[0287] S202: Represent the vector magnetic potential function A e at any point inside the triangular region Ω e as follows:

[0288]

[0289] In Equation (8):

[0290] x t and y t both represent the 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 the triangular region;

[0294] pi , q i and r i both represent the coefficients of the shape function N i , which are related to the vertex coordinates;

[0295]

[0296] In equation (9):

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

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

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

[0300] S203: Using the Galerkin method with natural boundary conditions, set the residue R e and the weighting function W e to have the following relationship:

[0301]

[0302] In equation (10):

[0303] R e represents the residue;

[0304] W e represents the weighting function;

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

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

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

[0308] S204: Set the residue R e = 0, and set the weighting function W e as the shape function for integration, and we can obtain:

[0309]

[0310] In equation (11):

[0311] m n represents the number of divided triangular regions;

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

[0313]

[0314] In Equation (12):

[0315] F c represents the current density matrix;

[0316] K c represents the coefficient matrix obtained by cumulative calculation of the matrix in Equation (11) according to the vertex numbers;

[0317] M c represents the coefficient matrix obtained by cumulative calculation of the matrix in Equation (11) according to the vertex numbers;

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

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

[0320]

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

[0322]

[0323] In Equation (14):

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

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

[0326] S coil represents 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 represents the number of vertices of the contour of the moving part;

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

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

[0332]

[0333] In Equation (16):

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

[0335] S2010: Using the axisymmetric model and according to the transformation in the cylindrical coordinate system, the vector magnetic potential and current density are perpendicular to the z-r plane. Considering the boundary conditions, we can obtain:

[0336]

[0337] In Equation (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 representing the z-axis direction;

[0341] r is a unit vector representing the r-axis direction;

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

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

[0344] S2011: Since after the axisymmetric transformation replaces 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 we can perform finite element modeling and solution according to Equation (16).

[0345] S3: Introduce the transmission line iteration method to decouple the high-order nonlinear equation of the electromagnetic finite element model in S2 into a high-order linear equation and several low-order nonlinear equations, and realize the parallel acceleration calculation of the low-order nonlinear equations of the electromagnetic finite element model;

[0346] Solving Equation (12) in S2 is the core of solving the finite element model. Due to the nonlinearity of the magnetic permeability, the coefficient matrix K in Equation (12) c will change nonlinearly with the magnetic induction intensity. Therefore, Equation (12) is an N-order nonlinear equation. The Newton-Raphson method is used to solve this equation in common commercial finite element software. However, 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 conductance and capacitance parameters as follows, Figure 1The vector magnetic potential difference between the vertices in [object] is the voltage value:

[0348]

[0349] In Equation (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 represents the conductivity of the component material;

[0360] S302: Add transmission lines to both ends of the nonlinear component containing the magnetic permeability v of the triangular region in Equation (18). The transmission lines are shown as the thick black solid lines in e (b); Figure 1 (b) as shown by the thick black solid lines;

[0361] S303: Perform Norton equivalent on the circuit containing the transmission lines. Through the resistance equivalent capacitance with a time step of Δt, the capacitance equivalent admittance parameters are as follows:

[0362]

[0363] In Equation (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 between vertex N and vertex K in the equivalent Norton circuit;

[0370] S304: Obtain a current source through capacitance equivalence to complete any triangular region Ω in the finite element model e such as Figure 1 (a) shows that it can be equivalent to the circuit as Figure 1 (b):

[0371]

[0372] In Equation (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 of vertex K in the equivalent Norton circuit;

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

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

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

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

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

[0382] S305: Based on the Euler's method and Equation (12), the initial value of the equivalent capacitance charge of the current source is obtained, Figure 1 (b) which is then transformed into Figure 1 (c). Figure 1 The characteristic admittance parameters of the transmission line are set as follows:

[0383]

[0384] In Equation (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 permeability. The closer the guessed permeability is to the true permeability, the faster the convergence rate of the transmission line iteration;

[0389] S306: Assume that the voltage between nodes during the n t -th transmission line iteration is U X (n t ). This voltage is jointly determined by the incident voltage U i (n t ) and the reflected voltage U r (n t ). Thus, the reflected voltage U t at the n r -th iteration can be obtained as shown in Equation (22): t U

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

[0391] U r (n t ) represents the reflected voltage during the n t -th transmission line iteration and is the input value;

[0392] U i (n t ) represents the incident voltage during the n t -th transmission line iteration;

[0393] S307: Calculate the n tIncident voltage U during the +1 - th transmission line iteration i (n t +1),

[0394] Since all triangular regions share the same magnetic permeability v of a triangular region e , when calculating the incident voltage, the nonlinear resistors in all triangular regions must simultaneously satisfy Equation (23), and Equation (24) needs to be solved simultaneously. However, parallel calculations can be performed between each triangular region.

[0395]

[0396] In Equation (23):

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

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

[0399] In Equation (24):

[0400] V x , V y , V z are the incident waves in the x, y, and z directions;

[0401] V x0 , V y0 , V z0 are the reflected waves in the x, y, and z directions;

[0402] S308: Solve Equation (24) by the Newton - Raphson method to calculate the derivative of the required conductance with respect to the incident wave:

[0403]

[0404] S309: Figure 1 (c) Calculate the current - source current based on the voltage difference across the node:

[0405]

[0406] In Equation (26):

[0407] I GKM (n t +1) represents the current of the current source between node K and node M during the (n t +1) - th transmission line iteration;

[0408] I GMN (n t +1) represents the current of the current source between node M and node N during the (n t +1) - th transmission line iteration;

[0409] I GNK (nt +1) represents the current of the current source between node N and node K in the n t +1-th transmission line iteration.

[0410] U KM (n t +1) represents the voltage difference between node K and node M in the n t +1-th transmission line iteration;

[0411] U MN (n t +1) represents the voltage difference between node M and node N in the n t +1-th transmission line iteration;

[0412] U NK (n t +1) represents the voltage difference between node N and node K in the n t +1-th transmission line iteration;

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

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

[0415] In equation (27):

[0416] K ctlm represents the original coefficient matrix K c linear part and the sum of the characteristic admittances in equation (21) according to the corresponding vertices;

[0417] M ctlm represents the sum of all capacitance admittances in equation (19) according to the corresponding vertices;

[0418] F tlm represents the sum of the current sources in the Norton equivalent circuit of the transmission line according to the corresponding vertices;

[0419] F Alast represents the sum of the current sources generated by the capacitance equivalence in equation (20) accumulated according to the vertices.

[0420] The establishment and solution process of equation (27) are as Figure 2 shown. Here, t max is the maximum solution time. The non-linear units in the figure can be solved in parallel, greatly improving the solution efficiency of the finite element method.

[0421] S4: Propose the relay improved intrinsic orthogonal decomposition model reduction parallel finite element method to reduce the order of the high-order linear equations decoupled from the electromagnetic finite element model in S3, and achieve the rapid calculation of the relay electromagnetic finite element model.

[0422] In the transmission line iteration method in S3, when the guessed magnetic permeability v eg has a large difference from the true magnetic permeability v e , the convergence speed will be very slow and more iteration times are required to obtain the final result. That is, the transmission line iteration method needs to solve the high-order linear equations (27) many times. Although the solution efficiency of the linear equations is higher than that of the non-linear equations, due to the increase in the number of iterations, this still consumes a lot of time. This problem makes the improvement of the finite element efficiency by the transmission line iteration method not fully satisfactory to users. To address this problem, the following improvements can be made:

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

[0424] (2) To accelerate the solution speed of the linear equations, since the characteristic admittance does not change when the guessed magnetic permeability v eg does not change, corresponding to the vertices and K ctlm and the capacitance admittance does not change corresponding to the vertices and M ctlm either, these two matrices can be processed in advance by LU decomposition and then the non-linear equations are solved, thereby improving the speed of solving the linear equations once. However, the time for LU decomposition also needs to be considered. At the same time, the high-order characteristic admittance corresponding to the vertices and K ctlm and the capacitance admittance corresponding to the vertices and M ctlm will significantly occupy the parallel device resources.

[0425] The above two methods are contradictory. The more times the guessed magnetic permeability v eg is updated, the better the convergence. It seems that the time for solving the linear equations can be reduced by reducing the number of times of solving the linear equations. However, at this time, the characteristic admittance corresponding to the vertices and K ctlm and the capacitance admittance corresponding to the vertices and M ctlm will change, and the number of times of LU decomposition required will also be more, which may instead increase the time for solving the linear equations. Therefore, from this perspective, this patent further improves the efficiency of the finite element method, proposes the relay improved intrinsic orthogonal decomposition model reduction parallel finite element method to solve the problem of time-consuming in the process of solving the linear equations of the electromagnetic finite element model, and the specific steps are as follows:

[0426] S401: Based on the proper orthogonal decomposition theory, after reducing the order of Equation (27), it is transformed into Equation (28).

[0427]

[0428] In Equation (28):

[0429] K ctlmr represents the sum of the reduced characteristic admittances corresponding to the vertices;

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

[0431] M ctlmr represents the sum of the reduced capacitance admittances corresponding to the vertices;

[0432] F tlmr represents the sum of the current sources in the reduced Norton equivalent circuit corresponding to the vertices;

[0433] F Alastr represents the sum of the current sources generated by the reduced capacitance equivalent accumulated by vertices;

[0434] Ψ P represents the orthogonal matrix;

[0435] represents the transpose matrix of the orthogonal matrix Ψ P ;

[0436] F r represents the transpose matrix of the reduced orthogonal matrix Ψ P and the product of the current density matrix F c ;

[0437] S402: The overall process of improving the reduced-order parallel finite element method of the proper orthogonal decomposition model is as Figure 3 shown.

[0438] As Figure 3 can be seen, during each transmission line iteration process, according to Equation (29), the unreduced vector magnetic potential A is calculated from the reduced vector magnetic potential A r ,

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

[0440] In Equation (29):

[0441] x p represents the high-dimensional space vector, x p ∈R n ;

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

[0443] The current density matrix F c and the current source generated by the capacitance equivalent are summed up according to the vertices as F Alast will change at each time step, so it is necessary to solve it again at each time step;

[0444] The current source in the Norton equivalent circuit is summed up according to the corresponding vertices as F tlm will change after each transmission line iteration, so it is necessary to recalculate after each transmission line iteration;

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

[0446] When the characteristic impedance of the transmission line does not change, that is, it is guessed that the magnetic permeability does not change, the characteristic admittance is summed up according to the corresponding vertices as K ctlm and the capacitance admittance is summed up according to the corresponding vertices as M ctlm are all invariant parameters, and a single order reduction calculation can be performed, which has a low impact on the overall computational workload.

[0447] Example 1:

[0448] Select a certain commercial relay (such as Figure 4 shown) for case analysis.

[0449] Solve the electromagnetic field of the relay by the present invention, the transmission line iteration finite element method and the Newton iteration finite element method, and compare the calculation accuracy and calculation time of the vector magnetic potential and the electromagnetic force.

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

[0451] When comparing the accuracy of the electromagnetic force, the electromagnetic force of the relay is actually measured by using the electromagnetic force test device for the static characteristics of the relay (prior art, not elaborated here) as the benchmark, and the electromagnetic force at the positions 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) Simulate the electromagnetic force when the relay armature is in the attracted position. 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 selected order for the simulation is 16th order. The simulation results are as Figure 5 shown, and the vector magnetic potential clouds obtained by the three methods Figure 1 are consistent. Based on the COMSOL calculation results, the average deviation of each node of the three algorithms is within 0.60%.

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

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

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

[0456]

[0457] The present invention reduces the order of the linear matrix, 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 electromagnetic characteristics of the present invention for the closing state is 8 times that of the Newton iterative finite element method and 2 times that of the transmission line iterative finite element method. The calculation efficiency of the electromagnetic characteristics of the released state is 7 times that of the Newton iterative finite element method and 2 times that of the transmission line iterative finite element method. While ensuring the calculation accuracy, the calculation efficiency is greatly improved.

[0458] For those skilled in the art, it is obvious that the present invention is not limited to the details of the above exemplary embodiments, and can be implemented in other forms without departing from the spirit or basic characteristics of the present invention. Therefore, from any point of view, the embodiments should be regarded as exemplary and non-restrictive. The scope of the present invention is defined by the appended claims rather than the above description. Therefore, all changes falling within the meaning and scope of the equivalent conditions of the claims are intended to be included in the present invention. Any reference signs in the claims should not be regarded as limiting the claimed rights.

[0459] In addition, it should be understood that although this specification is described according to the embodiments, not every embodiment only contains an independent technical solution. This narrative way of the specification is only for clarity. Those skilled in the art should regard the specification as a whole, and the technical solutions in each embodiment can also be appropriately combined to form other embodiments that can be understood by those skilled in the art.

Claims

1. A reduced-order parallel acceleration calculation method for the finite element model of a relay electromagnetic field, characterized in that: The method includes 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; The S2 includes the following steps: S201: Arbitrarily divide each solution region of the relay into multiple independent triangular regions. For any triangular region Ω e the three vertices are respectively numbered as K, M, and N; S202: Represent the vector magnetic potential function A of any point inside the triangular region Ω by means of linear interpolation e as follows e as follows: In formula (8): x t and y t both represent vertex coordinates; N i represents a shape function; A i The vector magnetic potential representing the vertex, i = K, M, N; Δ e represents the area of the triangular region; p i ,q i and r i Both represent shape function N i The coefficient of In formula (9): x K and y K respectively represent the x-axis and y-axis coordinates of point K; x M and y M respectively represent the x-axis and y-axis coordinates of point M; x N and y N respectively represent the x-axis and y-axis coordinates of the N points; S203: Set the residue R by using the natural boundary condition through the Galerkin method e and the weighted function W e The relationship is as follows: In formula (10): R e represents the residue; W e represents a weighting function; v e represents the magnetic permeability of the triangular region; σ se represents the conductivity of the triangular region; Js e represents the current density applied to the outside of the triangular region; S204: Set the residue R e = 0, and when the weighted function W e is set as the shape function for integration, we can obtain: In formula (11): m n represents the number of divided triangular regions; S205: According to the vertex numbers, organize formula (11) into the matrix form as follows: In formula (12): F c represents the current density matrix; K c represents the coefficient matrix obtained by cumulative calculation of the matrix in Equation (11) according to the vertex numbers; ​ M c denotes the coefficient matrix obtained by cumulative calculation of the matrix in Equation (11) according to the vertex numbers; ​ A column matrix representing the vector magnetic potential of each vertex, and ni represents the number of vertices to be solved; S3: Introduce the transmission line iteration method to decouple the high-order nonlinear equations of the electromagnetic finite element model in S2 into high-order linear equations and several low-order nonlinear equations, and realize the parallel acceleration calculation of the low-order nonlinear equations of the electromagnetic finite element model; The S3 includes the following steps: S301: Set the conductance and capacitance parameters as follows. The vector magnetic potential difference at 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 represents the conductivity of the component material; S302: Add transmission lines at both ends of the nonlinear element; S303: Perform Norton equivalent on the circuit containing the transmission line. Through the resistance equivalent capacitance with a time step Δt, the capacitance 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 shunt 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 shunt conductance of vertex M in the equivalent Norton circuit; G CN represents the shunt conductance of vertex N in the equivalent Norton circuit; S304: Obtain a current source through capacitance equivalence to complete the circuit equivalence of any triangular region Ω in the finite element model e : 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 of vertex K in the equivalent Norton circuit; I CM represents the current source current of vertex M in the equivalent Norton circuit; I CN represents the current source current of vertex N in the equivalent Norton circuit; Represents 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 formula (12), through the current source equivalent capacitance initial charge, set the characteristic admittance parameters of the transmission line as follows: 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 permeability. The closer the guessed permeability is to the true permeability, the faster the transmission line iterative convergence speed; S306: Assume that the voltage between nodes during the nth transmission line iteration is U t (n X ) t : U X (n t ) = U i (n t ) + U r (n t ) (22) U r (n t ) represents the reflected voltage during the n t -th transmission line iteration and is the input value; U i (n t ) represents the incident voltage at the n t -th transmission line iteration; S307: Calculate the incident voltage U during the (n + 1)-th transmission line iteration t (n + 1). i (n t + 1), Since the permeabilities v of all triangular regions share the same triangular region e , when calculating the incident voltage, the nonlinear resistors in all triangular regions must simultaneously satisfy Equation (23), and Equation (24) needs to be solved simultaneously: 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 are incident waves in the x, y, and z directions; V x0 , V y0 , V z0 are the reflected waves in the x, y, and z directions; S308: Solve formula (24) by the Newton 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 both ends of the node; In formula (26): I GKM (n t +1) represents the current of the current source between transmission line iteration nodes K and M in the (n t +1)-th transmission; I GMN (n t + 1) represents the current of the current source between transmission line iteration nodes M and N for the (n t + 1)-th time; I GNK (n t +1) represents the current of the current source between the n t +1-th transmission line iteration node N and node K; U KM (n t +1) represents the voltage difference between the n t +1-th transmission line iteration nodes K and M; U MN (n t + 1) represents the voltage difference between the transmission line iteration nodes M and N at the (n t + 1)-th time; U NK (n t +1) represents the voltage difference between the n t +1-th transmission line iteration node N and node K; S3010: Integrate the parameters in formulas (18) to (26) to achieve the purpose of solving formula (12). At this time, formula (12) becomes the linear equation formula (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 sum of the linear part and the characteristic admittance in Equation (21) for corresponding vertices; M ctlm Sum all the capacitive admittances in expression (19) according to the corresponding vertices; F tlm Indicates the sum of the current sources in the Norton equivalent circuit of the transmission line according to the corresponding vertices; F Alast The current sources generated by the equivalent capacitance in Expression (20) are accumulated and summed according to vertices; S4: Propose the relay improved proper orthogonal decomposition model reduction parallel finite element method to reduce the order of the high-order linear equations decoupled from the electromagnetic finite element model in S3, and realize the fast calculation of the relay electromagnetic finite element model; The S4 includes the following steps: S401: Based on the proper orthogonal decomposition theory, reduce the order of formula (27) and transform it into formula (28), In formula (28): K ctlmr Indicates that the characteristic admittance after order reduction is the sum of the corresponding vertices; A r represents the reduced vector magnetic potential; M ctlmr Indicates the sum of the reduced capacitance admittances corresponding to the vertices; F tlmr Indicates the sum of the current sources in the Norton equivalent circuit after order reduction according to the corresponding vertices; F Alastr Indicates that the current sources generated by the equivalent capacitance after order reduction are summed up by vertex; Ψ P represents an orthogonal matrix; denote the orthogonal matrix Ψ P is the transposed matrix of; F r represents the orthogonal matrix Ψ after rank reduction P and its transpose matrix is the product of the transpose matrix and the current density matrix F c ; S402: During each transmission line iteration, the non-reduced vector magnetic potential A is calculated from the reduced vector magnetic potential A according to Equation (29). r ​ 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 and the current sources generated by the equivalent capacitance are accumulated and summed at the vertices to obtain F Alast will change at each time step, so it is necessary to solve it again at each time step; The current source in the Norton equivalent circuit is according to the corresponding vertex and F tlm It will change after each transmission line iteration. Therefore, 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 it is guessed that the magnetic permeability does not change, the characteristic admittance is based on the corresponding vertices and K ctlm and the capacitance admittance is based on the corresponding vertices and M ctlm are all invariant parameters, and a single order reduction calculation can be performed.

2. A reduced-order parallel acceleration calculation method for a relay electromagnetic field finite element model according to claim 1, characterized in that: The S1 includes the following steps: S101: According to the working principle of the relay, use the finite element method to solve the vector magnetic potential A. The electromagnetic field equations of the relay are as follows: n·B = 0 (2) n×H = 0 (3) In formula (1): Denotes the symbol for partial derivative; v represents the magnetic permeability; A represents the vector magnetic potential; σ s represents conductivity; Js represents the externally applied current density; ▽ represents the vector differential operator, which is used to calculate the divergence of the magnetic induction intensity B; B r Represents the remanence of the permanent magnet; x is the unit vector, representing the x-axis direction; y is the unit vector, representing the y-axis direction; z is the unit vector, representing the z-axis direction; t represents time; Formulas (2) and (3) are the boundary conditions of the electromagnetic field of the relay, where: n is the unit vector, representing the normal direction of the interface, surface or boundary; B represents the magnetic induction intensity; H represents the magnetic field intensity; S102: According to the vector magnetic potential A, calculate the magnetic flux Ψ 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 for the electromagnetic force F1 on the armature within 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 for the dynamic characteristics of the armature according to the electromagnetic dynamic characteristic equations of the relay: Among them: The electromagnetic dynamic characteristic equations of the direct-acting relay are as follows: In formula (6): Ψ represents the magnetic flux linkage; U represents the coil excitation voltage; i represents the coil current; R represents the coil resistance; F l represents the electromagnetic force received by the armature; F f Indicates the reaction force of the contact spring system on the armature; m represents the mass of the armature; x1 represents the displacement of the armature; The electromagnetic dynamic characteristic equations of the rotating relay are as follows: In formula (7): T l Indicates the electromagnetic torque on the armature; T f Indicates the counter torque of the contact spring system on the armature; I x Indicates the moment of inertia of the armature; θ represents the rotation angle of the armature.

3. A reduced-order parallel acceleration calculation method for the electromagnetic field finite element model of a relay according to claim 2, characterized in that: The said S2 further includes the following steps: S206: After obtaining the magnetic vector potential solved according to Equation (12), the magnetic induction intensity B in the triangular region can be solved e : S207: Solve for the magnetic flux linkage according to formula (4) and formula (8): In formula (14): m coil represents the number of coil units; N c Indicates the number of turns of the coil; S coil represents the cross-sectional area of the coil; S208: Solve for the electromagnetic force on the relay armature according to formula (15): m ch represents the number of contour vertices of the moving part; S208: Simulate the electromagnetic performance of the relay based on the magnetic induction intensity, magnetic flux linkage, and electromagnetic force on the armature; S209: Since the vector magnetic potential and current density of the two-dimensional electromagnetic field only contain components in the z direction, therefore, formula (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 transformation, the vector magnetic potential and current density are perpendicular to the z-r plane. Then, considering the boundary conditions, we can obtain: In formula (17): A’ = rA represents the transformed vector magnetic potential; v’ = v / r represents the transformed magnetic permeability; z is a unit vector representing the z-axis direction; r is a unit vector representing the r-axis direction; Γ1 represents the first boundary condition of formula (17); Γ2 represents the second boundary condition of formula (17); S2011: Since after the axisymmetric transformation replaces the z-axis and r-axis of the cylindrical coordinate system with the x-axis and y-axis, formula (17) is the same as formula (16), so the finite element modeling and solution can be carried out according to formula (16).

Images