Energy Based Calculation of the Second-Order Levitation in Magnetic Fluid

: A permanent magnet immersed in magnetic ﬂuid experiences magnetic levitation force which is of the buoyant type. This phenomenon commonly refers to self-levitation or second-order buoyancy. The stable levitation height of the permanent magnet can be attained by numerical evaluation of the force. Various authors have proposed different computational methods, but all of them rely on force formulation. This paper presents an alternative energy approach in the equilibrium height calculation, which was settled on the minimum energy principle. The problem, involving a cylindrical magnet suspended in a closed cylindrical container full of magnetic ﬂuid, was considered in the study. The results accomplished by the proposed method were compared with those of the well-established surface integral method already veriﬁed by experiments. The difference in the results gained by both methods appears to be under 2.5%.


Introduction
In recent years, many researchers working on problems involving magnetic fluid were focused on the self-levitation phenomenon, which occurs when a permanent magnet (PM) is suspended in magnetic fluid (MF). This phenomenon was first reported by Rosensweig in 1966 [1,2]. A study of the forces acting on such a body has shown that, in addition to the Archimedes buoyancy, another type of buoyancy appears by virtue of the magnetic field. Such force is called a magnetic levitation force (MLF), also known as second-order buoyancy. This new kind of buoyancy is not limited only to the vertical direction, but acts as the restoring force for any displacement of the PM from the stable levitating position.
Because of this, a PM suspended in MF has found many practical applications, the most promising of which are dampers and accelerometers. The first damper of this kind, consisting of PM and MF, was the ferrofluid inertia damper proposed by Moskowitz et al. in 1978, designed to operate in the self-levitated condition [3]. Later on, the configuration was improved to achieve better performances, but the basic concept remained unchanged [4,5]. The passive ferrofluid dynamic vibration absorber presented in [6] is also a damping device based on the self-levitation phenomenon. Such an absorber is used to suppress small amplitude (<1 mm) and low-frequency vibrations (<1 Hz). Some modifications to the beforementioned absorber have been disclosed in [7]. For more recent theoretical and experimental studies on the ferrofluid dynamic vibration absorbers, the readers may refer to [8,9].
Another interesting and promising application that works on the self-levitation principle is the magnetic fluid-based accelerometer. The typical design of such a device, presented in [10], is similar to the one introduced with the dampers, but here the inductive coil is mounted over the housing to detect the PM movement. In the past few years, the magnetic liquid-based accelerometers have been researched intensively due to their relatively simple and reliable design [11,12]. Almost all the relevant publications in the field were provided by the research groups from China, where the topic of second-order buoyancy has been addressed both theoretically and experimentally. Yang et al. [13] calculated the magnetic self-levitation force of a cylindrical PM with the use of the magnetic charge image method and current image method, respectively. Results have shown a large difference of up to 62% between the calculated and measured results. On the other hand, by employing FEM on the problem the error did not exceed 9%. The same authors continued with research on the radial MLF in the following years [8,14]. The current image method and FEM were employed successfully in the optimization of a magnetic fluid accelerometer [11]. Based on the magnetic field solution obtained by the commercial FEM software, three effective and practical computational methods were proposed in [12]: the surface integral method (SIM), the magnetic force method, and equivalent magnetic force method. The results achieved by the aforementioned methods matched the experimental results within 7%. The experimental setup for measuring MLF involved a non-magnetic rod attached to the top of the PM which influenced the force acting on the PM. The detailed model based on the SIM and the influence of the measuring rod were studied in [15,16] There were also some research studies devoted to the development of micro magnetic fluid devices for driving micro machines. Here, the PM is coated with MF, which acts as a lubricant. The movement of the PM can be achieved in an alternating magnetic field [17,18].
This article uses two different approaches in calculating the equilibria states of the permanent magnet (PM) immersed in the magnetic fluid (MF). The first approach is based on the force formulation, where the stable magnetic state is calculated through the Surface Integral Method (SIM). This method is well covered, theoretically and experimentally, by various authors [12,15,16,19]. The main reason it appears in the current study is to use it as a benchmark for the proposed method, based on the energy formulation.
To the best of the authors' knowledge, the self-levitation equilibria of the PM have not been studied by virtue of the energy. To fill this gap in the field, and to illuminate the addressed topic from a different standpoint, a mathematical model grounded on the energy principle is presented in this article. Such a formulation is essential, since numerical issues arise when the SIM is implemented on the first-order finite element mesh [20]. Not going in deep with the error analysis associated with the numerical evaluation of the stress tensor, only a summary on the subject is given here. The mesh size on the contour path should be as fine as possible; the integration should be performed in the free space (air), or at least in a region with a constant permeability; the contour should not be defined at the interface between two different materials, but always a few elements away from any interfaces or boundaries [20,21]. In practical applications, whether the user is coding or just using the FEM-based software, considering the enlisted requirements could represent quite demanding work. Anyhow, when the self-levitation equilibrium is tackled with the energy formulation, the issues mentioned above could be avoided successfully.

Materials and Methods
The content of this chapter is devoted to the problem definition, and to introduce the reader to the basic, but still comprehensive, theoretical background. Even though the quite extensive mathematical employment used here could be omitted, as the majority of it can be found in [16], the authors have kept it here on purpose; in the first place, to point out some crucial differences between the existing approaches reported by other authors.

Problem Definition
Analysis of the self-levitation phenomena associated with a PM immersed in MF depends on the engaged geometry and material properties of the solids and liquids involved [11]. This study is focused on the simple cylindrical geometry, where both the PM and the container with MF are cylindrical. Onward, the PM is placed in a closed and fully filled container in such way that both axes coincide. The basic deployment is shown in Figure 1.

Surface Integral Method (SIM)
The mathematical expression of the magnetic fluid levitation force (MFLF) F l which acts on the PM immersed in MF could be derived from magnetic fluid stress tensor T m proposed by Rosensweig in 1966 [22]. The origin of the force is in the interaction between a permanent magnet's field and the magnetic liquid's magnetization. It drags a permanent magnet into the equilibrium state, where the total energy reaches its minimum.
Here, S p is the surface of the PM; n p is a unit vector normal to the PM's surface, pointing outward in the MF (see Figure 2); da is a surface element; p is the pressure in the MF in the absence of a magnetic field; v is the volume of MF per unit mass i.e., (1/ρ) where ρ stands for the MF's density; µ 0 is the permeability of the free space; M is the magnetization of the MF; B and B are the magnetic flux density in the vector notation and its absolute value, respectively; H and H are the magnetic field intensity in the vector and its corresponding absolute value, respectively; I is an identity matrix; and T is the absolute temperature. The product BH represents a tensor. All the physical quantities stated in Equation (1) and Equation (2) are, in general, functions of position, e.g., B = B(x, y, z), H = H(x, y, z), p = p(x, y, z), M = M(x, y, z), v = v(x, y, z), etc., but a shorter notation was used here for convenience.
Executing a derivation in the square brackets of the above expression at constant H and T gives When the fluid is at rest and fully enclosed within the container (there are no free surfaces), the pressure p represents the hydrostatic pressure caused by the gravity. To emphasize this, the subscript for the gravity is added to it, i.e., p = p g . The notation of Equation (3) is usually written in the following form [22], where the fluid-magnetic pressure p m and the magnetostrictive pressure p s are defined as In most of the practical cases where the fluid is considered incompressible, the latter term can be omitted from the further analysis, i.e., p s ≈ 0 [22]. Considering the MF equilibrium state with respect to the representation shown in Figure 2, the following equation can be given.
Here, S c is the surface of the container; n c is a unit vector normal to the surface S c , pointing inward to the MF (see Figure 2); V m f is the volume of the MF; g is the acceleration vector of gravity; and dV is a volume element. Equation (7) gives a balance of the forces in the MF; here, the term on the left side of the equation represents the gravity force of the MF, and the right side stands for the buoyancy force, which consists of magnetic and Archimedes components. The surface exposed to the MF is here indicated with S.
Expressing Equation (1) with Equation (7) results in a more propriate formulation of magnetic fluid levitation force F l [16].
Considering Equations (4)- (6) in Equation (8) a further development follows: By looking at Equation (9), it can easily be noticed that the first two terms on the right side correspond to the Archimedes buoyancy force F la . The others terms declare the magnetic contribution to the magnetic fluid levitation force F lm , i.e., F l = F la + F lm .
In the relevant literature, the force F lm also refers to the magnetic levitation force (MLF), the opposite value of which is known as second-order buoyancy F I I . In other words, [15,16,19].
The magnetic part of the second-order buoyancy force, or MLF (F lm ), can be calculated by integrating solely over the container's surface S c . The numerical evaluation of the second-order buoyancy force in this article, as well as in the research studies conducted in [15], is realized based on the statement given in Equation (12). The conditions are given illustratively in Figure 3; the p m acts in favor of the PM's stable position. . Visual representation of how magnetic fluid pressure contributes to the self-levitation mechanism: (a) In this position, a magnetic field is much higher at the bottom of the container than it is at the top of it, thus, the magnetic fluid pressure is also higher at the bottom, which results in PM upward movement; (b) In the stable levitation position, or in equilibrium height, the magnetic fluid pressure is balanced within the container; (c) If the PM tends to move up somehow, the pressure on the top increases due to the higher magnetic field, and the PM starts to move toward the stable position.
Maybe it is not evident at first, but Equation (12) gives a very useful mathematical tool. Namely, while all the physical significance of the interest, i.e., the force F I I , is reflected through the integration over the surface bounding the MF, despite changes inside the boundary. The following illustrative examples shown in Figure 4 should give a transparent understanding of Equation (12). In the first case, labeled by (a) in the Figure, where only a PM is floating inside the liquid, force F I I is obtained by integrating along the four Sections 1-2, 2-3, 3-4, and 4-1. Now, suppose some other object besides PM floats within the MF, as shown in Figure 4b. While the surface of the container is not affected by the additional object, the calculation of force F I I follows the same procedure. If the setup is given as in Figure 4c, where a non-magnetic rod is attached to the PM and passes through the container, then the integration over the rod's surface S r in the container's hole should be included in the force F I I calculation, due to the difference in the magnetic properties of the rod and the MF. The integration pertaining to the MF is performed through the paths 1-5, 6-2, 2-3, 3-4 and 4-1, while Sections 5 and 6 passes through the non-magnetic rod. A peculiarly interesting case relevant for an experimental verification or a practical application is the one marked by (c) in Figure 4. Here, the non-magnetic rod is attached to the PM at one side, while the other-upper-side goes to a dynamometer not shown in the figure. Because such an arrangement was used in the experiments conducted in [15], the same composition comes with the simulation model presented by the current study.

Energy Method
It was mentioned in the previous section that the static equilibrium is met when the energy of a physical system under consideration reaches its minimum. In classical physics, this is known as the principle of minimum energy. The principle has its background in thermodynamics, and could be applied prosperously to a wide range of problems involving magnetic fluids, such as magnetocaloric energy conversion and magnetic liquid free surface instabilities [22,23]. Not wishing to enter further into the topic, for the case study presented here, it suffices to consider the total energy as a sum of the gravitational and magnetostatic energy where W tot is the total energy; W g and W m are gravitational potential energy and magnetostatic energy, respectively.
Here, the gravitational energy W g is given as a function of a PM levitation height, as shown in Figure 5. Figure 5. Schematic representation of a magnetic field caused by the PM levitating in the MF.
A specific gravitational energy, or gravitational energy per unit of volume, can be stated as Here, ρ p and ρ m f are the density of the PM and MF, respectively; ∆ρ is the difference in material densities; V p is the PM volume; z 0 is the levitation height (see Figure 5); w g is the density of gravitational energy; and g is gravity acceleration. For simplicity, the MF is assumed as a homogeneous media, thus the parameters ρ p and ρ m f could be considered as constants. In reality, the concentration of magnetic nanoparticles is denser in the vicinity of the PM, while the magnetic nanoparticles are attracted toward the region with the higher magnetic field [11].
While the levitation height z 0 of the PM is subjected to change, the magnetic field will not be distributed symmetrically with respect to the z-axis, and, consequently, the magnetostatic energy density will vary with the position of the PM. To take into account the spatial variation of the energy as well as the MF nonlinear magnetization characteristic, the magnetostatic energy W m is written as a function of the density of magnetostatic energy w m (x, y, z) as follows, In the latter equation, B and B stand for the lower and upper integration limits, and represent the initial and the final fields established in the magnetizing process, respectively. The expression for the magnetostatic energy density given with Equation (17) holds in general for both the linear and nonlinear relations between quantities B(x, y, z) and H(x, y, z). With the linear magnetic conditions, this expression becomes w m (x, y, z) = BH/2.
The main idea is to find the PM levitation height z 0 at which the total energy announced by Equation (13) is minimal. The aforementioned task is not achieved easily without computer-aided numerical methods; after all, knowing the magnetic field distribution in the space is essential to the problem. However, still, with some simplifications applied to the model, the problem can be handled analytically to some extent. The following consideration of the equilibrium state of the PM immersed in an MF is addressed by qualitative rather than quantitative analysis.
Suppose the arrangement as shown in Figure 6, which is unbounded in both x and y directions. The MF with constant permeability µ occupies the space from the reference level at z = 0 to the height of the container at z = L. The space surrounding the MF is free space or air with the permeability µ 0 . Now, imagine having such a PM which has an infinitesimally small thickness and is extended extensively over the x − y plane at the height z = z 0 within the MF. As indicated in Figure 6, the computational domain is divided into four sections, 1 to 4. Section 1 corresponds to the region with the MF above the PM, Section 3 announces the free space above the PM, and a similar indication is carried out for the regions below the PM, where Sections 2 and 4 pertain to the MF and free space region below the PM, respectively. The magnetic field intensity and magnetic flux density at the PM's surface are H p and B p , respectively, which are specified by the PM's properties. Assuming the edges of the PM are distant from the segment of interest, then a z-component of a magnetic field will dominate near that fraction of the PM. To the first approximation, an exponential decrease in the magnetic field can be assumed with respect to the vertical distance from the piece of the PM in question. The rate of change in the field is determined by another PM's constant, α. Here, the levitation height of the PM is indicated with z 0 , and is taken as a parameter. The relation between B and H in the matter can be written in a particularly useful form when the medium has a constant permeability According to Figure 6, the magnetic field intensity in the MF's region above the PM in region 1, H 1 , is stated with Equation (20).
Applying Equation (19) to Equation (20) gives a magnetic flux density in region 1 At this point, it is appreciated to assert a total differential of B 1 (z) with respect to z, which appears in the expression for the density of magnetostatic energy in Equation (17) where B a corresponds to the magnetic flux density at the MF, i.e., the free space interface where z = L. Similarly, B p corresponds to the magnetic flux density at the PM's surface z = z 0 . w m1 is the density of magnetostatic energy in region 1. Similar expressions can be obtained for region 2 as follows: Here, B b corresponds to the magnetic flux density at the MF, i.e, the free space interface where z = 0. B p corresponds to the magnetic flux density at the PM's surface z = z 0 . w m2 (z) is the density of magnetostatic energy in region 2.
At the boundary between the MF and the free space (region 1 and region 3), the magnetic flux density remains unchanged, i.e., B 3 (z) = B 1 (z), Hence, the magnetic field intensity in the free space above the PM is described with the following expressions: where w m3 (z) stands for the magnetostatic energy density in region 3. The foregoing prosecution can be adapted similarly to the magnetic conditions in region 4.
Here, w m4 (z) corresponds to the density of the magnetostatic energy in region 4. Employment of the presented approach to the problem addressed in this article is indicated in Figure 7. Here, the column consisting of all four regions is virtually isolated from the rest of the space, and considered as an independent problem. The diameter D of the column equals the one of the PM, and also the other spatial extents should be selected so, to reflect the real situation. According to Figure 7, the magnetostatic energy accumulated in the column is the sum of the regions' energies, where S 0 is an area of the PM's base; a is an upper limit in region 1; c is a lower limit in region 4; L is the height of the container; and z 0 is the levitation height. For convenience, the dependence of w m on the vertical height z is not written explicitly in Equation (36). While the PM is abstracted to have the mass m p but no volume, the gravitational potential energy can then be written as The total energy is then, by Equation (13), the sum of magnetostatic and gravitational potential energy. The graphical representation of energies with respect to the levitation height z 0 is plotted in Figure 8. From the shape of the energy curves, the following deduction can be made. The minimum of the magnetostatic energy corresponds to the state where the second-order buoyancy force is equal to zero, which is in the center of the magnetic fluid region. Besides, the influence of the magnet's weight on the equilibrium height is evident in Figure 8. The stable height, or minimum in total energy, is moved toward lower values, which is somehow the expected effect. However, the aforementioned conclusion might be obtained intuitively as well, but here the effort was made for the reader's convenience. All the implementations of this section were referring to the basic model given in Figure 2, and could be extended to a model comprised of multiple bodies; in other words, the expression of total energy given by Equation (13) should be understood as the sum of all the gravitational potential and magnetostatic energies within the observed volume.

Simulation Setup
Based on the presented theory, the actual calculations were executed by the FEMbased computer program FEMM 4.2. Among others, the used FEM software is capable of handling two dimensional (2D) axisymmetric magnetostatic problems, and is a reasonable choice for launching such numerical effort. The dimensions and FEMM model of the problem in preprocessor mode are shown in Figure 9. The right side of the Figure shows the regions with finer mesh, and those regions pertain to the integration paths used with the SIM. The mesh density is chosen automatically by the FEMM 4.2 code, and it is constrained by the maximal mesh size defined by the user. Not shown in Figure 9, the computational domain is enclosed with a 100-mm radius of free space. The properties of the materials used with the calculation are listed in Table 1. The relative permeability µ r is a unitless number revealing how many times a magnetic material is more susceptible to magnetization of a free space, defined as µ r = µ/µ 0 . The initial susceptibility is defined as the ratio M/H at a low magnetic field intensity. The PM consists of a strong rare earth material that has quite linear B − H characteristics, with its slope very close to µ 0 , i.e., the relative permeability of such a magnet is very close to the unit. Strictly speaking, the actual value moves around 1.05, but here the idealized case was acquired to fit with the one stated in [15]. In particular, the values of the remanence B r and the coercivity H c of the magnet are given in Table 1.

Results and Discussion
Based on the presented numerical models introduced in the previous section, the simulations were performed for each of the proposed concepts. The concept employing the SIM method appears in two stages; firstly, the second-order buoyancy F I I is calculated with respect to the levitation height z 0, and, thereafter, the equilibrium height is determined as a point where the magnetic levitation force compensates the gravitational force of the floating PM. The alternative approach introduced here as a novelty used the principle of minimum energy as a criterion of stable levitation height.

Calculations Based on the Surface Integral Method (SIM)
The results of the first approach, in which the magnetic levitation force is calculated by applying the SIM, are presented in Figure 10. The measurements refer to the ones given in ref. [15] (p. 327- Figure 5b) for the case where the diameter of the rod is 2 mm. In the figure, the calculated values of F I I were evaluated at the same points as the measured values exposed in the aforementioned work of Yu et al. [15]. Numerical values and error analysis can be found in Table 2. Regardless of the deviations between measured and calculated values, the SIM can be considered as a suitable method when a numerical determination of the magnetic levitation force is in focus.
However, suppose that the influence of the rod attached to the top of the PM can somehow be eliminated from the measurements while it is a part of the measurement device which measures the force (dynamometer). Then, by knowing the second-order buoyancy force F I I , the equilibrium height or stabile levitation state of a permanent magnet is determined by finding an intersection point between the gravitational force of the PM in the MF, F g , and F I I . With respect to the calculated values, the equilibrium height of the PM for the case study corresponds to z 0 = 7.64 mm. The intersection is drawn in Figure 10 (the intersection between F g and F I I ). Gravitational force is calculated as the weight of the PM subtracted by the weight of the displaced fluid for a given example F g ≈ 0.094 N.   Figure 5b) for the case where the diameter of the rod is 2 mm.

Calculations Based on the Energy Method (EM)
Obeying the minimum energy principle, the equilibrium height of the levitation object can be found at the point where the total energy is minimal. Based on the theoretical conclusions brought in Section 2.1, the PM is expected to levitate in the lower half of the container. Magnetostatic energy calculations were executed in the FEM-based software FEMM 4.2. The properties of the model used with the calculations are given in Figure 9a and Table 1. The computational strategy for the problem was selected as follows. The total energy was evaluated at evenly spaced distances from the container's bottom in steps of 0.1 mm. The attained discrete energy values were thereafter approximated by the sixth degree polynomial curve. The equilibrium height, z 0 , was determined at the minimal energy value on the curve.
The results obtained by numerical analysis, which are shown in Figure 11, fit well with the theory at this point. The figure is drawn with two energy axes, due to the large difference between the magnetostatic W m and gravitational potential energy W g ; the left axis refers to the magnetostatic energy (in Joules), while the right applies to the gravitational potential energy (in milli-Joules). A careful analysis of the energy curve W tot shows that the levitation height coincides with the minimum of the total energy at z 0 = 7.58 mm. At first look, the result is in good agreement with the one attained by the SIM, which was verified experimentally. However, the results presented so far are based on a single case, and no general conclusion can be given on the subject. For that, a comprehensive parametric analysis is brought in the section that follows.

Parametric Analysis
To verify the proposed EM method, an additional parametric analysis of magnetic levitation was performed. The following quantities were selected as parameters, and varied with respect to the reference case. Five cases, marked A, B, C, D, and E, were included in the analysis; a single parameter was varied in each of the cases. The reference case was the one with the geometry introduced in Figure 9a and with the properties stated in Table 1 (i.e., r = 2 mm, B r = 1.24 T, µ r = 1.31, d 0 = 10 mm, h 0 = 20 mm, D = 20 mm, and L = 40 mm). The container's vol-ume and the PM's volume were taken as constants, and did not change with the variation of geometry. Due to the large difference in energy values, each of the energy curves in the figures were normalized with its maximum value and were stated on a per unit (p.u) scale. A minimum value on the energy curve is indicated with a dot. The error analysis was attained with help of the following formula ∆ (%) = z 0,SI M −z 0,EM z 0,SI M ·100. Here, ∆ is the relative error expressed in %, z 0,SI M is the equilibrium levitation height computed with SIM, and z 0,EM is the equilibrium levitation height computed with EM.

Case A (Variation of the Rod's Radius-r)
The impact of the rod's radius (r) on the PM's equilibrium height is presented in this sub-section. The parameter was varied from 1 mm to 3mm in steps of 1 mm. The results were obtained from Figures 12-15 and collected in Table 3. In this case, both methods can be compared with the experimental results deduced from [15]. Other relevant parameters to the analysis are given as B r = 1.24 T, µ r = 1.31, d 0 = 10 mm, h 0 = 20 mm, D = 20 mm, and L = 40 mm. The mesh parameters were the same as the ones stated in Table 1, and the mesh density was selected automatically by FEMM 4.2.  The force was calculated by the SIM. An intersection between the gravitational F g and mag-netic levitation force F I I is shown in the magnified window. Measured data were acquired from ref. [15]. Figure 13. The measured and calculated curves of magnetic levitation force F I I for the case r = 2 mm. The force was calculated by the SIM. An intersection between the gravitational F g and mag-netic levitation force F I I is shown in the magnified window. Measured data were acquired from ref. [15].

Figure 14.
The measured and calculated curves of magnetic levitation force F I I for the case r = 3 mm. The force was calculated by the SIM. An intersection between the gravitational F g and mag-netic levitation force F I I is shown in the magnified window. Measured data were acquired from ref. [15].

Case B (Variation of the MF's Relative Permeability-µ r )
This subsection exposed the impact of the MF's relative permeability (µ r ) on the PM's equilibrium height z 0 . Three values of µ r were considered in the calculations: µ r = 1.16, µ r = 1.31, and µ r = 1.46, respectively. The results are given in Table 4. The solutions are presented graphically in Figures 16 and 17. Other relevant parameters used with the analysis were r = 2 mm, B r = 1.24 T, d 0 = 10 mm, h 0 = 20 mm, D = 20 mm, and L = 40 mm. The mesh parameters were the same as the ones stated in Table 1.   Variation of the PM's remanent flux density (B r ) and its influence on the PM's equilibrium height z 0 is studied here. The values of B r were set as follows: B r = 0.7 T, B r = 1 T, and B r = 1.24 T. Other relevant parameters used with the analysis were r = 2 mm, µ r = 1.31, d 0 = 10 mm, h 0 = 20 mm, D = 20 mm, and L = 40 mm. The mesh parameters were the same as in the ones stated in Table 1. The results of the analysis are collected in Table 5 and presented in Figures 18 and 19.   The results presented in Table 6 and Figures 20 and 21 demonstrate the influence of the PM's geometry proportions on its equilibrium height z 0 . The analysis was performed with three different values of d 0 : 9 mm, 10 mm, and 11 mm. The geometry of the PM was conditioned with the assumption of constant PM volume. Other relevant parameters used with the analysis were r = 2 mm, µ r = 1.31, B r = 1.24 T, D = 20 mm, and L = 40 mm. The mesh parameters were the same as those stated in Table 1.   The geometry of the container was varied similarly as in the previous case (case D). The container's diameter (D) was varied at a constant container volume. The results presented in Table 7 and Figures 22 and 23 show how the container's proportions impact the PM's equilibrium height z 0 . Other relevant parameters used with the analysis were r = 2 mm, µ r = 1.31, B r = 1.24 T, d 0 = 10 mm, and h 0 = 20 mm. The mesh parameters were the same as those stated in Table 1.

Conclusions
In the current study, the numerical investigation of the equilibrium height of a PM immersed in MF was carried out in two different concepts: the surface integral method (SIM) and the employment of the energy method (EM). The latter approach brings novelty to the field. To verify the proposed method, a comprehensive parametric analysis was carried out on the subject. The analysis encompassed five cases, where the properties of both MF and the PM were varied. For each of the five cases, the PM's equilibrium height was calculated with both methods, i.e., SIM and EM. In addition to the numerical results and analysis, the reader is also provided with the essential theoretical background on which the energy concept is based.
The following can be summarized, focusing on the EM, presented as an alternative and novel method to calculate the equilibrium levitation height of a PM drowned in MF. In both concepts based on EM and SIM, magnetic field distributions in space are essential data which are nowadays, almost without exception, obtained numerically, most commonly with FEM. From that point of view, both methods are practically equivalent in the computational time consumption. The main advantage of the EM over the SIM is that the adoption of the Maxwell stress tensor surface integral can be avoided. As mentioned in the Introduction, a numerical problem arises when evaluating this integral on a finite element mesh made of first-order triangles [20].
The following remarks can be stated based on the parametric analysis.
• By increasing the radius of the measuring rod (r), the PM's levitation height will tend to reach a higher level. The reason is evident; while the rod usually consists of non-magnetic material, the fluid-magnetic pressure at the PM's top face is going to produce a lower force than at the PM's bottom face. The difference in this pressure acts upward, hence, the magnet will levitate at a higher level. For this case, the error between SIM and EM was less than 1%.

•
The magnetic properties of the MF are reflected in the value of relative permeability µ r which is, here, assumed as a constant. Truly speaking, this is only a rough approximation. Namely, in real magnetic fluids, the magnetizing characteristic are non-linear, and the concentration of the particles also decreases with the distance from the magnet. However, even with this simplification, the conclusion can be declared that the PM immersed in a magnetic fluid with lower permeability will experience a weaker magnetic buoyancy force. Consequently, the levitation height will be lower to some extent.

•
The same can be derived by considering the problem in the light of the energy. In highly permeable MF, the PM would rather be centered in the fluid, because its magnetostatic energy is minimal there. By decreasing the MF's permeability the magnetostatic energy increases and its variation with respect to the levitation height becomes smaller, thus, the influence of the gravitational energy comes to the fore, and moves the minimum of the total energy curve towards a lower value. The results' comparison obtained by the SIM and EM for the case shows that the error between the methods is less than 4%. The results fit better at higher values of the MF's permeabilities. • A remanent magnetic flux density B r is a property of a PM, and from the presented results, it is evident that higher B r leads to a higher levitation height. The remanence B r is related to the origin of the magnetic force in the MF, and, therefore, it is reasonable to expect such an effect. Again, a good agreement was achieved between the methods. The maximum error value was about 5%. • Variation in the PM's geometry shows that the levitation height is proportional to the PM's diameter defining the top and the bottom surface area. The theoretical explanation is straightforward. The magnetic levitation force appears mainly as a consequence of a fluid-magnetic pressure acting on those surfaces. An increase in the surface area will result in a rise in the magnetic levitation force. The numerical results attained by the SIM and EM showed the same tendency, with the maximum error of 12.2%. • Change in the container's geometry acts oppositely to the change in the PM's geometry as discussed previously. Namely, in this case, the magnetic force acting on the PM reduces proportionally to the ratio of the diameters of the PM and container. The error in the analysis between the methods was less than 1%.
Regarding the presented numerical results, both methods, the SIM and EM, predicted similar results. The maximum deviation between the results obtained by the methods was 12.2%, but for most of the cases, the expected error should not exceed a few percent. However, the SIM for Case A was verified experimentally, where the maximum error in the lower half of the container appeared to be 15% [15]. The authors of the text believe that the error can be reduced significantly by an improved MF model that apprehends the magnetic nonlinearity and inhomogeneity of the particle distribution. The researchers in the field may find this a topic to focus on.
In the end, a final remark on the subject can be made as follows: when the numerical evaluation of the PM levitation height is in question, satisfactory results can be attained by the proposed EM approach.

Conflicts of Interest:
The authors declare no conflict of interest.