- Email: [email protected]
Numerical simulation of normal-field instability in the static and dynamic case
ARTICLE IN PRESS
Journal of Magnetism and Magnetic Materials 289 (2005) 346–349 www.elsevier.com/locate/jmmm
Numerical simulation of normal-field instability in the static and dynamic case Gunar Matthiesa, Lutz Tobiskab, a
Ruhr-Universita¨t Bochum, Universita¨tsstraX e 150, D-44780 Bochum, Germany Otto-von-Guericke-Universita¨t, Postfach 4120, D-39016 Magdeburg, Germany
b
Available online 28 November 2004
Abstract The normal-field instability is studied numerically by a finite element approach based on the 3D equations for the magnetostatic potential and the 3D Navier–Stokes equations for the velocity and pressure fields of the magnetic fluid. The finite element method is able to determine the critical value of the magnetic field intensity for the onset of the instability and to describe the transition from the undisturbed surface to the fully developed pattern in the dynamic case. r 2004 Elsevier B.V. All rights reserved. PACS: 75.50.Mn; 02.60.Cb; 47.20.Ma Keywords: Normal-field instability; Ferrofluid; Finite element solution
1. Introduction In the normal-field instability, a uniform magnetic field applied perpendicular to a layer of ferrofluid produces a spontaneous surface deformation when the field exceeds a critical value. This stationary wave pattern was first observed by Cowley and Rosensweig [1]. It is caused by the interaction of gravitational, capillary and magnetic forces. Additionally, in the transition from the undisturbed surface to the fully developed pattern, hydrodynamic forces have to be taken into consideration. Stability analysis of this phenomenon is often restricted to determine the critical magnetic field intensity in the case of small perturbations [2]. A numerical simulation based on the nonlinear system of coupled Corresponding author. Tel.: +49 391 6718650;
fax: +49 391 6718073. E-mail address: [email protected] (L. Tobiska).
Maxwell and Navier–Stokes equations together with the Young–Laplace equation, representing the force balance at the unknown free surface, opens the way for a deeper theoretical insight. Our main objective is to describe different numerical algorithms to handle both the static and the dynamic case. Simplified numerical models for the static case have been already studied in [3,4].
2. Governing equations The Maxwell equations for a non-conducting fluid are given by curl H ¼ 0;
div B ¼ 0;
in O;
with the magnetic field strength H and the magnetic induction B satisfying the constitutive relation B ¼ m0 ðM þ HÞ in OF ðtÞ; B ¼ m0 H outside OF ðtÞ;
0304-8853/$ - see front matter r 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.jmmm.2004.11.098
ARTICLE IN PRESS G. Matthies, L. Tobiska / Journal of Magnetism and Magnetic Materials 289 (2005) 346–349
where M is the magnetization, m0 ¼ 4p 107 Vs=Am is the permeability constant, OF ðtÞ is the subdomain occupied by the ferrofluid at time t. The magnetization M is assumed to be parallel to the magnetic field and to follow the Langevin law, i.e., 1 H M ¼ M S cothðgjH jÞ gjHj jHj with the saturation magnetization M S ; the Langevin parameter g ¼ 3w0 =M S ; and the initial susceptibility w0 : The hydrodynamic properties are described by the timedependent incompressible Navier–Stokes equations rðut þ ðu rÞuÞ ¼ div sðu; p; HÞ þ f
in OF ðtÞ;
div u ¼ 0
in OF ðtÞ;
where sðu; p; HÞ is the stress tensor with s ¼ Zðru þ ruT Þ ðp þ
m0 jH j2 Þ þ BH: 2
Here, r denotes the density, Z the dynamic viscosity, u the velocity, and p the pressure. Finally, we have the force balance and the kinematic condition at the free surface ½sn ¼ aKn;
u n ¼ VG
on GF ðtÞ;
with the sum of the principal curvatures K, the coefficient of surface tension a; the unit normal vector n; and the normal velocity V G of the free surface GF :
3. Problem splitting and discretisation The coupled problem is iteratively split into subproblems which depend on whether the static or dynamic case is considered. In the static case we have u ¼ 0 and the pressure p can be determined directly from the Navier–Stokes equations. Thus, in this case we solve a 3D magnetostatic problem inside and outside of the fluid for a given OF and update OF by solving the 2D Young–Laplace equation for the graph representing GF : For the dynamic case we incorporate the force balance on GF ðtÞ; i.e. the Young–Laplace equation, in the flow calculation which leads to the following iteration: Given OF ðtn Þ we solve first a 3D magnetostatic problem inside and outside of the fluid, then a 3D Navier–Stokes problem in a time-dependent domain by the ALE (Arbitrary-Lagrangian–Eulerian) method, and finally use the information of u n at GF ðtn Þ to find the position for GF ðtnþ1 Þ: For approximating the solutions in space finite element methods were applied. In the static case we used piecewise bilinears for representing the graph of GF and piecewise trilinears for the magnetostatic potential. In the dynamic case we used a fractional step W-scheme
347
of second order for the time discretisation of the Navier–Stokes equations and piecewise triquadratics for both the magnetostatic problem and the velocity while the pressure was approximated by discontinuous piecewise linear functions. The inf–sup stability condition of this second-order pair of finite elements has been proven in [5]. Of course we could also use continuous pressure approximations satisfying the inf–sup stability condition, like the Taylor–Hood element consisting of piecewise triquadratics for the velocity and continuous piecewise trilinears for the pressure. However, our experience is that the use of discontinuous pressure approximation lead to a better mass conservation.
4. Multilevel solution In the static case we have to solve two nonlinear systems of algebraic equations which correspond to the 3D problem for the magnetostatic potential and the 2D problem for the graph of the unknown free surface. In both cases a fixed-point iteration has been applied where the arising linear systems of equations were solved by a multilevel algorithm in each step of the iteration. In the 3D problem for the magnetostatic potential a geometric multilevel method has been applied based on a family of successively refined 3D hexahedral meshes. When solving the Young–Laplace equation over a 2D domain we have to find the projection of the surface mesh onto the plane. However, this results in quite a difficult data structure of the projected family of surface meshes. As an alternative we applied an algebraic (instead of a geometric) multilevel method for the solution of the linearised Young–Laplace equation. For this type of coupled 2D/3D approach special data structures had to be developed and were implemented in the programme package MooNMD [6]. In the dynamic case the situation is different. Instead of the Young–Laplace equation on a 2D domain we have to solve the Navier–Stokes equations on a 3D timedependent domain. This results in each time-step in a saddle point problem of the form ! ! u f AðuÞ BT ¼ p 0 B 0 for the velocity and pressure, respectively. This nonlinear system is linearised by a simple fixed point iteration replacing AðuÞ by Aðuold Þ where uold is a previous iterate of the velocity. The linear systems are again solved by a multilevel method. It should be mentioned that in contrast to the problem for the magnetostatic potential now a simple SSOR smoother cannot be used. We apply instead a Vanka-type smoother which essentially is a block Gauss–Seidel iteration where the blocks are built by the degrees of
ARTICLE IN PRESS 348
G. Matthies, L. Tobiska / Journal of Magnetism and Magnetic Materials 289 (2005) 346–349
freedom which are connected with one mesh cell. As already mentioned the information of normal velocity is then used to define a new position of the free surface in the next time-step. There are several mesh-update algorithms in the literature, for example, algebraic, elastic and complete remeshing algorithms are used. In our case a simple algebraic mesh-update has been used where the mesh-points within and outside the fluid are stretched according to the actual peak height.
5. Results of numerical simulations Our numerical simulation starts from a horizontally unbounded ferrofluid layer of a finite depth. For applied fields with a magnetic field strength slightly above the critical value a hexagonal pattern of peaks occurs [2,7]. The wavelength of this pattern at the onset of the instability can be obtained by a linear stability analysis [2]. We restrict our calculations to a hexagonal prism of finite height where the base was chosen such that it will contain exactly one peak. Furthermore, we assume that the upper and lower boundary have a sufficiently large distance from the ferrofluid such that there is no interaction between the interface and the magnetic field strength at the upper and lower boundary. We carried out numerical simulations for the ferrofluid EMG 909. From its physical parameters (density 1020 kg=m3 ; surface tension coefficient 0:0265 kg=s2 ; saturation magnetisation 16 kA/m, initial susceptibility 0.85) we obtained by formulas from Ref. [2] a critical magnetic field strength of 14.85 kA/m. Fig. 1 shows the peak height as a function of the applied field for different refinement levels of the mesh of the underlying finite element method. The given peak height is the height difference between the midpoint (highest point) and a vertex (lowest point) of the hexagon. We see that the qualitative behaviour is reproduced even on very coarse meshes. It is obvious
Fig. 2. Time-dependent peakheight for different applied fields switched on at t ¼ 0:
that one gets higher peaks on finer meshes. Furthermore, on the finest mesh we obtain numerically a value of the critical magnetic field strength which is very close to the theoretical value. In addition to the results of simulations of the static case one graph for simulations of the dynamic case is given in Fig. 1. The peak height after 2 s is shown. We see that the curve from the dynamic simulation is very close to the curve of static calculation on the finest mesh. The critical magnetic field strength is approximated very well. Fig. 2 presents the time-dependent peak height for different values of the applied magnetic field strength. The pictured peak height is the difference of the height of the midpoint of the hexagon and the height of one vertex of the hexagon. This is the reason for the negative heights which occur for small values of the magnetic field strength since the midpoint moves below the height of a vertex due to inertia. For all simulations of the dynamic case the applied magnetic field was switched on suddenly at time t ¼ 0: We see that one gets damped oscillations for values of the magnetic field strength above the critical value. The amplitude of the oscillations increases with the field strength. If the field strength is smaller than the critical value the peak height decreases to zero. Furthermore, the damping for higher values of the field strength is less than for lower values. In all cases the oscillations have almost disappeared for time larger than 1 s.
Acknowledgements The authors like to thank the German Research Foundation (DFG-FOR 301) for supporting this work.
References Fig. 1. Peak height depending on the applied field for different refinement levels of the mesh in the static case compared to the dynamic case on the finest mesh.
[1] M.D. Cowley, R.E. Rosensweig, J. Fluid Mech. 30 (1967) 671.
ARTICLE IN PRESS G. Matthies, L. Tobiska / Journal of Magnetism and Magnetic Materials 289 (2005) 346–349 [2] R.E. Rosensweig, Dover Publications, New York, 1997. [3] V.G. Bashtovoi, O.A. Lavrova, V.K. Polevikov, L. Tobiska, J. Magn. Magn. Mater. 252 (2002) 299. [4] O.A. Lavrova, G. Matthies, T. Mitkova, V.K. Polevikov, L. Tobiska, Lect. Notes Comput. Sci. Eng. 35 (2003) 160.
349
[5] G. Matthies, L. Tobiska, Computing 69 (2002) 119. [6] V. John, G. Matthies, Comput. Visual. Sci. 6 (2004) 163. [7] A. Gailitis, J. Fluid Mech. 82 (1977) 401.
Journal of Magnetism and Magnetic Materials 289 (2005) 346–349 www.elsevier.com/locate/jmmm
Numerical simulation of normal-field instability in the static and dynamic case Gunar Matthiesa, Lutz Tobiskab, a
Ruhr-Universita¨t Bochum, Universita¨tsstraX e 150, D-44780 Bochum, Germany Otto-von-Guericke-Universita¨t, Postfach 4120, D-39016 Magdeburg, Germany
b
Available online 28 November 2004
Abstract The normal-field instability is studied numerically by a finite element approach based on the 3D equations for the magnetostatic potential and the 3D Navier–Stokes equations for the velocity and pressure fields of the magnetic fluid. The finite element method is able to determine the critical value of the magnetic field intensity for the onset of the instability and to describe the transition from the undisturbed surface to the fully developed pattern in the dynamic case. r 2004 Elsevier B.V. All rights reserved. PACS: 75.50.Mn; 02.60.Cb; 47.20.Ma Keywords: Normal-field instability; Ferrofluid; Finite element solution
1. Introduction In the normal-field instability, a uniform magnetic field applied perpendicular to a layer of ferrofluid produces a spontaneous surface deformation when the field exceeds a critical value. This stationary wave pattern was first observed by Cowley and Rosensweig [1]. It is caused by the interaction of gravitational, capillary and magnetic forces. Additionally, in the transition from the undisturbed surface to the fully developed pattern, hydrodynamic forces have to be taken into consideration. Stability analysis of this phenomenon is often restricted to determine the critical magnetic field intensity in the case of small perturbations [2]. A numerical simulation based on the nonlinear system of coupled Corresponding author. Tel.: +49 391 6718650;
fax: +49 391 6718073. E-mail address: [email protected] (L. Tobiska).
Maxwell and Navier–Stokes equations together with the Young–Laplace equation, representing the force balance at the unknown free surface, opens the way for a deeper theoretical insight. Our main objective is to describe different numerical algorithms to handle both the static and the dynamic case. Simplified numerical models for the static case have been already studied in [3,4].
2. Governing equations The Maxwell equations for a non-conducting fluid are given by curl H ¼ 0;
div B ¼ 0;
in O;
with the magnetic field strength H and the magnetic induction B satisfying the constitutive relation B ¼ m0 ðM þ HÞ in OF ðtÞ; B ¼ m0 H outside OF ðtÞ;
0304-8853/$ - see front matter r 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.jmmm.2004.11.098
ARTICLE IN PRESS G. Matthies, L. Tobiska / Journal of Magnetism and Magnetic Materials 289 (2005) 346–349
where M is the magnetization, m0 ¼ 4p 107 Vs=Am is the permeability constant, OF ðtÞ is the subdomain occupied by the ferrofluid at time t. The magnetization M is assumed to be parallel to the magnetic field and to follow the Langevin law, i.e., 1 H M ¼ M S cothðgjH jÞ gjHj jHj with the saturation magnetization M S ; the Langevin parameter g ¼ 3w0 =M S ; and the initial susceptibility w0 : The hydrodynamic properties are described by the timedependent incompressible Navier–Stokes equations rðut þ ðu rÞuÞ ¼ div sðu; p; HÞ þ f
in OF ðtÞ;
div u ¼ 0
in OF ðtÞ;
where sðu; p; HÞ is the stress tensor with s ¼ Zðru þ ruT Þ ðp þ
m0 jH j2 Þ þ BH: 2
Here, r denotes the density, Z the dynamic viscosity, u the velocity, and p the pressure. Finally, we have the force balance and the kinematic condition at the free surface ½sn ¼ aKn;
u n ¼ VG
on GF ðtÞ;
with the sum of the principal curvatures K, the coefficient of surface tension a; the unit normal vector n; and the normal velocity V G of the free surface GF :
3. Problem splitting and discretisation The coupled problem is iteratively split into subproblems which depend on whether the static or dynamic case is considered. In the static case we have u ¼ 0 and the pressure p can be determined directly from the Navier–Stokes equations. Thus, in this case we solve a 3D magnetostatic problem inside and outside of the fluid for a given OF and update OF by solving the 2D Young–Laplace equation for the graph representing GF : For the dynamic case we incorporate the force balance on GF ðtÞ; i.e. the Young–Laplace equation, in the flow calculation which leads to the following iteration: Given OF ðtn Þ we solve first a 3D magnetostatic problem inside and outside of the fluid, then a 3D Navier–Stokes problem in a time-dependent domain by the ALE (Arbitrary-Lagrangian–Eulerian) method, and finally use the information of u n at GF ðtn Þ to find the position for GF ðtnþ1 Þ: For approximating the solutions in space finite element methods were applied. In the static case we used piecewise bilinears for representing the graph of GF and piecewise trilinears for the magnetostatic potential. In the dynamic case we used a fractional step W-scheme
347
of second order for the time discretisation of the Navier–Stokes equations and piecewise triquadratics for both the magnetostatic problem and the velocity while the pressure was approximated by discontinuous piecewise linear functions. The inf–sup stability condition of this second-order pair of finite elements has been proven in [5]. Of course we could also use continuous pressure approximations satisfying the inf–sup stability condition, like the Taylor–Hood element consisting of piecewise triquadratics for the velocity and continuous piecewise trilinears for the pressure. However, our experience is that the use of discontinuous pressure approximation lead to a better mass conservation.
4. Multilevel solution In the static case we have to solve two nonlinear systems of algebraic equations which correspond to the 3D problem for the magnetostatic potential and the 2D problem for the graph of the unknown free surface. In both cases a fixed-point iteration has been applied where the arising linear systems of equations were solved by a multilevel algorithm in each step of the iteration. In the 3D problem for the magnetostatic potential a geometric multilevel method has been applied based on a family of successively refined 3D hexahedral meshes. When solving the Young–Laplace equation over a 2D domain we have to find the projection of the surface mesh onto the plane. However, this results in quite a difficult data structure of the projected family of surface meshes. As an alternative we applied an algebraic (instead of a geometric) multilevel method for the solution of the linearised Young–Laplace equation. For this type of coupled 2D/3D approach special data structures had to be developed and were implemented in the programme package MooNMD [6]. In the dynamic case the situation is different. Instead of the Young–Laplace equation on a 2D domain we have to solve the Navier–Stokes equations on a 3D timedependent domain. This results in each time-step in a saddle point problem of the form ! ! u f AðuÞ BT ¼ p 0 B 0 for the velocity and pressure, respectively. This nonlinear system is linearised by a simple fixed point iteration replacing AðuÞ by Aðuold Þ where uold is a previous iterate of the velocity. The linear systems are again solved by a multilevel method. It should be mentioned that in contrast to the problem for the magnetostatic potential now a simple SSOR smoother cannot be used. We apply instead a Vanka-type smoother which essentially is a block Gauss–Seidel iteration where the blocks are built by the degrees of
ARTICLE IN PRESS 348
G. Matthies, L. Tobiska / Journal of Magnetism and Magnetic Materials 289 (2005) 346–349
freedom which are connected with one mesh cell. As already mentioned the information of normal velocity is then used to define a new position of the free surface in the next time-step. There are several mesh-update algorithms in the literature, for example, algebraic, elastic and complete remeshing algorithms are used. In our case a simple algebraic mesh-update has been used where the mesh-points within and outside the fluid are stretched according to the actual peak height.
5. Results of numerical simulations Our numerical simulation starts from a horizontally unbounded ferrofluid layer of a finite depth. For applied fields with a magnetic field strength slightly above the critical value a hexagonal pattern of peaks occurs [2,7]. The wavelength of this pattern at the onset of the instability can be obtained by a linear stability analysis [2]. We restrict our calculations to a hexagonal prism of finite height where the base was chosen such that it will contain exactly one peak. Furthermore, we assume that the upper and lower boundary have a sufficiently large distance from the ferrofluid such that there is no interaction between the interface and the magnetic field strength at the upper and lower boundary. We carried out numerical simulations for the ferrofluid EMG 909. From its physical parameters (density 1020 kg=m3 ; surface tension coefficient 0:0265 kg=s2 ; saturation magnetisation 16 kA/m, initial susceptibility 0.85) we obtained by formulas from Ref. [2] a critical magnetic field strength of 14.85 kA/m. Fig. 1 shows the peak height as a function of the applied field for different refinement levels of the mesh of the underlying finite element method. The given peak height is the height difference between the midpoint (highest point) and a vertex (lowest point) of the hexagon. We see that the qualitative behaviour is reproduced even on very coarse meshes. It is obvious
Fig. 2. Time-dependent peakheight for different applied fields switched on at t ¼ 0:
that one gets higher peaks on finer meshes. Furthermore, on the finest mesh we obtain numerically a value of the critical magnetic field strength which is very close to the theoretical value. In addition to the results of simulations of the static case one graph for simulations of the dynamic case is given in Fig. 1. The peak height after 2 s is shown. We see that the curve from the dynamic simulation is very close to the curve of static calculation on the finest mesh. The critical magnetic field strength is approximated very well. Fig. 2 presents the time-dependent peak height for different values of the applied magnetic field strength. The pictured peak height is the difference of the height of the midpoint of the hexagon and the height of one vertex of the hexagon. This is the reason for the negative heights which occur for small values of the magnetic field strength since the midpoint moves below the height of a vertex due to inertia. For all simulations of the dynamic case the applied magnetic field was switched on suddenly at time t ¼ 0: We see that one gets damped oscillations for values of the magnetic field strength above the critical value. The amplitude of the oscillations increases with the field strength. If the field strength is smaller than the critical value the peak height decreases to zero. Furthermore, the damping for higher values of the field strength is less than for lower values. In all cases the oscillations have almost disappeared for time larger than 1 s.
Acknowledgements The authors like to thank the German Research Foundation (DFG-FOR 301) for supporting this work.
References Fig. 1. Peak height depending on the applied field for different refinement levels of the mesh in the static case compared to the dynamic case on the finest mesh.
[1] M.D. Cowley, R.E. Rosensweig, J. Fluid Mech. 30 (1967) 671.
ARTICLE IN PRESS G. Matthies, L. Tobiska / Journal of Magnetism and Magnetic Materials 289 (2005) 346–349 [2] R.E. Rosensweig, Dover Publications, New York, 1997. [3] V.G. Bashtovoi, O.A. Lavrova, V.K. Polevikov, L. Tobiska, J. Magn. Magn. Mater. 252 (2002) 299. [4] O.A. Lavrova, G. Matthies, T. Mitkova, V.K. Polevikov, L. Tobiska, Lect. Notes Comput. Sci. Eng. 35 (2003) 160.
349
[5] G. Matthies, L. Tobiska, Computing 69 (2002) 119. [6] V. John, G. Matthies, Comput. Visual. Sci. 6 (2004) 163. [7] A. Gailitis, J. Fluid Mech. 82 (1977) 401.