CC BY
45
Electrical Machines and Apparatus
UDC 621.316+621.316.53 doi: 10.20998/2074-272X.2020.5.01
E.I. Baida, |B.V. Klymenko|, M.G. Pantelyat, Yu.A. Yelanskyi, D. Trichet, G. Wasselynck
CHALLENGES OF DYNAMIC SIMULATION OF HIGH-SPEED ELECTROMAGNETIC VALVES OF GAS DISTRIBUTION DEVICES
High-speed electromagnetic valves of gas distribution devices are used in modern missile and space technology as jet micromotors of the executive elements of missile stabilization systems, as well as to control the movement of spacecrafts in space. The problem of creating such valves which are simple and reliable in the operation is relevant. In this work, it is proposed at the development and design stage to perform computer modelling of mutually coupled electromechanical processes, such as: distribution of transient electromagnetic field, transients in an electric circuit, and movement of an electromagnet armature. Besides, the calculation of the force with which the compressed gas acts on the corresponding structural elements of the valve is proposed to be performed by solving the system of Navier-Stokes equations. All problems are solved by numerical methods in axisymmetrical formulation with the corresponding initial and boundary conditions. Improvement of the accuracy of electromagnetic calculations and taking into account the movement of the armature of an electromagnet in the process of multiphysics numerical simulation is achieved using so-called tunable elastic meshes. The paper presents a comparative analysis of the numerical results obtained for several designs of electromagnets. The features of the dynamics of high-speed electromagnets of gas distribution valves during on and off operations are analyzed, the corresponding dynamic characteristics calculated using the proposed technique are presented. References 17, table 1, figures 10. Key words: high-speed electromagnets, dynamics, Finite Element Method, multiphysics, elastic mesh.
Быстродействующие электромагнитные клапаны газораспределительных устройств применяются в современной ракетно-космической технике в качестве реактивных микродвигателей исполнительных органов систем стабилизации ракет, а также для управления движением космическими летательными аппаратами в пространстве. Проблема создания простых в эксплуатации и надежных в работе указанных клапанов является актуальной. В настоящей работе предлагается на стадии разработки и проектирования выполнять компьютерное моделирование взаимосвязанных электромеханических процессов, таких как: распределение нестационарного электромагнитного поля, переходные процессы в электрической цепи, движение якоря электромагнита. При этом расчет силы, с которой сжатый газ действует на соответствующие конструктивные элементы клапана, предлагается выполнять путем решения системы уравнений Навье-Стокса Все задачи решаются численными методами в осесимметричной постановке с соответствующими начальными и граничными условиями. Повышение точности электромагнитных расчетов и учет движения якоря электромагнита в процессе мультифизического численного моделирования достигается благодаря использованию так называемых перестраиваемых упругих сеток. В работе приведен сравнительный анализ численных результатов, полученных для нескольких конструкций электромагнитов. Проанализированы особенности динамики быстродействующих электромагнитов газораспределительных клапанов при выполнении операций включения и отключения, приведены соответствующие динамические характеристики, рассчитанные по предложенной методике. Библ. 17, табл. 1, рис. 10. Ключевые слова: быстродействующие электромагниты, динамика, метод конечных элементов, мультифизика, упругая сетка.
Introduction. The use of jet micro-motors as the executive elements of missile stabilization systems [1] imposes stringent requirements on their dynamic characteristics. The range of control forces of such motors is much smaller than the forces created by missile engines designed to bring the spacecraft to a given flight path. To ensure the operation of micro-motors, a special power system and automation units (valves, regulators, starting devices) are required. Micro-motors with a power system form an autonomous system to obtain a control force due to the reaction of gas outflow from the nozzles of the motors.
Among the main problems of improving the gas-jet executive elements described above, it is necessary to note the problem of creating easy-to-operate and reliable in operation high-speed electromagnetic valves of gas distribution devices. As mentioned above, such valves are used, in particular, in modern missile and space technology in the system of executive devices to control the movement of spacecrafts in space (orientation, stabilization, orbit correction, maneuver, etc.) [1, 2].
In the vast majority of DC electromagnets, including the electromagnets of the drive of high-speed gas distribution valves, shunt windings [3] are used, which are
connected in parallel to the control circuit power sources.
At the stage of development and design of objects with DC electromagnets, including those used in space and missile technology, it is necessary to carry out numerical solution of the multiphysics problem including computer simulation of the following coupled phenomena and processes: distribution of transient electromagnetic field, transient current changes in the winding, and the movement of the armature of an electromagnet. Here, numerical solution of the problem of calculation of the force with which compressed gas acts on the corresponding structural elements of the valve should be obtained by solving the system of Navier-Stokes equations for the valve model. Despite the large number of publications on the calculation of DC electromagnets, there are no publications on complex computer modelling of these interrelated multiphysical processes.
A feature of the dynamics of electromagnets used in micro-motors of control systems for orienting spacecrafts is a small (about 0.3 mm) change in the size of the working gap at the displacement of its movable system, high speed of the order of 10-15 ms and short travel time of the order of 1-1.5 ms. Therefore, improvement of the
© E.I. Baida, |B.V. Klymenko|, M.G. Pantelyat, Yu.A. Yelanskyi, D. Trichet, G. Wasselynck
accuracy of electromagnetic calculations and reducing the time required to perform calculations in the process of multiphysical modelling in such cases can be achieved by the use of so-called tunable elastic meshes. This should also be implemented when developing an algorithm for numerical solution of the problem.
The goal of the paper is to develop a mathematical model of coupled multiphysics processes in high-speed DC electromagnets, to construct a numerical algorithm for calculation of the dynamics of high-speed DC electromagnets based on the technique of tunable elastic meshes as well as to apply the developed mathematical model and algorithm for numerical analysis of the dynamics of high-speed electromagnets of gas distribution valves of micro-motors of missile and space technology.
Problem definition. This paper presents a technique and results of the calculation of the dynamics of a highspeed armored-type electromagnet used in missile and space technology, in particular, in micro-motors of control systems for the orientation of rockets [3, 4].
Figure 1 shows a sketch of one of the designs of an electromagnet gas distribution valve with a switching ball as a locking element. Compressed high-pressure gas (7.5 MPa), which enters the valve through the inlet orifice, tends to lift the ball and direct the gas flow into the Laval nozzle, which provides the creation of reactive force, but this is prevented by a compression spring located in the electromagnet. This spring, with one side resting on a fixed stop, with the other side acts on the armature and, through the rod, pressed into the armature, on the ball. The spring is designed so that the force of its initial compression surpasses the force with which the compressed gas acts on the ball, which is pressed against the locking saddle, as a result of which the gas does not enter the nozzle.
switching ball
valves, controlled by short-stroke high-speed electromagnets, is a few milliseconds. In the process of designing valves, it is necessary to calculate a variety of design options. For example, in the original version (Fig. 2,a) there is a zone of possible saturation of the magnetic core, so it is important to compare the results of the calculation of the dynamics in the original version and in the modified one (Fig. 2,b), where such a zone is absent. In the process of calculations, a different combination of materials (non-corrosive steel / «Armco» steel) is tested, as well as various conditions during the switch-on operation (cold winding / heated winding, maximum allowable voltage / minimum allowed voltage, etc.). For such calculations, it is necessary to use numerical methods for calculating dynamics that adequately reflect real processes in the electromagnet. The basis of such techniques is mathematical computer codes aimed at solving the problems of electromagnetic fields calculations.
nrn
Laval nozzle (gas outlet)
Fig. 1. Sketch of the electromagnetic valve
If a current is passed through the electromagnet winding, the armature is attracted to the stop, the rod moves together with the armature and stops exerting pressure on the ball, which, under the action of compressed gas, moves to the drainage (upper) saddle, releasing the locking (lower) saddle and the gas begins to flow into the Laval nozzle, creating the necessary reactive force.
One of the key parameters for electromagnetic valves is its switch-on speed. The operation time of these
saturation
zone «Armco» steel
non-corrosive steel / «Armco» non-magnetic steel winding
a b
Fig. 2. Sketches of geometrical models of electromagnet for gas
distribution valve: on the left (a) - the initial version, on the right (b) - the modified version
In gas distribution valves, armor electromagnets with penetrating armature are usually used, which with a high degree of adequacy can be considered as 2D axisymmetrical objects in a cylindrical coordinate system [5, 6], which greatly simplifies the technique and speeds up the calculation.
A technique proposed.
Geometry description (development of a geometric model). The use of computer codes in which the Finite Element Method is used requires careful preparatory work related to the description of the problem being solved, namely, the description of the geometry of the electromagnet on the r-z plane passing through the axis of its symmetry; a dividing of the cross-section of the electromagnet and the environment approximated to it (in the aggregate - the calculation domain) into subdomains that can be considered homogeneous in terms of physical properties; specifying the boundaries of these subdomains and differential equations describing the processes in these subdomains; as well as specifying the initial and boundary conditions.
The geometry description (development of a geometric model) of an electromagnet, taking into account the axisymmetrical nature of the calculation problem, is carried out on a plane passing through its axis of symmetry. This axis is located vertically, forming the z
air
axis in a cylindrical coordinate system, and its positive direction must coincide with the direction of movement of the armature of the electromagnet. One of the points on the z axis should be designated as the origin of the calculation domain, then the straight line passing through the point (0, 0) perpendicular to the z axis is the r axis in this coordinate system.
On a specified plane in the preprocessor environment of the corresponding code (FEMM or another code that uses the Finite Element Method), it is necessary to specify the coordinates of the nodal points of the «imprint» of a half axial section of the electromagnet and connect them with straight lines or arcs. The result should be a combination of closed figures that correspond to different parts of the electromagnet. For each of these closed figures, it is necessary to designate the material of which the corresponding part of the electromagnet is made, and in the material library of the corresponding computer code to prescribe the physical properties of this material to be used in the calculations.
The process of preparing a geometric model of an electromagnet described above is quite labor-intensive, however, it can be significantly accelerated if using one CAD software (AutoCAD, SolidWorks, etc.) to draw the mentioned «imprint» without dimensioning, shading, pouring, etc. In this case, this drawing must be placed on a computer «board» so that the origin of coordinates on the drawing and on the computer «board» coincide. If such a fragment of the drawing is saved in the .dxf format, then the resulting file can be imported into the appropriate mathematical computer code, which uses the Finite Element Method. One use this technique of preparing geometric models.
A mathematical model and problem solution.
The defined task is solved as a multiphysics problem involving computer simulation of the following coupled phenomena: transient electromagnetic field distribution in axisymmetrical formulation, transients of current change in the winding, and the movement of the armature of the electromagnet. Attempts of generalization and techniques of various multiphysics problems numerical solution are presented in [7-13].
Electromagnetic calculation is carried out by numerical solution of a system of partial differential equations describing transient magnetic field within the computational domain. The number of equations in the system is equal to the number of subdomains (including the outer space). One of these subdomains, namely the subdomain that corresponds to the cross-section of the winding, differs from the others in that it is influenced by external electric field. All subdomains filled with non-ferromagnetic and non-conductive materials have almost the same electric and magnetic characteristics as air or vacuum. If such subdomains are located nearby, they can be combined into one subdomain. In subdomains filled with non-ferromagnetic, but conductive materials, eddy currents arise, which ones depend on the specific resistance of the material of which the corresponding part is made. Therefore, each of these subdomains must be described separately.
When calculating transient magnetic field, the system of partial differential equations in terms of the
magnetic vector potential A is solved, which, provided that all subregions in the calculation region remain fixed, has the following form [14]:
CA _1
°k • + curl((^Q • ^rk) • curl A) = ck • Ek, Ct
k = 1,2,..., n,
(1)
where k is the number of the subdomain; ak is the electrical conductivity of the material of the k-th subdomain; n is the total number of subdomains, including the outer space; |0 is the magnetic constant; /urk(H) is the relative magnetic permeability of the material of the k-th subdomain depending on the magnetic field strength H; Ek is the vector of the external electric field strength created by an external field source in the k-th subdomain; Jk is the current density vector in the k-th subdomain; t is time.
In electromagnets, external sources act only on the winding, in which the vectors of the external electric field strength and current density have only azimuthal components, therefore the right-hand side of (1) for the subdomain of the winding does not depend on the electrical conductance of its conductive part
Ok ■ Ek = Jk = (i ■ N / S) ■ 10, (2)
where i is the current in the winding; N is the number of turns of the winding; S is the area of the winding window indicated in the initial data: the winding's filling rate (the ratio of the cross-sectional area of copper conductors and the total cross-sectional area of the winding, i.e. the total cross-sectional area of the conductors and insulation) of the winding under consideration is 0.65; 1g is the azimuthal ort - a unit vector directed perpendicular to the plane of the calculation domain.
Since the winding consists of many insulated coils and its electrical conductance in the directions r and z is equal to zero, the equivalent electrical conductance of the winding is usually assigned a zero value and equation (1) for the subdomain of the winding acquires the following form:
curl (curl A) = |0 ■ (i ■ N/ S) ■ 1ff . (3)
In other subdomains, external sources of electric field strength are absent, and in them Ek = 0, and in subdomains occupied by air, also ak = 0, so the equations for sections of the magnetic core are the following: CA _
ak ■^7 + curl ((<"0 ■ "rk ) ■curl A) =0; (4)
Ct , (4)
k = 1, 2, ..., nm and in the subdomains occupied by air or insulation:
curl (curl A) = 0, k = 1, 2, ..., na , (5)
where nm and na are the numbers of subdomains of the magnetic core and air, respectively.
Since vector A has only an azimuthal component A = Aff = A, equations (3)-(5) can be represented in scalar form, but the corresponding mathematical expressions are very cumbersome and are not given here.
One assume that by the time of the beginning of the transient (t = 0 s) the electromagnet is at rest. This means that there is some initial air gap ¿¡0 between the armature and the stop, the winding is not powered and there is no
magnetic field at all points of the computational domain, which corresponds to the zero initial condition for the magnetic vector potential:
A(r , z)| t = 0 = A(r , z, 0) = 0. (6)
To obtain a solution to the system of equations (3)-(5), it is also necessary to formulate conditions at the boundaries of the calculation domain. A natural geometric boundary is the axis of symmetry, at which the magnetic vector potential is zero. Considering the electromagnet as an object located in infinite space, then in solving such problems, one can use the Kelvin transform [15]. At the same time, experience shows that as the distance from the electromagnet moves away, the field attenuates very quickly, so with minor error this problem can be considered with a closed external area, limited in the r-z plane in a semicircle centered at the origin of coordinates. The radius of the semicircle should be 1.5-2 times the distance from the origin of coordinates to the farthest point of the «imprint» of the cross-section of the electromagnet. Thus, the calculation domain is a closed area, divided into n closed subdomains corresponding to individual parts of the electromagnet and the surrounding space within the calculation domain. Due to the fact that the magnetic field is very rapidly weakened, approaching zero with distance from the electromagnet with a practically closed magnetic core, the nature of the boundary conditions practically does not affect the calculation results, therefore zero conditions can be set at all points of the boundary r (see Fig. 2) of the calculation domain for the azimuthal component of the magnetic vector potential A (Dirichlet zero conditions) [14]:
A| r = 0. (7)
Provided that all parts of the electromagnet remain fixed, for some given law of current i variation in time, using well-known numerical methods, it is possible to obtain a solution of the system of equations (2)-(4) taking into account the initial condition (5) and the boundary condition (6). In this work, the problem is solved by the Finite Element Method, and for discretization in time domain the implicit multistep Adams-Moulton method of the 4th order with adaptive selection of time step values is used. On each time step the system of equations (2)-(4) and the problem of dynamics of motion descrbed below are solved sequentially with the same numerical value of the time step for both indicated problems. The numerical solution is obtained in the form of a table of values of the magnetic vector potential at different points in time at the nodal points of the mesh, covering the of calculation domain, however, using various methods of interpolation, numerical differentiation and integration, it is possible to build a graphical picture of the force lines at different temporal points, as well as to calculate the values of a number of physical quantities, in particular, the values of magnetic flux density, magnetic flux, the force acting on individual subdomains, etc. at these times.
The solution of the system of equations (3)-(5) for a given law of current variation in the winding corresponds to the problem of calculation of the electromagnets with serial windings, which are used, for example, in current relays. In the absolute majority of electromagnets, including in the driving electromagnets of gas distribution
valves, shunt windings are used, which are connected in parallel to the power sources of control circuits. When solving the problem of calculating the transient magnetic field of such electromagnets, provided that all its parts remain fixed, the system of equations (2)-(4) should be supplemented with a non-linear differential equation describing the transient of current change in the winding:
U + e = R ■ i + L ■ —, (8)
dt
where U is the winding power source voltage; R is the resistance of the winding circuit; L is the inductance of the external circuit; e is the back EMF arising in the winding:
„ N fdA J J e = -2 ■ n-----r ■ dr • dz,
s idt
(9)
When solving the system of equations (2)-(4), (7), taking into account dependence (8), zero boundary conditions (initial conditions (5) and conditions (6) on the boundary of the calculation domain) are used for the azimuthal component of the magnetic vector potential, and for the current i in the winding zero initial condition should be set:
i| , = o = i (0) = 0. (10)
If the power source of the shunt winding is a capacitor with a capacitance C, pre-charged to the voltage UC0, then in equation (7) the value U must be replaced by the current value of the voltage uC on the capacitor:
uC + e = R • i + L •—, (11)
dt
and the system (2)-(4), (10) with regard to (8) must be supplemented by additional differential equation:
C ■ ^ = -. dt
(12)
In the process of solving the system (3)-(5), (8) taking into account dependence (9) or system (3)-(5), (11), (12) taking into account (9), it is necessary to control the value of the traction force F, acting on the fixed armature of the electromagnet, and compare it with the initial value of the opposing force Fr0. Until the force F is less than the opposing force Fr0, the electromagnet is at rest and the armature is remain fixed.
The gas distribution valve is a precision device - the full stroke of its armature is 0.26 mm, and the stroke of the switching ball, from the stop saddle to the drainage one, is 0.19 mm. The characteristic of opposing forces acting on the electromagnet armature is shown in Fig. 3. The initial force of compression of the spring is equal to Fs0 = 40 N, and the initial opposing force is only Fr0 = 21 N, since the force created by the compressed gas and acting on the ball «pushes» it upwards and reduces the force opposing the movement of the armature. Thus, the valve will not operate if the electromagnetic force is less than 21 N in the presence of compressed gas at the valve inlet or when the force is less than 40 N without gas.
The numerical solution of the problem of calculating the force with which the compressed gas acts on the ball at its different positions is obtained by the Finite Element Method based on solving the Navier-Stokes system of equations for the valve model (see below). Compressed air heated to temperature of 1200 °C is considered as a moving
medium. The boundary conditions are the zero speed on the surface of the ball, the rod and on the inner walls of the valve, as well as the value of air pressure at the valve inlet (7.5 MPa) and exit to the atmosphere (0.1 MPa).
Fr, N 45
40
35
30
25
20
15
10
5
Fs
Fs0
Fr
Fr0 Fb__
flow channel, the axial symmetry conditions, the inlet temperature which is assumed to be 1200 °C, and the open boundary condition at the outlet:
k-W-3 = P-Cp ■ u3. (16)
The force acting on the switching ball is defined as the surface integral of the forces caused by pressure, viscosity and speed pressure:
(
Fz =
( duz dur - nzP + nrri\ — + — ^ or oz
+ 2nzr
OUz dz
■ dS,
(17)
0 0,05 0,10 0,15 0,20 0,25 ¿mm
Fig. 3. Characteristics of the opposing forces acting on the electromagnet armature: Fs is the opposing force of the spring;
Fb is the force with which the compressed gas acts on the switching ball; Fr is the total opposing force that occurs when gas is supplied to the valve inlet (in the absence of gas at the valve inlet, only the spring with the initial force Fs0 opposes the movement of the armature)
The system of Navier-Stokes equations in the absence of volumetric forces has the following form in the matrix representation:
p\u ■ V)-u = V■ [-p ■ I + u + (V-u)T)-
(2 ^ ( ) (13)
-[ 3-?J-(v-u)-i];
V(pu)= 0; p = p(p,3); 7 = 7(3), (14)
where p is the gas density; u is the gas speed vector matrix; p is the pressure; I is the unit matrix; 7 is the gas viscosity; 3 is the temperature (here and above units of all quantities are used in the Si system).
Equation (13) is the equation for the steady flow of a gas stream, and expressions (14) are conditions for continuity of flow.
The boundary conditions for equations (13), (14) are: zero speed values on the walls of the elements, the condition of axial symmetry is that the normal speed component is zero, and the pressure values on inlet (7.5 MPa - compressed gas pressure) and outlet (0.1 MPa - atmospheric pressure).
The thermophysical characteristics of a gas (density and viscosity) substantially depend on temperature, therefore the system of equations (13), (14) must be supplemented with the equation of heat conduction for a gas in a stationary mode. Since gas is a compressible medium, it is necessary to take into account pressure work in the heat conduction equation - the ability of a gas to produce work by pressure, which significantly affects the temperature of the medium. The corresponding heat conduction equation in this case has the following form: V-(-V-3) = p■Cp u-V-3 +
+№i+-v-p) (15)
where 2 is the thermal conductivity of gas; Cp is the heat capacity at constant pressure.
The boundary conditions for equation (15) are: the absence of heat transfer along the boundaries of the gas
where n is the ort of appropriate direction.
A sketch of the valve system and its geometric model are shown in Fig. 4. The outlet in the valve has the shape of a cylinder, and in the model it has the shape of a disk with the same cross-sectional area in the direction of movement of the compressed gas (in the model - in the radial direction).
A specific feature of this problem solution is that the sought quantities (the temperature 3, the speed u and pressure p at all nodes of the mesh covering the calculation domain) are complex interdependent functions, and the compressed gas characteristics (the specific heat conductivity 2, the density p, the viscosity 7, and the specific heat Cp) in individual nodes of the mesh can differ from each other by several orders of magnitude. The use of simple iterative procedures for solution of problems in which the convergence of the iterative process strongly depends on the initial values of the sought quantities, as a rule, does not lead to success. Our problem belongs to the category of exactly such ones: at unsuccessfully selected initial values of the sought values, the iterations do not converge, and at successfully selected ones they converge. Therefore, in the calculations we applied more complex iterative procedures, which, finally, led to the convergence of iterations. The essence of these procedures is that the solution to the problem is divided into two stages: 1) preliminary calculations to determine the initial values of the sought quantities in the iterative process of joint solution of equations (13)-(15), and 2) the iterative process of joint solution of equations (13)-(15) with initial the values of the sought quantities found in step 1). in turn, stage 1) is divided into two sub-stages: 1a) solving the problem of calculating the distribution of the speed and pressure in the calculation domain at a given initial temperature distribution, and 1b) solving the problem of calculating the temperature distribution in the calculation domain at the found distribution of the speed and pressure. With this approach to solving this problem, the calculation algorithm looks like this:
1. We set a fixed temperature distribution (for example, the temperature at all mesh nodes in the computational domain is taken equal to the ambient temperature) and solve equations (13), (14), as a result of which values of the temperature and speed are obtained at all mesh nodes which are considered as new initial values for the next iteration, etc. Here, the results at the previous and subsequent iterations are compared and the calculations continue until the changes in the pressure and speed in the module in all the calculation nodes become smaller than some small positive predefined values.
+
0
2. Having fixed the obtained values of the pressure and speed, the heat equation (15) is solved using the iteration method, where at each iteration at each node the compressed gas characteristics are determined (the specific thermal conductivity 2, the density p, and specific heat Cp) as a function of temperature values, found in the previous iteration. The results of calculating the temperature at all nodes of the mesh are considered as the initial values for the next iteration and the calculation is repeated until the temperature changes in the module in all the calculation nodes become less than some small positive values given in advance.
IN
3. Having accepted the obtained values of the pressure, speed and temperature at all nodes of the mesh as the initial ones, we solve the equations of gas flow motion and thermal conductivity (13)-(15) jointly using the iteration method, where at each iteration in each node values of compressed gas characteristics (the thermal conductivity 2, the density p, the viscosity j, and the specific heat Cp) are determined as a function of pressure, speed and temperature values found at the previous iteration, until the changes in the pressure, speed, and temperature in the module in all calculation nodes become smaller than some small positive values preset in advance.
/> = 0.1 MPa
ta
S ®
II
p = 7.5 MPa
a b
Fig. 4. Sketch of the valve system (a) and its geometric model (b)
The results of the calculation of the forces acting on the ball in its various positions are shown in Table 1.
Table 1
The results of the calculation of the forces acting on the ball
s, mm 0.00 0.05 0.10 0.15 0.19
F, N 18.25 21.46 21.13 23.39 23.45
Notes:
1. The path s is measured from the initial position of the ball.
2. The force F has only the axial component.
Figures 5, 6 illustrate some intermediate results of the calculations: a picture of the pressure field of the compressed gas on the ball and a picture of the speed field of the compressed gas in the nozzle at the outlet of the valve.
When performing calculations of the dynamics with a moving armature for the subdomain corresponding to the armature, equation (4) and expression (9) should be replaced by expressions that take into account the movement of the armature:
dA
dt
• M-rar )
'm
S i\~0t
- varz x curlA | +
(18)
-1
+v„,
curlA) = 0:
dA ^
• r • dr • dz,
(19)
2 N fl e = -2 ■ n---I , varz
•U dt ^
where ct^, ^raj are the electrical conductance and relative magnetic permeability of the material (steel) of which the armature is made; varz is the speed vector of the armature in the direction of z coordinate.
p, MPa 7
4
Fig. 5. Picture of pressure field of the compressed gas on the ball at s = 0.19 mm (final ball position)
v, m/s 1000
0
Fig. 6. Picture of the speed field of the compressed gas in the nozzle at the valve outlet at s = 0.10 mm (close to the average position of the ball)
6
5
When calculating the dynamics in the process of motion, the system of equations describing the magnetic field and current in the winding is supplemented by the equations of motion, which, with a constant moving mass m, has the form:
m ■
dVa
dt ds dt
■ = F - Fr
= Va
(20)
(21)
with zero initial conditions:
Varz(t0) = 0; s(t0) = 0. (22)
Here, the following notation is used: F is the electromagnetic force that ensures the movement of the armature; Fr is the force opposing the movement of the armature; varz is the module of the vector varz; s is the path traversed by the armature after the start of the movement; t0 is the point in time corresponding to the beginning of the movement.
The fixed and moving subdomains appear to be located in a certain elastic space covered with a mesh.
When the position of the moving subdomain changes, the meshes belonging to the fixed and moving parts remain unchanged.
Similar algorithms based on the use of deformed elastic finite element meshes are used in the numerical solution of hydrodynamic problems [16, 17]. The mesh belonging to the elastic subdomain is deformed in accordance with the law of motion of the armature, which is defined by equations (20), (21). With a significant change in the gap, the angles opposite the bases of some mesh cells become close to 180°, which reduces the accuracy of the calculations.
In this work, an algorithm is used that involves monitoring of the angles of the deformable mesh cells and stopping the computation when one of the angles becomes greater than a certain predetermined value (in our case, 170 °C). After stopping, the code rearranges the mesh and continues the calculation.
Figure 7, as an example, shows a fragment of computational subdomains covered with a triangular mesh. As can be seen, the meshes on moving and fixed subdomains, corresponding to different parts of the electromagnet magnetic core, are not deformed during movement of the armature due to the invariance of the shape of these subdomains, and the mesh covering the subdomain surrounding the armature is elastic and significantly deformed.
»A4«
moving armature
T\ i \ / \ /xi AAPmsI
a
moving armature
I \l I / I —■—
b
Fig. 7. A fragment of the calculation domain covered with a triangular mesh: a - armature in the initial position; b - armature in the end position
Numerical results and their analysis.
Figure 8 shows the results of the calculation of a modified version of an electromagnet corresponding to the geometric model shown in Fig. 2,b. in Fig. 8,a, in addition to the picture of the field distribution, the main dimensions of the electromagnet are shown: the gap S varies from 0.3 mm in the initial position to 0.04 mm in the final position, thus the armature stroke is 0.26 mm. in Fig. 8,b, in addition to the field distributions, the main subdomains are indicated: 1 - electrical steel, 2 - winding, 3 -insulation («air»), 4 - non-magnetic steel, as well as a scale of values of magnetic flux density B is shown.
Magnetic field distributions presented in Fig. 8, show the process of penetration of the magnetic field into the magnetic core along the entire perimeter, and the magnetic flux density in the area of the working gap reaches large values (1.8 T) within 5 ms after the start of the process.
Figures 9,a,b show the results of the calculation of the basic dynamic characteristics of this electromagnet under the following operating conditions: the armature is
made of electrical steel, the number of windings is 1155, the winding is heated to 112 °C, the resistance of the heated winding is 55 Q, the supply voltage is minimal allowable: 24 V. Figure 9,a shows the calculation results for the case when gas is supplied to the valve inlet, and Fig. 9,b - when the gas is not supplied. in the first case, the operating time is 9 ms, and in the second case 13.5 ms. it can be seen that the modified electromagnet operates in both cases. in the original version of the electromagnet with a saturation zone corresponding to the geometric model shown in Fig. 2,a, the operating time when there is a supply of compressed gas to the valve inlet is 9.5 ms, and without gas supply the electromagnet does not operate.
The requirements for the speed of gas distributors are put forward not only when the winding is switched on, but also when it is switched off. Switching off the winding with the help of contact switching elements is accompanied by their sparking, which can cause radio interference and a negative effect on the electronic
devices of the control object. Switching off the winding using semiconductor switching elements requires
overvoltage protection. Typically, such protection is
implemented by shunting the winding with an «opposite» diode that does not conduct current when the winding is energized, and automatically opens when the switching element stops energizing the winding. After this moment, the current in the winding, and, consequently, the force of attraction of the armature to the stop, begin to decrease,
but for some time the armature remains in the attracted state until the force of attraction exceeds the opposing force. At the moment when the force of attraction is compared with the opposing force, the movement of the armature begins in the opposite direction. This movement continues until the armature returns to its original position. The time interval from the moment when the winding supply ceases to the moment when the armature movement ends is called the release time.
I
4,4
J17
B,T
1,6
1,4
1,2
1,0
0,8
0,6
0,4
0,2 0,0
F,N 100 80
i, A 0,6 60
0,4 40
0,2 20 0 0
(H5,6N)
F
s
À . i
//
F, N 100
i 80
^ i, A 0,3 0,6 60
0,2 0,4 40
0,1 0,2 20
0 0 0
(fî 5,6 N
F
s
i
0,3 0,2 0,1 0
0 5 10 15 20 25 t, ms 0 5 10 15 20 25 t, ms
a b a b
Fig. 8. Distribution of the magnetic flux density Fig. 9. Basic dynamic characteristics of a modified version of an electromagnet module in 1 ms after the beginning of the when performing the switching on operation: a - when supplying compressed gas transient process (a); the same in 5 ms (b) to the valve inlet; b - the same without gas supply
Calculations of dynamic processes during the operation of switching off the winding, shunted by an «opposite» diode, are performed for the case when the
operation of switching off occurs when there is a supply of compressed gas to the valve inlet. In this case, the release time is 16 ms (Fig. 10,a).
F, N 100
80
60
i, A 0,4 40
0,2 20
0 0
F
V s
i s
0,3
F, N 100
80
60
i, A 0,2 0,4 40
0,1 0,2 20
0
10
0
15 t, ms
0 0
■ F
s
\ \
i
0,3 0,2 0,1 0
0
5 10 15 t, ms
a b
Fig. 10. Dynamic characteristics of a modified version of an electromagnet when performing the switching off operation of the winding, shunted by the «opposite» diode: a - without additional resistor; b - with a resistor with resistance of 56 Q
8
It is possible to reduce the release time connecting an additional resistor in series with the diode. This reduces the time constant of the winding circuit, and accelerates the process of reducing current, magnetic flux and force. This is confirmed by the results of the calculations presented on Fig. 10,b: in the presence of an additional resistor with resistance of 56 Q, the release time is reduced from 16 ms to 9.5 ms which corresponds to experimental results.
Conclusions.
A numerical technique and results of the calculation of the dynamics of a high-speed armored-type electromagnet used in micro-motors of control systems for the orientation of rockets are presented. The multiphysics computer simulation involves numerical analysis of the following coupled phenomena: transient electromagnetic field distribution in axisymmetrical formulation, transients of current change in the winding, and the movement of the
armature of the electromagnet. All corresponding differential equations are solved with regard to appropriate initial and boundary conditions as well as taking into account nonlinear material properties for all structural parts of electromagnets. Electromagnetic field calculations are carried out by the Finite Element Method involving the use of so-called tunable elastic meshes. Peculiarities of the numerical algorithms used are described in detail. Comparative analysis of numerical results obtained for a few electromagnet designs is presented.
REFERENCES
1. Olejnik V.P., Yelanskyi Yu.A., Kaluger L.G. Mathematical modelling of a gas distributor of the carrier rocket gas-jet control system. Space Technology. Missile Weapons, 2017, iss. 1 (113), pp. 59-66. (Rus).
2. Beliaev N.M., Belik N.P., Uvarov E.i. Jet Control Systems for Spacecrafts. Moscow, Mechanical Engineering Publ., 1979. 232 p. (Rus).
3. Bajda Ye.i., Klymenko B.V., Pantelyat M.G., Korol O.G., Yelanskyi Yu.A. Peculiarities of calculating forced electromagnets shunt windings heating in transient modes. Proceedings of the 18th International IGTE Symposium on Numerical Field Calculation in Electrical Engineering. Graz, Austria, September 2018, p. 31-36.
4. Bajda Ye.i., Klymenko B.V., Pantelyat M.G., Yelanskyi Yu.A., Trichet D., Wasselynck G. Peculiarities of calculating the dynamics of high-speed electromagnets using tunable elastic meshes. Proceedings of the 18th International Symposium on Electromagnetic Fields in Mechatronics, Electrical and Electronic Engineering (ISEF'2019). Nancy, France, August 2019. 6 p.
5. Zhang J.Z., Cai C.H., Wu C.G. Design and analysis of a new permanent magnet actuator for medium voltage vacuum circuit breakers. Applied Mechanics and Materials, 2013, vol. 313-314, pp. 20-26. doi: 10.4028/www.scientific.net/amm.313-314.20.
6. Bissal A., Magnusson J., Salinas E., Engdahl G. Multiphysics modeling and experimental verification of ultra-fast electromechanical actuators. International Journal of Applied Electromagnetics and Mechanics, 2015, vol. 49, no. 1, pp. 5159. doi: 10.3233/jae-140176.
7. Bajda Ye.i., Klymenko B.V., Pantelyat M.G., Trichet D. Wasselynck G. Electromagnetic and thermal transients during induction heating of cylindrical workpieces. Acta Technica, 2018, vol. 63, no. 5, pp. 657-682.
8. Pantelyat M.G. Multiphysical numerical analysis of electromagnetic devices: state-of-the-art and generalization. Electrical Engineering & Electromechanics, 2013, no. 3, pp. 29-35.
9. Pantelyat M.G. Multiphysics in electromagnetic devices simulation and design: an attempt of generalization. Acta Technica, 2012, vol. 57, no. 2, pp. 127-142.
10. Pantelyat M.G., Shulzhenko N.G., Matyukhin Yu.i., Gontarowsky P.P., Dolezel i., Ulrych B. Numerical simulation of electrical engineering devices: magneto-thermo-mechanical coupling. COMPEL: The International Journal for Computation and Mathematics in Electrical and Electronic Engineering, 2011, vol. 30, no. 4, pp. 1189-1204. doi: 10.1108/03321641111133127.
11. Pantelyat M.G., Feliachi M. Magneto-thermo-elastic-plastic simulation of inductive heating of metals. The European Physical Journal Applied Physics, 2002, vol. 17, no. 1, pp. 2933. doi: 10.1051/epjap:2001001.
12. Podoltsev O.D., Kucheriava i.M. Multiphysics modeling of electrotechnical devices. Technical Electrodynamics, 2015, no. 2, pp. 3-15 (Rus).
13. Podoltsev O.D., Kucheriava i.M. Multiphysics modeling in electrical engineering. Kyiv, The institute of Electrodynamics of the National Academy of Sciences of Ukraine Publ., 2015. 305 p. (Rus).
14. Stratton J.A. Electromagnetic Theory. NJ, Wiley, 2007. 640 p.
15. Meeker D. improvised open boundary conditions for magnetic finite elements. IEEE Transactions on Magnetics, 2013, vol. 49, iss. 10, pp. 5243-5247. doi: 10.1109/tmag.2013.2260348.
16. Cizmas P.G.A., Gargoloff J.i. Mesh generation and deformation algorithm for aeroelasticity simulations. Journal of Aircraft, 2008, vol. 45, no. 3, pp. 1062-1066. doi: 10.2514/1.30896.
17.Dwight R.P. Robust mesh deformation using the linear elasticity equations. Computational Fluid Dynamics, 2006, pp. 401-406. doi: 10.1007/978-3-540-92779-2_62.
Received 20.04.2020
E.I. Baida1, Doctor of Technical Science, Associate Professor, B.V. Klymenko1, Doctor of Technical Science, Professor,
M.G. Pantelyat, Candidate of Physics and Mathematics,
Associate Professor,
Yu.A. Yelanskyi2,
D. Trichet3, Professor,
G. Wasselynck3, Assistant Professor,
1 National Technical University «Kharkiv Polytechnic Institute»,
2, Kyrpychova Str., Kharkiv, 61002, Ukraine,
e-mail: baida.kpi@gmail.com, m150462@yahoo.com
2 Yuzhnoye State Design Office,
3, Krivorozhskaya Str., Dnipro, 49008, Ukraine, e-mail: yuyelans@gmail.com
3 IREENA, Polytech'Nantes, University of Nantes, 37 boulevard de l'université BP406,
44600 Saint-Nazaire, France, e-mail: didier.trichet@univ-nantes.fr, guillaume.wasselynck@univ-nantes.fr
How to cite this article:
Baida E.i., |Klymenko B.V.|, Pantelyat M.G., Yelanskyi Yu.A., Trichet D., Wasselynck G. Challenges of dynamic simulation of high-speed electromagnetic valves of gas distribution devices. Electrical engineering & electromechanics, 2020, no. 5, pp. 3-11. doi: 10.20998/2074-272X.2020.5.01.
