VISCOUS AND MAGNETIC FORCE ACTING IN A HIGH SPEED SOLENOID VALVE

Viscous and the magnetic force acting on the armature of a high speed solenoid valve (SV) during its closure movement were identified numerically using a novel computational method. The numerical solutions were compared to the analytical ones. The derivation steps of the analytical formulas were presented. The convergence of the numerical solution was shown in figures. The deviation between the numerically and analytically obtained forces was listed in tables. The distribution of the calculation error was presented in figures by means of the magnitude of the magnetic flux density and the viscous stress.


INTRODUCTION
SVs are very widely set in the automotive industry.Also the domestic washing machines and the dishwashers use SVs to control water entry into the machine.They can be used for a wide array of industrial applications, including general on-off control, calibration and test stands, pilot plant control loops and process control systems.There are lots of reports regarding finite-element approaches for example [2,19] and simulation models [13] of SVs as well as investigation of temperature distribution and thermal deformations inside the [1].Also effects like the estimation of the inductive power caused by the pulse width modulation of the coil current and its influence on the time dependent temperature distribution in the armature of a SV were recently reported by in [7].In [19] the author presented a research of key parameters influencing the driving electromagnetic force of a SV.The driving force causes the movement of the SV armature and decreases the axial magnetic gap.In conventional linear solenoids the armature slides during its movement on the inner surface of the sleeve (see Fig. 1).The sleeve is usually made of a paramagnetic material e.g. an aluminum alloy.The objective of the us-age of a paramagnetic material is the insertion of the second magnetic gap in the magnetic circuit.This gap is called the radial magnetic gap.In most studies on SVs the thickness of the radial gap is kept constant independent from the circumferential angle.In other words researchers take an assumption that the armature is placed concentrically in the sleeve.This assumption implies the homogeneously distribution of radial magnetic forces over the armature circumference.Because of e.g.manufacturing imperfections the armature is in fact positioned eccentrically in the sleeve.In this case the distribution of radial magnetic forces over the armature circumference is not homogeneous and the resulting radial magnetic force attracts the armature toward the inner side of the sleeve.Thus, an increased transverse force acting on the armature exists and causes friction between the armature and the sleeve.Friction with these components degrades the performance of the solenoid and causes wear.The investigation of the transverse magnetic force has been rarely reported.However, one can find some inventions having an objective to minimize the transverse force that occurs from armature eccentricity.In [3] the inventor described an idea of a reduction of the transverse force by the introduction of segmenta-tion in the armature member.Another inventor [6] had an idea of minimizing the transverse force by an insertion of small radial slots in the armature side surface.The objective of the current study is the estimation of the transverse magnetic force acting on the armature with the cylindrical shape which is positioned eccentrically in the sleeve.For this purpose a sophisticated computational method was developed and published in [23] and [9].This method considers a mathematical way of transforming the geometry of the radial air gap into rectangular computation domain.In this computation domain the Laplace's differential equation can be solved for finding the distribution of the magnetic scalar potential and finally the transversal magnetic force.Exactly the same mathematical algorithm can be used for solving of the Poisson's differential equation having an objective of the estimation of the distribution of the oil velocity in the area between the armature and the sleeve (see Fig. 1).For a known velocity distribution the viscous force acting on the armature during its axial movement can be computed.

DEFINITION OF THE PROBLEM
The object of the investigation is the determination of forces acting on the armature of a SV, which simplified layout was shown in Figure 1.The armature is not positioned concentrically in the sleeve but it is shifted toward the sleeve by the distance e.This distance is called the armature eccentricity.The sleeve is made of a paramagnetic material.The room between the armature and the sleeve is filled with oil having the temperature T. The armature moves in the z-direction with the constant velocity u 0 (T).
During this movement acts on the armature the viscose friction force which can be estimated by the integration of the viscous stress tensor over the armature side surface [16]: The function dA in (1) is the vectorial surface element of the armature side surface.For the viscous stress tensor  yields [10]: The parameter η in (2) is the dynamic viscosity, u is the velocity of oil and is the unit matrix.The oil velocity needed in (2) can be estimated by solving the stationery Navier-Stokes differential equation [20]: The vector function f in (3) described the mass forces, ρ is the oil density, p is the pressure and υ is the kinematic viscosity of the oil.In the investigated SV the mass forces were neglected and it was assumed that the oil density and the pressure drop are constant.Furthermore, the armature of the SV was approximated with an infinite long cylinder.In such of an approximated case (3) simplifies to the Poisson's differential equation [12]: The constant C in (4) equals: (5) The parameter L h is the axial length of the armature (see Fig. 1).Because of the taken assumptions equation (1) simplifies to: (6) where the function τ in ( 6) is the viscous stress vector on the side surface of the armature, for which yields: Except of the viscous force that acts on the armature, also an electromagnetic force, which consists of the driving component F mz e z and the transversal component F mx e x is present.In the considered investigation however only the transversal component is relevant.This component can be estimated using the integration [22]: (8) The function B in ( 8) is the magnetic flux density in the radial air gap and μ 0 is the permeability of the vacuum.In a magneto-static case the magnetic field density existing in the domain r∈[ρ 1 , R 3 ] can be represented by the gradient filed according to [5]: The function ψ is the magnetic scalar potential of the magnetic field intensity.Because in the investigated SV the total drop of magneto-motoric force θ was assumed to consist of the sum of drops in the radial and in the axial air gap, the equation ( 9) gets with the use of the Hopkinson's law [15]: (10) The parameters G 1 and G 2 in (10) are the magnetic permeances of the axial and the radial air gaps and Ψ is the normalized magnetic scalar potential.After the use of the Gauß rule [21] the divergence of the magnetic field density is equal zero and one obtains the Laplace differential equation which is a special case of the Poisson's differential equation [14]: (11) The permeance of the axial air gap can be found using the analytic formula (12) for the room between two parallel surfaces [14]: The parameter h b in (12) is the length of the axial air gap and A d is the cross section of the armature.For the permeance of the radial air gap yields [22]:

NUMERICAL MODEL
The equation ( 4) was solved numerically using the method of finite differences.The co-ordinate system in each this equation was solved is the curvilinear system a, α.
In this co-ordinate system the distance from the armature contour (described in the universal case by the function ζ(φ)) to the inner side of the sleeve (described in the universal case by the function ζ(φ)+h(φ)) is constant independently from the circumferential position φ (see.Fig. 2).The function h is here the local thickness of the oil film.The transformation from the polar coordinate system r, φ to the a, α -co-ordinate was performed using the relations [9,23]  Laplace operator transformed to the a, α -co-ordinate system has been derived in [23].The derivation was performed using a shoal of differential function.It equals: (16) where: The computation method in the a, α -co-ordinate system is a relatively novel computation method [23,9].It has this advantage in comparison to the Finite Difference Method (FEM) that regions with small magnetic or hydraulic gap in which both magnetic flux density and viscous shear stress are the highest can be automatically mesh more densely than regions with big gaps.This method of discretizing allows an incensement of computation precision with a simultaneous reduction of the mesh node numbers.The Laplace operator ( 16) can also be derived using its generalized definition that allows the operations on functions defined on surfaces in Euclidean space.This operator of a scalar function in any curvilinear co-ordinate system is called the Laplace-Beltrami operator.In the considered case this operator simplifies to [4,17,18]: (18) where: g is the contravariant metric tensors which can be found using [4]: (19) X in (19) is the position vector: (20) The derivatives of (20) are: (21) (22) where: Setting ( 21) and ( 22) into (19) yields: The determinant of the metric tensor (19) equals: (29) Setting ( 25)-( 29) into (18) and performing differentiations yields and Laplace operator identical to (16).The full derivation using this method has been recently presented in [9].One of possible ranges of validity for the variables a, α is 0 ≤ n ≤ 1 and 0 ≤ α ≤ 2π.In the investigated case the location a = 0 indicates points on the armature side surface and a = 1 points on the inner side of the sleeve.After having defined the Laplace operator the equation ( 4) has been solved with the boundary conditions: The parameter u 0 in (30) is the axial velocity of the armature.In order to estimate the viscous share force in (6) one needs the nabla operator [9]: and the vectorial side surface of the armature [9]:  Similar simplifications can be also done for (32) and (33).The local thickness of the oil film can be described by the function [23]: The differentiation of (35) yields: Exactly the same computation algorithm was used for the estimation of the transversal component of the magnetic force.Solved was the equation ( 11) with the boundary conditions: The only difference in the use of the Laplace operator (34) is the function h.In the magneto-static computation the function h represents the circumference dependent radial air gap.For this computation the inner radius of the sleeve R 2 in equations ( 35) to (37) must be replaced with the outer radius of the sleeve R 3 .

ANALYTICAL ESTIMATION OF THE VISCOUS AND THE MAGNETIC FORCE
For e≪R 2 the local oil film thickness (35) can be simplified to: The parameters δ h and ε h in (40) are the hydraulic clearance and the relative hydraulic eccentricity which are defined as: In the case of thin oil film δ h ≪ R 1 the derivative of the oil velocity in radial direction is dominant over the derivative in the peripherical direction: r -1 ∂u z /∂φ ≪ ∂u z /∂r so it can be assumed that: Under the assumption (43) the Poisson's equation ( 4) simplifies to: The boundary conditions of (44) are: Setting the relations (61), ( 58) and ( 56) into (8) yields: (62) where: (63) The parameter L m (62) in is the axial length of the yoke pole (see Fig. 1).The integration (63) and the integration of (62) was presented in [8].For the transversal component of the magnetic force yields: (64)

CALCULATION RESULTS
The numerical and analytical calculations of both fluid mechanical and magneto-static forces were performed for the parameters showed in Table 1.The calculations of viscous shear force were done for two operating temperatures of the SV.The calculations were performed for three different configurations of the SV indicated in Table 2 with: P 1 , P 2 , P 3 .The nominal air gap δ m = 300 μm according to (59) in the case of the configuration P 1 is identical to the configuration P 2 .The nominal hydraulic clearance δ h = 10 μm according to (59) in the case of the configuration P 1 is identical to the configuration P 3 .The eccentricity e in Table 2 is different for each evaluation point.The current eccentricity is a result of an equilibrium position of the armature.The equilibrium occurs when roughness peaks of the armature stay in contact with roughness peaks of the sleeve.The pressure existing in this contact is a function of the eccentricity and can be calculated using the model presented in [11].
The increase of the eccentricity results in the increase of the contact pressure.The resulting contact force is the reaction to the transversal component of the magnetic force.However, the presented investigations do not concern the estimation of the armature eccentricity, so the value e ought to be considering as a given calculation input.The armature cross section A d (used in (12)) is less than the area of the top of the cylinder shown in Figure 1.This is due to the fact, that the calculation was perfumed for SV in which the armature has some axial holes used for the mass reduction.An example of such the holes can be seen in [6].
In the first computation step the convergence of the viscous share force in all evaluation points and at both the investigated temperatures was examined.The computed viscous force depends on the number of mesh points in both radial N r and pheriperical direction N p as well as on the number of numerical iteration N ν .The numerical iteration was run as long as the maximum value of the residuum of the oil velocity R s (a, α) has converged and reached the value smaller than 1E-16 (see   For this converged velocity distribution the complete 2-dimensional convergence matrix F hz (N r , N p ) have been examined.In all evaluation points and at both the investigated temperatures the dependence of the viscous share force on the number of mesh points in radial direction was much smaller than its dependence in pheriperical direction: ∂F hz /∂N r ≪ ∂F hz /∂N p .An example of the convergence of the viscous share force was shown in the upper plot of Figure 3.The force estimated by the analytical solution is plotted with the straight line while the force obtained by numerical computation in plotted with square points.The number of mesh points in pheripherical direction was increased from ten to ten thousand.In this convergence example the number of mesh point in radial direction was equal six.An almost identical convergence line was obtained for natural numbers of N r ∈(4, 10).The maximal residual error falls below 1E-16 already for N ν = 64.After about N ν = 100 the variable max(R s ) has converged.
In the further computation step the convergence of the transversal magnetic force in all evaluation points was examined.Also a full mesh independace study was performed.In the considered model there is no temperature dependency of the magnetic force.An example of the convergence of the transversal magnetic force was shown in Figure 4.In his example the maximal residual error falls below 1E-16 for N ν = 92.After about N ν = 140 the value max(R s ) has converged.
In Figure 5 percentage error distributions of the magnetic flux density and the viscous stress between analytical and numerical solutions in the evaluation points: P 1 , P 2 , P 3 was shown.In the case of the magnetic flux density the maximal value of the percentage error in the evaluation point P 2 is the biggest.In this evaluation point the relative magnetic eccentricity ε m according to (60) is much higher than in the other evaluation points.In the case of big relative magnetic eccentricity the accuracy of the formula (61) is getting less.However the percentage error of the  3).
In the case of the viscous stress, the maximal value of the percentage error in the evaluation point P 2 is also the biggest at both T min and T max (see Fig. 6 and Fig. 7).In this evaluation point the hydraulic clearance δ h according to (41) is much In the evaluation points P 2 and P 3 the percentage errors of the viscous force are less than 1%.

CONCLUSIONS
1.The viscous share force at T min in the evaluation point P 2 is about three times lower than in the evaluation points P 1 and P 3 .However, this force at T max it is about 50% higher than in P 1 and P 3 .2. The magnetic force in the evaluation point P 1 is the lowest.It reaches its biggest value in the evaluation point P 2 which it is about 24-times higher than in P 1 .The magnetic force in the  evaluation point P 3 is about 3-times higher than in the evaluation point P 1 .3. The deviation between numerical computations and analytical approximations in the case of the magnetic force lies under 1% for all the evaluation points.In the case of the viscous share force stays the error under 1% only in P 1 and P 3 .For P 2 (in each the eccentricity is the biggest) the error equals about 7%.

Fig. 3
Fig.3lower plot).For this converged velocity distribution the complete 2-dimensional convergence matrix F hz (N r , N p ) have been examined.In all evaluation points and at both the investigated temperatures the dependence of the viscous share force on the number of mesh points in radial direction was much smaller than its dependence in pheriperical direction: ∂F hz /∂N r ≪ ∂F hz /∂N p .An example of the convergence of the viscous share force was shown in the upper plot of Figure3.The force estimated by the analytical solution is plotted with the straight line while the force obtained by numerical computation in plotted with square points.The number of mesh points in pheripherical direction was increased from ten to ten thousand.In this convergence example the number of mesh point in radial direction was equal six.An almost identical convergence line was obtained for natural numbers of N r ∈(4, 10).The maximal residual error falls below 1E-16 already for N ν = 64.After about N ν = 100 the variable max(R s ) has converged.

Fig. 3 .Fig. 4 .
Fig. 3. Convergence of the viscous share force for the evaluation point P 3 at T min (upper plot) and the course of the residual error

Fig. 5 .
Fig. 5. Percentage error distributions of the magnetic flux density and the viscous stress

Fig. 6 .
Fig. 6.Percentage error distributions of the viscous stress at T min

Table 2 .
Analyzed evaluation points

Table 3 .
Numerically computed forces and their percentage error regarding to analytical calculation