CN107609274B - Two-dimensional static magnetic field parallel finite element method based on transmission line and level scheduling method - Google Patents
Two-dimensional static magnetic field parallel finite element method based on transmission line and level scheduling method Download PDFInfo
- Publication number
- CN107609274B CN107609274B CN201710827716.2A CN201710827716A CN107609274B CN 107609274 B CN107609274 B CN 107609274B CN 201710827716 A CN201710827716 A CN 201710827716A CN 107609274 B CN107609274 B CN 107609274B
- Authority
- CN
- China
- Prior art keywords
- matrix
- node
- triangular
- nonlinear
- network
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Landscapes
- Complex Calculations (AREA)
Abstract
A two-dimensional static magnetic field parallel finite element method based on a transmission line and level scheduling method belongs to the field of electric appliance numerical calculation, and mainly solves a two-dimensional nonlinear static electromagnetic field, including two-dimensional plane and two-dimensional axial symmetry. The invention has the advantages that: the transmission line iteration method and the level scheduling method are adopted to carry out finite element iteration solution, the global matrix Y can be kept unchanged in the iteration solution process, the LU decomposition method is adopted in the matrix solution process, LU decomposition is only needed to be carried out in the first calculation step, and the LU decomposition generally occupies about 95% of the time of the matrix solution, so that by adopting the method, the LU decomposition process of the global matrix does not need to be carried out again in each iteration step, and 95% of the time can be saved. Meanwhile, the level scheduling method is applied to the matrix triangular solving process after LU decomposition, and the algorithm can effectively accelerate the triangular solving process.
Description
Technical Field
The invention belongs to the field of electric appliance numerical calculation, and particularly relates to a two-dimensional static magnetic field finite element parallel acceleration solving method combining a transmission line iteration method and a level scheduling method.
Background
The finite element method is the most common numerical calculation method in industrial design, is adopted by a plurality of commercial simulation software, and is widely applied. However, with the increasing complexity of the solution model and the increasing number of the sub-grid units, the nonlinear finite element solution method using the traditional newton iteration method as the core faces the problem of time-consuming solution, which is directly related to the speed and efficiency of product development.
The core of the finite element problem solution lies in solving a linear equation set, for a nonlinear problem, a new iteration result is used for regenerating a global matrix of a finite element model in each step of a traditional Newton iteration method, the dimension of the global matrix is continuously increased along with the continuous increase of the model sub-network, the time consumed by LU decomposition and the like of each step of the matrix is correspondingly increased, and the total solution time may be increased in a geometric form along with the densification of the sub-network.
Therefore, a new iterative method needs to be researched to solve the problems of long solving time and low efficiency caused by solving the finite element nonlinearity problem by the newton iterative method.
Disclosure of Invention
The invention aims to solve the problems of long solving time and low efficiency caused by solving the finite element nonlinear problem by the existing Newton iteration method, and provides a method for realizing finite element parallel accelerated solving on a two-dimensional nonlinear static magnetic field model by combining a transmission line iteration method and a level scheduling method.
In order to achieve the purpose, the technical scheme adopted by the invention is as follows:
a two-dimensional static magnetic field parallel finite element acceleration method based on a transmission line and level scheduling method comprises the following specific steps:
the method comprises the following steps: establishing a two-dimensional plane coordinate system and establishing a geometric model of the static magnetic field problem;
step two: for a control equation and boundary conditions in a two-dimensional nonlinear static magnetic field, a set of differential equations is obtained, wherein the control equation is as follows:
wherein A is the magnetic potential of the variable to be solved, mu0For air permeability, M is the magnetization vector, αmIs the included angle between the positive directions of the M and the x axis, and J is the current density; the boundary conditions are as follows:
Γ1：A＝0，Γ1representing the distribution of magnetic potential a over the boundary;
for a static magnetic field with a two-dimensional axisymmetric structure, a group of differential equations is obtained, and the control equation is as follows:
let μ ═ x μ, a ═ xA, to give
Wherein A 'is the magnetic potential of the variable to be solved, mu' is the air permeability, M is the vector of the magnetization, αmIs the included angle between the positive directions of the M and the x axis, and J is the current density; the boundary conditions are as follows:
Γ1：A'＝0，Γ1representing the distribution of magnetic potential a over the boundary;
the two-dimensional and two-dimensional axial symmetry problems can be solved by adopting the following steps,
step three: carrying out rough net division on a solving area, carrying out discrete net division on the solving area by using triangular units to obtain a finite element network comprising a plurality of triangular units, wherein the total number of the triangular units in the finite element network is N1, the total number of nodes is M1, and the triangular units and the nodes are respectively numbered from 1 to N1 and from 1 to M1;
step four: in each triangular cell, a cell coefficient matrix [ Y ] is calculatede]Matrix of excitation current source cells [ J ]e]Matrix of magnetization [ M ]e]Wherein [ Y ]e]Is a 3 × 3 matrix, [ Je]And [ Me]Are all 3 × 1 matrixes, K, L,N is the number of three nodes of the triangular unit according to the anticlockwise direction; wherein, in the two-dimensional case, the matrix [ Y ]e]Each of the elements of (1) is
(i, j ═ K, L, N), i, j are subscripts of variables, b, c are variable expressions, and the concrete form is
Delta is the area of the triangle, mu is the magnetic conductivity of the material, and x and y are the horizontal and vertical coordinate values of the point;
in the case of two-dimensional axial symmetry, the matrix [ Y ]e]Each of the elements of (1) is
(i, j ═ K, L, N), i, j are subscripts of variables, b, c are variable expressions, and the concrete form is
Delta is the area of the triangle, mu is the magnetic conductivity of the material, and x and y are the horizontal and vertical coordinate values of the point;
step five: according to the obtained unit coefficient matrix [ Y ] of each triangular unite]And a matrix of excitation current source cells [ Je]And a magnetization matrix [ M ]e]Carrying out finite element assembly on N1 triangular units to obtain global matrixes Y, J and M, wherein Y is an M1 multiplied by M1 matrix, and J, M is an M1 multiplied by 1 matrix;
step six: solving a nonlinear equation system YA ═ J + M according to the global matrix Y, J, M obtained in the step five to obtain the magnetic potential A of each node in the two-dimensional axisymmetric nonlinear static magnetic field, wherein A is a node magnetic potential matrix of M1 multiplied by 1, and A ═ A1A2LAM1]T；
Step seven: calculating the magnetic induction intensity B of each triangular unit according to the following formula according to the node magnetic potential matrix A obtained by calculation in the step six,
wherein, BxThe component of the magnetic induction B in the x-direction, ByComponent of magnetic induction B in the y-direction, AK、AL、ANRespectively are magnetic potential values of three points of delta KLN;
step eight: calculating the magnetic permeability mu of each triangular unit according to the B-H curve of the ferromagnetic material and the magnetic induction intensity B of each triangular unit calculated in the step seven;
step nine: carrying out fine triangulation on the solution domain based on the triangulation result in the third step to obtain a finite element network with the total number of triangular units being N2 and the total number of nodes being M2, and numbering the triangular units and the nodes from 1 to N2 and from 1 to M2 respectively;
step ten: according to the method in the fourth step, the finite element network obtained in the ninth step is subjected to the calculation again to obtain the element coefficient matrix [ Y ] of each triangular elemente]1And a matrix of excitation current source cells [ Je]1Magnetization matrix [ M ]e]1；
Step eleven: converting the finite element network into a circuit model, and obtaining a unit coefficient matrix [ Y ] in the step tene]1Viewed as an admittance matrix of the circuit, exciting a matrix of current source cells [ Je]1And a magnetization matrix [ M ]e]1Regarding as a current source matrix between each node and the ground, establishing an equivalent circuit network for each triangular unit in the finite element network;
step twelve: an assembly circuit, which is used for connecting the equivalent circuit network corresponding to each triangular unit established in the step eleven through nodes to assemble a complete nonlinear circuit network, wherein the nonlinear circuit network is equivalent to a circuit comprising a linear network and a plurality of nonlinear elements to be solved;
step thirteen: for the nonlinear circuit network obtained in the step twelve, in order to use a transmission line iteration method to solve, a section of transmission line needs to be added between the nonlinear element and the linear network, and due to the delay effect of the transmission line on signal transmission, the nonlinear solving process of the circuit comprises an incident stage and a reflection stage,
in the incident stage, the voltage signal of the nonlinear circuit element is incident to the linear network, which is equivalent to the parallel connection of the transmission line admittance and the virtual current source,
in the reflection stage, the voltage signal is transmitted to the nonlinear element by the linear network, the nonlinear element is solved, the incident stage and the reflection stage are iterated continuously in this way until the circuit reaches a steady state,
adding a section of transmission line between a linear part and a nonlinear element, wherein the admittance of the transmission line is calculated by the following method:
(1) determining the estimated value of the magnetic permeability mu of each triangular unit, checking the triangular unit of the first sub-network corresponding to the gravity center of the triangular unit obtained after the ninth sub-network, setting the magnetic permeability of the corresponding triangular unit of the first sub-network as the magnetic permeability of the triangular unit,
(2) the admittance of the nonlinear element is a variable related to the magnetic permeability mu, the mu value obtained in the last step is substituted into the nonlinear element expression, the obtained result is used as the admittance value of the corresponding transmission line,
(II) when the iteration starts, the voltage of each node is 0, and when the voltage signal of the nth node is the incident voltage VinWhen reflecting to the linear network, each nonlinear element to be solved is equivalent to a parallel circuit of an admittance and a current source, wherein the admittance is the corresponding transmission line admittance YnThe current value in the current source is 2VinYnSolving the equivalent circuit, wherein the solved expression is YV (J + M + I)c+2VinYnIn which IcFor the voltage-controlled current source in the circuit, Y is the admittance matrix of the circuit, and Y is kept unchanged in the iterative process, when the voltage value V of each node is obtainedin，YV＝J+M+Ic+2VinYnThe matrix Y is kept unchanged in the iteration process, so that the matrix Y is decomposed into LUs in the first step of the iteration, namely Y is equal to LU, LUs are respectively a lower triangular matrix and an upper triangular matrix, and then the LUV is directly triangulated to J + M + I in each solvingc+2VinYnSolving is carried out by adopting a parallel level scheduling method;
thirdly, calculating and updating the admittance value of the nonlinear element according to the voltage value of each node by using a relational expression between the nonlinear element and the voltage,
(IV) calculating the voltage value V of each node incident to the nonlinear elementrnV at node nrn＝Vn-Vin，
(V) incident process, each nonlinear element to be solved is equivalent to a parallel circuit of an admittance and a current source, wherein the admittance is the corresponding transmission line admittance YnThe current value in the current source is 2VrnYnObtaining the voltage V across each nonlinear elementNEnSince the calculation of each element can be calculated independently, the step adopts a parallel algorithm to calculate,
sixthly, calculating the voltage value V of each node incident to the linear networkinV at node nin＝VNEn-Vrn，
And (seventhly) repeating the steps from the second step to the sixth step until the voltage value V obtained in the step (two) in two adjacent iterationsnWhen the preset convergence error is reached, the calculated voltage value V of each node is obtainednThe voltage value is the calculated voltage value;
fourteen steps: and drawing a magnetic potential cloud picture in the two-dimensional nonlinear static magnetic field according to the voltage value of each node.
Compared with the prior art, the invention has the beneficial effects that: the transmission line iteration method and the level scheduling method are adopted to carry out finite element iteration solution, the global matrix Y can be kept unchanged in the iteration solution process, the LU decomposition method is adopted in the matrix solution process, LU decomposition is only needed to be carried out in the first calculation step, and the LU decomposition generally occupies about 95% of the time of the matrix solution, so that by adopting the method, the LU decomposition process of the global matrix does not need to be carried out again in each iteration step, and 95% of the time can be saved. Meanwhile, the level scheduling method is applied to the matrix triangular solving process after LU decomposition, and the algorithm can effectively accelerate the triangular solving process.
Drawings
FIG. 1 is a lower triangular matrix and level diagram;
FIG. 2 is a diagram of a lower triangular matrix sorted by level;
FIG. 3 is a block diagram of a lower triangular matrix of 4096 × 4096 in rank order;
FIG. 4 is a static magnetic field geometric model diagram;
FIG. 5 is a graph of the calculation time of each iteration step using Newton's iteration to calculate the number of CPU cores and the number of sub-grid units;
FIG. 6 is a graph of the calculation time of each iteration step for calculating the number of CPU cores and the number of sub-grid units by using the method of the present invention;
FIG. 7 is a graph of total computation time for computing the number of CPU cores and the number of grid elements by Newton's iteration;
FIG. 8 is a graph of the total computation time for the number of CPU cores and the number of the sub-grid units.
Detailed Description
The technical solution of the present invention is further described below with reference to the accompanying drawings, but not limited thereto, and any modification or equivalent replacement of the technical solution of the present invention without departing from the spirit and scope of the technical solution of the present invention shall be covered by the protection scope of the present invention.
The first embodiment is as follows: the embodiment describes a two-dimensional static magnetic field parallel finite element acceleration method based on a transmission line and level scheduling method, which comprises the following specific steps:
the method comprises the following steps: establishing a two-dimensional plane coordinate system and establishing a geometric model of the static magnetic field problem, as shown in FIG. 4;
step two: for a control equation and boundary conditions in a two-dimensional nonlinear static magnetic field, a set of differential equations is obtained, wherein the control equation is as follows:
wherein A is the magnetic potential of the variable to be solved, mu0For air permeability, M is the magnetization vector, αmIs the included angle between the positive directions of the M and the x axis, and J is the current density; the boundary conditions are as follows:
Γ1：A＝0，Γ1representing the distribution of magnetic potential a over the boundary;
for a static magnetic field with a two-dimensional axisymmetric structure, a group of differential equations is obtained, and the control equation is as follows:
let μ ═ x μ, a ═ xA, to give
Wherein A 'is the magnetic potential of the variable to be solved, mu' is the air permeability, M is the vector of the magnetization, αmIs the included angle between the positive directions of the M and the x axis, and J is the current density; the boundary conditions are as follows:
Γ1：A'＝0，Γ1representing the distribution of magnetic potential a over the boundary;
the two-dimensional and two-dimensional axial symmetry problems can be solved by adopting the following steps,
step three: carrying out rough net division on a solving area, carrying out discrete net division on the solving area by using triangular units to obtain a finite element network comprising a plurality of triangular units, wherein the total number of the triangular units in the finite element network is N1, the total number of nodes is M1, and the triangular units and the nodes are respectively numbered from 1 to N1 and from 1 to M1;
step four: in each triangular cell, a cell coefficient matrix [ Y ] is calculatede]Matrix of excitation current source cells [ J ]e]Matrix of magnetization [ M ]e]Wherein [ Y ]e]Is a 3 × 3 matrix, [ Je]And [ Me]All the three nodes are 3 multiplied by 1 matrixes, and K, L, N are numbered as three nodes of a triangular unit according to the anticlockwise direction; wherein, in the two-dimensional case, the matrix [ Y ]e]Each of the elements of (1) is
(i, j ═ K, L, N), i, j are subscripts of variables, b, c are variable expressions, and the concrete form is
Delta is the area of the triangle, mu is the magnetic conductivity of the material, and x and y are the horizontal and vertical coordinate values of the point;
in the case of two-dimensional axial symmetry, the matrix [ Y ]e]Each of the elements of (1) is
(i, j ═ K, L, N), i, j are subscripts of variables, b, c are variable expressions, and the concrete form is
Delta is the area of the triangle, mu is the magnetic conductivity of the material, and x and y are the horizontal and vertical coordinate values of the point;
step five: according to the obtained unit coefficient matrix [ Y ] of each triangular unite]And a matrix of excitation current source cells [ Je]And a magnetization matrix [ M ]e]Carrying out finite element assembly on N1 triangular units to obtain global matrixes Y, J and M, wherein Y is an M1 multiplied by M1 matrix, and J, M is an M1 multiplied by 1 matrix;
step six: solving a nonlinear equation system YA ═ J + M according to the global matrix Y, J, M obtained in the step five to obtain the magnetic potential A of each node in the two-dimensional axisymmetric nonlinear static magnetic field, wherein A is a node magnetic potential matrix of M1 multiplied by 1, and A ═ A1A2LAM1]T；
Step seven: calculating the magnetic induction intensity B of each triangular unit according to the following formula according to the node magnetic potential matrix A obtained by calculation in the step six,
wherein, BxThe component of the magnetic induction B in the x-direction, ByComponent of magnetic induction B in the y-direction, AK、AL、ANRespectively are magnetic potential values of three points of delta KLN;
step eight: calculating the magnetic permeability mu of each triangular unit according to the B-H curve of the ferromagnetic material and the magnetic induction intensity B of each triangular unit calculated in the step seven;
step nine: carrying out fine triangulation on the solution domain based on the triangulation result in the third step to obtain a finite element network with the total number of triangular units being N2 and the total number of nodes being M2, and numbering the triangular units and the nodes from 1 to N2 and from 1 to M2 respectively;
step ten: according to the method in the fourth step, the finite element network obtained in the ninth step is subjected to the calculation again to obtain the element coefficient matrix [ Y ] of each triangular elemente]1And a matrix of excitation current source cells [ Je]1Magnetization matrix [ M ]e]1；
Step eleven: converting the finite element network into a circuit model, and obtaining a unit coefficient matrix [ Y ] in the step tene]1Viewed as an admittance matrix of the circuit, exciting a matrix of current source cells [ Je]1And a magnetization matrix [ M ]e]1Regarding as a current source matrix between each node and the ground, establishing an equivalent circuit network for each triangular unit in the finite element network;
step twelve: an assembly circuit, which is used for connecting the equivalent circuit network corresponding to each triangular unit established in the step eleven through nodes to assemble a complete nonlinear circuit network, wherein the nonlinear circuit network is equivalent to a circuit comprising a linear network and a plurality of nonlinear elements to be solved;
step thirteen: for the nonlinear circuit network obtained in the step twelve, in order to use a transmission line iteration method to solve, a section of transmission line needs to be added between the nonlinear element and the linear network, and due to the delay effect of the transmission line on signal transmission, the nonlinear solving process of the circuit comprises an incident stage and a reflection stage,
in the incident stage, the voltage signal of the nonlinear circuit element is incident to the linear network, which is equivalent to the parallel connection of the transmission line admittance and the virtual current source,
in the reflection stage, the voltage signal is transmitted to the nonlinear element by the linear network, the nonlinear element is solved, the incident stage and the reflection stage are iterated continuously in this way until the circuit reaches a steady state,
adding a section of transmission line between a linear part and a nonlinear element, wherein the admittance of the transmission line is calculated by the following method:
(1) determining the estimated value of the magnetic permeability mu of each triangular unit, checking the triangular unit of the first sub-network corresponding to the gravity center of the triangular unit obtained after the ninth sub-network, setting the magnetic permeability of the corresponding triangular unit of the first sub-network as the magnetic permeability of the triangular unit,
(2) the admittance of the nonlinear element is a variable related to the magnetic permeability mu, the mu value obtained in the last step is substituted into the nonlinear element expression, the obtained result is used as the admittance value of the corresponding transmission line,
(II) when the iteration starts, the voltage of each node is 0, and when the voltage signal of the nth node is the incident voltage VinWhen reflecting to the linear network, each nonlinear element to be solved is equivalent to a parallel circuit of an admittance and a current source, wherein the admittance is the corresponding transmission line admittance YnThe current value in the current source is 2VinYnSolving the equivalent circuit, wherein the solved expression is YV (J + M + I)c+2VinYnIn which IcFor the voltage-controlled current source in the circuit, Y is the admittance matrix of the circuit, and Y is kept unchanged in the iterative process, when the voltage value V of each node is obtainedin，YV＝J+M+Ic+2VinYnThe matrix Y is kept unchanged in the iteration process, so that the matrix Y is decomposed into LUs in the first step of the iteration, namely Y is equal to LU, LUs are respectively a lower triangular matrix and an upper triangular matrix, and then the LUV is directly triangulated to J + M + I in each solvingc+2VinYnSolving is carried out by adopting a parallel level scheduling method;
thirdly, calculating and updating the admittance value of the nonlinear element according to the voltage value of each node by using a relational expression between the nonlinear element and the voltage,
(IV) calculating the voltage value V of each node incident to the nonlinear elementrnV at node nrn＝Vn-Vin，
(V) incident process, each nonlinear element to be solved is equivalent to a parallel circuit of an admittance and a current source, wherein the admittance is the corresponding transmission line admittance YnThe current value in the current source is 2VrnYnObtaining the voltage V across each nonlinear elementNEnDue to the fact thatThe calculations of the various elements can be calculated separately, this step being calculated using parallel algorithms,
sixthly, calculating the voltage value V of each node incident to the linear networkinV at node nin＝VNEn-Vrn，
And (seventhly) repeating the steps from the second step to the sixth step until the voltage value V obtained in the step (two) in two adjacent iterationsnWhen the preset convergence error is reached, the calculated voltage value V of each node is obtainednThe voltage value is the calculated voltage value;
fourteen steps: and drawing a magnetic potential cloud picture in the two-dimensional nonlinear static magnetic field according to the voltage value of each node.
The second embodiment is as follows: in a two-dimensional static magnetic field parallel finite element method based on the transmission line and level scheduling method according to the first embodiment, in step eleven, an equivalent circuit network is established as follows:
will the cell matrix [ Ye]The elements on the diagonal are considered to be self-conductive, the elements on the off-diagonal are considered to be mutually conductive,
for elements on the off-diagonal, if
Representation matrix [ Ye]And the elements of the r-th row and the s-th column are provided with a voltage-controlled current source between a node r and a node s in the equivalent circuit network corresponding to the triangular unit, and the current in the controlled current source is UrsYrsThe direction is from node r to node s, where UrsIs the magnetic potential difference between node r and node s,
for elements on the off-diagonal, if
A pure resistor is arranged between the node r and the node s in the equivalent circuit network corresponding to the triangular unit, and the admittance of the pure resistor is
A current source is arranged between each node and the ground, and the current magnitude in the current source between the node K, the node L, the node N and the ground is JK、JL、JNThe direction of current flow is from ground to the node.
The third concrete implementation mode: in a two-dimensional static magnetic field parallel finite element method based on the transmission line and the level scheduling method according to the first embodiment, in step thirteen (second), the level scheduling method is established as follows:
taking the solution of the following triangular matrix Lx ═ b as an example, L is a lower triangular matrix with the size of n × n, x and b are both n × 1 matrices, and the calculation steps are as follows:
A. calculating the level (i) of the ith variable x (i) in the matrix L, and if the i-th row of the matrix L contains all the elements L in the i-th (i-1, 2, …, n) th column, the j-th (j-1, 2, …, n) th row of the matrix LijNot zero, then level (i) ═ 1+ maxlevel (j) is updated, until all elements are traversed,
B. after the level of each variable is obtained through calculation, performing row-column transformation on the matrix L according to the level of the variable x, so that the levels of the variables are in a sequence from small to large, the newly sequenced matrix is changed into L ', x and b are also transformed, and L' x 'is obtained as b', as shown in FIGS. 1 and 2, the matrix is sequenced according to the level, and the number in a small box is the level of the variable x;
C. as shown in fig. 3, a lower triangular matrix with a size of 4096 × 4096 is a matrix sorted according to the level size, wherein the small black dots represent that the position element is nonzero, variables with different levels are separated by dotted lines, the variables at the same level can be solved simultaneously, in level 1, the values of the variables are calculated by x ' (j) ═ b ' (j)/L ' (j), in the same level, the solution of each variable is put into different calculation units to realize parallel calculation,
the variable values in the rectangular area under the square area on the diagonal in level 1 have been solved in level 1, in the equation L ' x ' b ', the solved x ' is substituted, then the constant is moved to the right, b ' is updated,
as in the above calculation method, the variable value in the square area in level 2 is calculated next, the variable value is solved using the formula x '(j) ═ b' (j)/L '(j), the calculated x' is substituted in the equation L 'x' ═ b ', then the constant is shifted to the right, b' is updated again,
the subsequent level calculation is the same as the above, until the calculation is finished, the solution of the unknown variable x is completed;
the solving method of the upper triangular matrix Ux ═ b is similar to the solving process of L, but the solving process is performed from the opposite direction.
From the comparison of time in fig. 5, 6, 7 and 8, it is found that the calculation time used in the present invention is much more advantageous than the newton iteration method used in the conventional finite element software.
Claims (2)
1. A two-dimensional static magnetic field parallel finite element method based on a transmission line and level scheduling method is characterized in that: the method comprises the following specific steps:
the method comprises the following steps: establishing a two-dimensional plane coordinate system and establishing a geometric model of the static magnetic field problem;
step two: for a control equation and boundary conditions in a two-dimensional nonlinear static magnetic field, a set of differential equations is obtained, wherein the control equation is as follows:
wherein A is the magnetic potential of the variable to be solved, mu0For air permeability, M is the magnetization vector, αmIs the included angle between the positive directions of the M and the x axis, and J is the current density; the boundary conditions are as follows:
Γ1：A＝0，Γ1representing the distribution of magnetic potential a over the boundary;
Γ2：
Γ2represents the rate of change of the magnetic potential a along the outer normal direction of the boundary;
for a static magnetic field with a two-dimensional axisymmetric structure, a group of differential equations is obtained, and the control equation is as follows:
let mu' x mu0A' ═ xA, to give
Wherein A 'is the magnetic potential of the variable to be solved, mu' is the air permeability, M is the vector of the magnetization, αmIs the included angle between the positive directions of the M and the x axis, and J is the current density; the boundary conditions are as follows:
Γ1：A'＝0，Γ1represents the distribution of the magnetic potential A' on the boundary;
Γ2：
Γ2represents the rate of change of the magnetic potential a' along the outer normal direction of the boundary;
the two-dimensional and two-dimensional axial symmetry problems can be solved by adopting the following steps,
step three: carrying out rough net division on a solving area, carrying out discrete net division on the solving area by using triangular units to obtain a finite element network comprising a plurality of triangular units, wherein the total number of the triangular units in the finite element network is N1, the total number of nodes is M1, and the triangular units and the nodes are respectively numbered from 1 to N1 and from 1 to M1;
step four: in each triangular cell, a cell coefficient matrix [ Y ] is calculatede]Matrix of excitation current source cells [ J ]e]Matrix of magnetization [ M ]e]Wherein [ Y ]e]Is a 3 × 3 matrix, [ Je]And [ Me]All the three nodes are 3 multiplied by 1 matrixes, and K, L, N are numbered as three nodes of a triangular unit according to the anticlockwise direction; wherein, in the two-dimensional case, the matrix [ Y ]e]Each of the elements of (1) is
i ═ K, L, N and j ═ K, L, N, i, j are subscripts of variables, b, c are variable expressions of specific form
Delta is the area of the triangle, mu is the magnetic conductivity of the material, and x and y are the horizontal and vertical coordinate values of the point;
in the case of two-dimensional axial symmetry, the matrix [ Y ]e]Each of the elements of (1) is
i ═ K, L, N and j ═ K, L, N, i, j are subscripts of variables, b, c are variable expressions of specific form
Delta is the area of the triangle, mu is the magnetic conductivity of the material, and x and y are the horizontal and vertical coordinate values of the point;
step five: according to the obtained unit coefficient matrix [ Y ] of each triangular unite]And a matrix of excitation current source cells [ Je]And a magnetization matrix [ M ]e]Carrying out finite element assembly on N1 triangular units to obtain global matrixes Y0, J0 and M0, wherein Y0 is an M1 multiplied by M1 matrix, and J0 and M0 are M1 multiplied by 1 matrices;
step six: solving a nonlinear equation system Y0A ═ J0+ M0 according to the global matrixes Y0, J0 and M0 obtained in the step five to obtain a magnetic potential A in the two-dimensional axisymmetric nonlinear static magnetic field, wherein A is a node magnetic potential matrix of M1 multiplied by 1, and A ═ is1A2…AM1]T，A1，A2，AM1The magnetic potential values of the 1 st, 2 nd and M1 th nodes are respectively;
step seven: calculating the magnetic induction intensity B of each triangular unit according to the following formula according to the node magnetic potential matrix A obtained by calculation in the step six,
wherein, BxThe component of the magnetic induction B in the x-direction, ByComponent of magnetic induction B in the y-direction, AK、AL、ANRespectively are magnetic potential values of three points of delta KLN;
step eight: calculating the magnetic permeability mu of each triangular unit according to the B-H curve of the ferromagnetic material and the magnetic induction intensity B of each triangular unit calculated in the step seven;
step nine: carrying out fine triangulation on the solution domain based on the triangulation result in the third step to obtain a finite element network with the total number of triangular units being N2 and the total number of nodes being M2, and numbering the triangular units and the nodes from 1 to N2 and from 1 to M2 respectively;
step ten: according to the method in the fourth step, the finite element network obtained in the ninth step is subjected to the calculation again to obtain the element coefficient matrix [ Y ] of each triangular elemente]1And a matrix of excitation current source cells [ Je]1Magnetization matrix [ M ]e]1；
Step eleven: converting the finite element network into a circuit model, and obtaining a unit coefficient matrix [ Y ] in the step tene]1Viewed as an admittance matrix of the circuit, exciting a matrix of current source cells [ Je]1And a magnetization matrix [ M ]e]1Regarding as a current source matrix between each node and the ground, establishing an equivalent circuit network for each triangular unit in the finite element network;
step twelve: an assembly circuit, which is used for connecting the equivalent circuit network corresponding to each triangular unit established in the step eleven through nodes to assemble a complete nonlinear circuit network, wherein the nonlinear circuit network is equivalent to a circuit comprising a linear network and a plurality of nonlinear elements to be solved;
step thirteen: for the nonlinear circuit network obtained in the step twelve, in order to use a transmission line iteration method to solve, a section of transmission line needs to be added between the nonlinear element and the linear network, and due to the delay effect of the transmission line on signal transmission, the nonlinear solving process of the circuit comprises an incident stage and a reflection stage,
in the incident stage, the voltage signal of the nonlinear circuit element is incident to the linear network, which is equivalent to the parallel connection of the transmission line admittance and the virtual current source,
in the reflection stage, the voltage signal is transmitted to the nonlinear element by the linear network, the nonlinear element is solved, the incident stage and the reflection stage are iterated continuously in this way until the circuit reaches a steady state,
adding a section of transmission line between a linear part and a nonlinear element, wherein the admittance of the transmission line is calculated by the following method:
(1) determining the estimated value of the magnetic permeability mu of each triangular unit, checking the triangular unit of the first sub-network corresponding to the gravity center of the triangular unit obtained after the ninth sub-network, setting the magnetic permeability of the corresponding triangular unit of the first sub-network as the magnetic permeability of the triangular unit,
(2) the admittance of the nonlinear element is a variable related to the magnetic permeability mu, the mu value obtained in the last step is substituted into the nonlinear element expression, the obtained result is used as the admittance value of the corresponding transmission line,
(II) when the iteration starts, the voltage of each node is 0, and when the voltage signal of the nth node is the incident voltage VinWhen reflecting to the linear network, each nonlinear element to be solved is equivalent to a parallel circuit of an admittance and a current source, wherein the admittance is the corresponding transmission line admittance YnThe current value in the current source is 2VinYnTo the equivalent circuitSolving the expression YV as J + M + Ic+2VinYnIn which IcFor the voltage-controlled current source in the circuit, Y is the admittance matrix of the circuit, and Y is kept unchanged in the iterative process, when the voltage value V of each node is obtainedin，YV＝J+M+Ic+2VinYnThe matrix Y is kept unchanged in the iteration process, so that the matrix Y is decomposed into LUs in the first step of the iteration, namely Y is equal to LU, LUs are respectively a lower triangular matrix and an upper triangular matrix, and then the LUV is directly triangulated to J + M + I in each solvingc+2VinYnSolving is carried out by adopting a parallel level scheduling method, and V is the node voltage in the circuit;
thirdly, calculating and updating the admittance value of the nonlinear element according to the voltage value of each node by using a relational expression between the nonlinear element and the voltage,
(IV) calculating the voltage value V of each node incident to the nonlinear elementrnV at node nrn＝Vn-Vin，VnIs the voltage value of the nth node;
(V) incident process, each nonlinear element to be solved is equivalent to a parallel circuit of an admittance and a current source, wherein the admittance is the corresponding transmission line admittance YnThe current value in the current source is 2VrnYnObtaining the voltage V across each nonlinear elementNEnSince the calculation of each element can be calculated independently, the step adopts a parallel algorithm to calculate,
sixthly, calculating the voltage value V of each node incident to the linear networkinV at node nin＝VNEn-Vrn，
And (seventhly) repeating the steps from the second step to the sixth step until the voltage value V obtained in the step (two) in two adjacent iterationsnWhen the preset convergence error is reached, the calculated voltage value V of each node is obtainednThe voltage value is the calculated voltage value;
fourteen steps: and drawing a magnetic potential cloud picture in the two-dimensional nonlinear static magnetic field according to the voltage value of each node.
2. A two-dimensional static magnetic field parallel finite element method based on transmission line and level scheduling method as claimed in claim 1, wherein: in the eleventh step, the method for establishing the equivalent circuit network comprises the following steps:
will the cell matrix [ Ye]The elements on the diagonal are considered to be self-conductive, the elements on the off-diagonal are considered to be mutually conductive,
for elements on the off-diagonal, if
Representation matrix [ Ye]And the elements of the r-th row and the s-th column are provided with a voltage-controlled current source between the node r and the node s in the equivalent circuit network corresponding to the triangular unit, and the current in the voltage-controlled current source is UrsYrsThe direction is from node r to node s, where UrsIs the magnetic potential difference between node r and node s,
for elements on the off-diagonal, if
A pure resistor is arranged between the node r and the node s in the equivalent circuit network corresponding to the triangular unit, and the admittance of the pure resistor is
A current source is arranged between each node and the ground, and the current magnitude in the current source between the node K, the node L, the node N and the ground is JK、JL、JNThe direction of current flow is from ground to the node.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710827716.2A CN107609274B (en) | 2017-09-14 | 2017-09-14 | Two-dimensional static magnetic field parallel finite element method based on transmission line and level scheduling method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710827716.2A CN107609274B (en) | 2017-09-14 | 2017-09-14 | Two-dimensional static magnetic field parallel finite element method based on transmission line and level scheduling method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107609274A CN107609274A (en) | 2018-01-19 |
| CN107609274B true CN107609274B (en) | 2020-03-13 |
Family
ID=61063873
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710827716.2A Active CN107609274B (en) | 2017-09-14 | 2017-09-14 | Two-dimensional static magnetic field parallel finite element method based on transmission line and level scheduling method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107609274B (en) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109408927B (en) * | 2018-10-13 | 2022-01-18 | 哈尔滨工业大学 | Two-dimensional static magnetic field parallel finite element acceleration method based on black box transmission line model |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106227982A (en) * | 2016-08-27 | 2016-12-14 | 国网冀北电力有限公司电力科学研究院 | A kind of electromagnetic relay static characteristic computational methods and device |
| CN106294894A (en) * | 2015-05-15 | 2017-01-04 | 南京理工大学 | Quickly analyze the finite element boundary integral method of non-homogeneous electromagnetic characteristic of scattering |
| CN106649939A (en) * | 2016-09-28 | 2017-05-10 | 哈尔滨工业大学 | Transmission line iteration-based solving method for 2D axial symmetric nonlinear magnetostatic field model |
-
2017
- 2017-09-14 CN CN201710827716.2A patent/CN107609274B/en active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106294894A (en) * | 2015-05-15 | 2017-01-04 | 南京理工大学 | Quickly analyze the finite element boundary integral method of non-homogeneous electromagnetic characteristic of scattering |
| CN106227982A (en) * | 2016-08-27 | 2016-12-14 | 国网冀北电力有限公司电力科学研究院 | A kind of electromagnetic relay static characteristic computational methods and device |
| CN106649939A (en) * | 2016-09-28 | 2017-05-10 | 哈尔滨工业大学 | Transmission line iteration-based solving method for 2D axial symmetric nonlinear magnetostatic field model |
Non-Patent Citations (2)
| Title |
|---|
| 基于云计算的电磁问题并行计算方法;金亮 等;《电工技术学报》;20161130;第31卷(第22期);第5-11页 * |
| 电磁场计算中的四边形组合有限元法;路勇 等;《哈尔滨电工学院学报》;19950630;第18卷(第2期);第181-185页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107609274A (en) | 2018-01-19 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Dulmage et al. | Two algorithms for bipartite graphs | |
| Einstein et al. | The gravitational equations and the problem of motion | |
| Taylor et al. | An algorithm for computing Fekete points in the triangle | |
| JP5426716B2 (en) | Distribution network power flow analysis system and method | |
| CN106646645B (en) | A kind of gravity forward modeling accelerated method | |
| CN110457790A (en) | The discontinuous golden finite element method of gal the Liao Dynasty of near field dynamics for malformation analysis | |
| CN111062610B (en) | Power system state estimation method and system based on information matrix sparse solution | |
| CN111709150B (en) | Simulation method for magnetic field spatial distribution of magnet in any shape | |
| CN107069696A (en) | A kind of parallel calculating method of Power system state estimation | |
| CN107038308B (en) | A kind of regular grid terrain modeling method based on linear interpolation | |
| CN103632046A (en) | Power grid load flow calculation method | |
| CN107609274B (en) | Two-dimensional static magnetic field parallel finite element method based on transmission line and level scheduling method | |
| CN104809258B (en) | The adaptive amending method of dough sheet normal vector in electromagnetic scattering simulation modeling | |
| CN111339688A (en) | Method for solving rocket simulation model time domain equation based on big data parallel algorithm | |
| Kulikov et al. | The Principles of Discrete Modeling of Rod Constructions of Architectural Objects | |
| CN105205299B (en) | The quick Dimension Reduction Analysis method of TV university electromagnetic characteristic of scattering | |
| CN106649939A (en) | Transmission line iteration-based solving method for 2D axial symmetric nonlinear magnetostatic field model | |
| CN107239586A (en) | The domain decomposition parallel method effective to unconditionally stable FDTD method | |
| Dumont et al. | The best of two worlds: The expedite boundary element method | |
| CN115329975A (en) | Simulation method, device, equipment and storage medium | |
| CN109408927B (en) | Two-dimensional static magnetic field parallel finite element acceleration method based on black box transmission line model | |
| CN112685936B (en) | Modeling method for shell mother-of-pearl microstructure finite element analysis | |
| CN115310227A (en) | Method for constructing gear power propagation digital twin model | |
| CN114741956A (en) | Carbon emission flow tracking method based on graph neural network | |
| CN109038588A (en) | A kind of distribution power system load flow calculation method |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant |