Harvesting energy from the sloshing motion of ferrofluids in an externally excited container: Analytical modeling and experimental validation

Vibratory energy harvesting is a process by which, otherwise, wasted mechanical energy is transformed into useful electricity via the ability of active materials and electromechanical coupling mechanisms to generate an electric potential in response to mechanical stimuli.^1–5 This novel power generation technology has been recently utilized to design scalable power generators which provide low-power levels to maintain remote sensors and spatially distributed sensor networks that are, otherwise, hard to access and maintain. Further evolutions in this field of technology has the potential of lowering our dependence on batteries, which are known to have limited storage capacity, low-energy density, and require regular replacement and/or recharging.^6–10 

Typically, vibratory energy harvesters incorporate a solid energy transduction element. For instance, piezoelectric and magnetostrictive vibratory energy harvesters transform the strain produced in a solid piezoelectric beam or a magnetostrictive rod into electric charge or magnetic.^11–13 Electromagnetic energy harvesters transform mechanical motions into electricity by setting a solid magnet in motion relative to a stationary coil.^6,7,14–16 The solid nature of the transduction mechanism in typical vibratory energy harvesters can place limitations on their capabilities especially in terms of conformability to different shapes and sensitivity to external excitations. This has motivated the authors to explore liquid-state materials to transduce mechanical motions directly into electricity.^17–19 

To this end, this paper exploits liquid-state materials, namely, ferrofluids, as the transduction element in electromagnetic vibratory energy harvesters. Ferrofluids consist of stable ferrous nanoparticles in colloidal suspension forming nanoscale permanent magnetic dipoles.²⁰ The magnetic dipoles are coated with a surfactant to prevent conglomeration, such that, in the absence of an external magnetic field, the magnetic dipoles are randomly oriented in the carrier fluid as shown in Fig. 1(a). When an external field is applied, the dipoles rotate in the direction of the magnetic field producing a net magnetic moment as shown in Fig. 1(b). When a container carrying the magnetized ferrofluid is subjected to external mechanical stimuli, the sloshing motion of the magnetized ferrofluid generates a time-varying magnetic flux, which can be used to induce an electromotive force in a coil placed adjacent to the container. This process generates an electric current in the coil, and therewith, can be used to transduce external vibrations into electric energy as shown in Fig. 1(c).

In addition to conformability, a key advantage of this approach lies in the fact that the fluid column has a large number of closely spaced modal frequencies. This allows the fluid to respond to a wide range of excitation frequencies resulting in a broadband response. The presence of a large number of closely spaced modal frequencies also facilitates the activation of nonlinear modal interactions between the different vibration modes.²¹ This can lead to energy exchange among the interacting modes resulting in large-amplitude responses over a wide range of frequencies which broadens the voltage-frequency response curves even further.

In 2012, the authors introduced the concept of vibratory energy harvesting via the sloshing motion of ferrofluids.^18,19 An experimental investigation performed using a cylindrical container was presented to illustrate the feasibility of the proposed concept. Nevertheless, results also demonstrated that the bandwidth of the harvester is very small. This investigation was followed by an experimental case study performed by Viswas et al.²² in 2013 on a similar system. Subsequently, Dae Woong et al.²³ and Kim²⁴ also adopted the same concept and showed that the change of magnetic flux can be improved by utilizing a ferromagnetic yoke.

In a recent paper,¹⁹ we presented experimental results illustrating that the container carrying the ferrofluid can be designed with specific dimensions and fluid heights that make the modal frequencies of the fluid column nearly commensurate. This serves to activate a nonlinear energy transfer mechanism between the commensurate modes further improving the steady-state bandwidth of the harvester. In order to understand the influence of the important design parameters on the output power in the vicinity of the internal energy pump, this article formulates an analytical model which describes the electro-magnetohydrodynamics of the harvester (Section II). The availability of this model is essential to understand the role of the design parameters on the response and to optimize the power with respect to these parameters. This will help further increase the power density of the device and cannot be easily done by carrying a large number of experiments.

The presented model will be able to capture, with reasonable accuracy, the magnetohydrodynamic behavior which deals with the motion of the magnetized fluid in the container, and the electromagnetic induction which deals with characterizing how the fluid motion is induced into electrical energy. Using the developed model, the article first investigates variation of the modal frequencies with the height-to-width ratio of the fluid column and identifies the regions where internal resonances can be activated (Section III). Using the method of multiple scales, the article presents an approximate analytical solution of the harvester’s response to a primary resonance excitation of the first mode in the presence of a two-to-one internal resonance between the first and second modes (Section IV). An experimental setup is also constructed and used to validate the resulting model (Section V). Finally, conclusions regarding the validity of the model are presented (Section VI).

We consider the two-dimensional finite-amplitude sloshing dynamics of an irrotational, incompressible ferrofluid in a rectangular container of width L. The ferrofluid whose density is denoted by ρ is assumed to be of height h. Our goal is to characterize the dependence of the output voltage of the harvester on the design parameters. To this end, we consider the system shown in Fig. 2 with a rotating coordinate (⁠$x̄,z̄$⁠) located at point O. As the fluid starts to move due to external perturbations, surface waves of height, $η̄(x̄,t̄)$⁠, arise. The equations and associated boundary conditions governing the two-dimensional motion of the ferrofluid can be written as follows:

Here, the subscripts denote partial derivatives with respect to the independent variables.

Equation (1a) invokes the irrotationality assumption on the velocity field by expressing the two-dimensional velocity field, u, as the gradient of a scalar potential, $ϕ̄$⁠. Equation (1b) is a consequence of the incompressibility assumption for which the continuity equation requires the Laplacian of the velocity potential, $ϕ̄$⁠, to vanish. Equation (1c) states that the velocity normal to the sidewalls and bottom wall vanishes. Equation (1d) is the kinematic boundary condition at the surface which states that the velocity of a fluid particle on the surface must be equal to the velocity of the surface itself. Finally, Equation (1e) represents the dynamic boundary condition at the surface obtained by enforcing the unsteady Bernoulli’s equation.

In Equation (1e), the first term accounts for the unsteadiness in the velocity field; the second term represents the kinetic energy of the fluid; the third term accounts for the potential energy; where g is the gravitational acceleration; the fourth term represents the force due to surface tension, where σ represents the surface tension coefficient; and the fifth term represents forces exerted on the surface due to a harmonic base acceleration in the $x̄$ direction. Here, X₀ and ω are the amplitude and frequency of excitation, respectively. Finally, the sixth term represents the magnetic moment exerted on the surface, where $M̄$ and $H̄$ are the magnetization and magnetic field, respectively; μ₀ represents the permeability of vacuum, ρ is the fluid density, and $C̄$ is a constant.

Since the influence of surface tension on the nonlinear sloshing waves is considered, a contact angle should be enforced at the sidewalls. However, unless the contact angle is significantly different from 90°, i.e., horizontal surface, the contact angle has a very little influence on the dynamics of the surface. As such, to facilitate the analytical analysis, it is assumed that the contact angle is 90° and that the contact line is free to slip on the container’s surface.²⁵ Furthermore, since the fluid loses energy due to electric induction, an electric damping term should be accounted for in the unsteady Bernoulli’s equation. However, we assume that this term is very small and has negligible influence on the fluid dynamics.

In this study, the applied magnetic field is static and varies only along the $x̄$-axis, i.e., $H̄(x̄,ȳ,z̄,t̄)=H̄(x̄)$⁠. The magnetic field is also assumed to have a linear relationship with the magnetization, $M̄(x̄)=χmH̄(x̄)$⁠, that is,²⁰ 

where χ_m is the magnetic susceptibility. Equation (1) can be further non-dimensionalized by introducing the following dimensionless quantities:

where $ω0=πgLδ1$ is the first modal frequency of the fluid column in the absence of surface tension and magnetic field, δ₁ = tanh(πh/L), and $H0=Lω02ρμ0$⁠. This yields the following non-dimensional equations:

where β = σ/(ρgL²), also known as the inverse of the bond number, is the ratio between surface tension and gravitational forces.

Since maximum energy transfer from the base excitation to the fluid occurs near one of the modal frequencies of the fluid column, it is important to characterize variation of the modal frequencies with the design parameters. To this end, we linearize Equations (4c) and (4d) about η(x, t) = 0, and set the external forcing x₀ = 0 to obtain

Differentiating Equation (5a) once with respect to time and substituting Equation (5b) in the resulting equation yields

Equation (6) in conjunction with Equation (4b) admit a solution of the form

Figure 3(a) depicts variation of the first five modal frequencies with the height-to-width ratio, h/L, of the fluid column for a container of width L = 10 cm and ferrofluid which has a surface tension coefficient, σ = 0.03 N/m, and density ρ = 1420 kg/m³. It is evident that the modal frequencies are closely spaced in the considered parameter space. As a result, the bandwidth of frequencies where the harvester cannot respond to the external excitation is very narrow, which could result in a wideband response behavior.

Figure 3(a) also depicts variation of some multiple integers of the lowest three modal frequencies with h/L. Thus, these diagrams permit finding the values of h/L for which the modal frequencies are commensurate. In the vicinity of these points, the sloshing conditions are such that nonlinear interactions among the modes possessing commensurate frequencies can be activated. Consequently, when the nth mode is directly excited via a primary resonance and that mode is in internal resonance with another mode, say the mth, the response will exhibit contributions from both modes even when the mth mode is not directly excited.

It is evident that nonlinear modal interactions are possible near several values of h/L. For instance, when h/L ≈ 0.02; ω₂ = 2ω₁, ω₃ = 3ω₁, and ω₄ = 2ω₂. Furthermore, when h/L ≈ 0.14, ω₄ = 3ω₁, and when h/L ≈ 0.2, ω₅ = 3ω₁. It is worth noting however that, commensuration between two frequencies is necessary but not sufficient; that is, two modes can be commensurate but not interacting. For instance, we have shown previously in Ref. 26 that, internal resonances of the two-to-one type can only occur when ω_m ≈ 2ω_n and m = 2n. In other words, even when the condition ω₃ = 2ω₁ is satisfied when h/L ≈ 0.32, it does not lead to a two-to-one internal resonance between the first and third modes. This is because the activation of a two-to-one internal resonance between any two modes is dependent on the orthogonality of the homogeneous solution of the linear unforced problem to the quadratic coupling terms in the equation of motion. If the mode shapes of the interacting modes are such that the homogeneous solution of the linear unforced problem is orthogonal to the quadratic coupling terms, then the effect of the nonlinear coupling vanishes and energy transfer is not possible. For the problem at hand, the solution of the linear unforced problem is orthogonal to the quadratic coupling terms unless the interacting modes, say the nth and mth, are commensurate and the condition m = 2n is satisfied.

To validate the theoretical model and study the influence of the design parameters, a solution of the nonlinear equations of motion should be obtained. To this end, it can be easily shown that, Laplace’s equation, Equation (4a) subject to the static boundary conditions, Equation (4b) admits a general solution of the following form:

where a₀(t), a_k(t), b_k(t), b_m(t), are unknown functions of time that should be determined by enforcing the kinematic and dynamic boundary conditions at the surface, i.e., Equations (4c) and (4d).

Since Equations (4c) and (4d) are nonlinear due to the advection term ϕ_xη_x in Equation (4d) as well as the kinetic energy and surface tension terms in Equation (4e), an exact solution cannot be easily found. To overcome this issue, an approximate analytical solution of the equations is obtained using the method of multiple scales.²⁷ To this end, the time dependence in the equation is expanded into multiple time scales in the form

where ε = (x₀)^1/2 is a scaling parameter. Based on the definition of the new time scales, the time derivatives can be expressed as

where $Dn=∂∂Tn$⁠.

The dependent variables ϕ, and η, as well as the unknown constant, C are also expanded in the following forms:

The analytical solution is obtained considering the primary resonance behavior of the first mode; in other words, when the frequency of acceleration, Ω, is close to the fundamental frequency; that is Ω = ω₁ + εσ₁, where σ₁ is a frequency detuning parameter. Furthermore, we consider the case when that first mode is in a two-to-one internal resonance with the second mode. Consequently, we express the nearness of ω₂ to twice ω₁, by letting ω₂ = 2ω₁ + εσ₂, where σ₂ is also a small frequency detuning parameter.

Since the dynamic and kinematic boundary conditions are evaluated at the surface η(x, t), which is still unknown, the dependence of ϕ on η is expanded in a Taylor series around η = 0. In other words, we let ϕ(η) ≈ ϕ(0) + εϕ_z(0) η + ε²/2ϕ_zz(0) η² + O(ε³). Note that this assumption is accurate as long as the surface waves are finite but sufficiently small.

Upon substituting Equations (9)-(11) into Equations (4c) and (4d), and collecting terms of like powers of ε, we obtain the following cascade of linear partial differential equations:

O(ε¹):

$D0ϕ1+1πδ1(η1−βη1xx)+χmH2(x)=0,$

(12a)

 

$D0η1−ϕ1z=0,$

(12b)

O(ε²):

$D0ϕ2 +1πδ1(η2−βη2xx)=−12(ϕ1z2+ϕ1x2)−D1ϕ1 +η1D0ϕ1x+xΩ2sin(ΩT0)+C2,$

(13a)

 

$D0η2−ϕ2z=−D1η1+η1ϕ1zz−η1xϕ1x.$

(13b)

In the vicinity of the two-to-one internal resonance, the solution of Equations (12a) and (12b) can be expressed in the form

where C_k(x) = cos(kπx), S_k(x) = sin(kπx), $Chk(z)=coshkπ(z+h/L)coshkπ(h/L)$⁠, A_k(T₁, T₂) is a complex-valued functions that will be obtained by enforcing the solvability conditions at the second stage of the perturbation analysis, and $Āk(T1,T2)$ is its complex conjugate.

Substituting Equation (14) into Equation (13a) yields

where the constants are given by

The secular terms are then eliminated from Equation (15) by enforcing the right-hand side of Equation (15) to be orthogonal to the homogeneous solution, i.e., Equation (14), which yields

where $ρ12=−h̄112h11$⁠, and $ρ11=−h31h12$⁠.

To solve Equation (16), we express the unknown complex functions in the polar form A_n = 1/2a_ne^iβ_n, $Ān=1/2ane−iβn$ and separate the real and imaginary parts to obtain

where γ₁ = σ₂t + β₂ − 2β₁, γ₂ = σ₁t − β₁, $F0=2δ1Ω2π2$⁠, and μ₁, μ₂ are modal damping added to represent viscous damping effects in the fluid. Actual values of modal damping were obtained experimentally using the widely celebrated quality-factor method where the linear frequency response near the first two modal frequencies is fit to experimental findings under very low levels of excitation.

To obtain the steady-state solutions, we set the time derivatives in Equation (17) to zero and solve the resulting algebraic system of equations analytically for the steady-state amplitude, (a₁₀, a₂₀), and phase (γ₁₀, γ₂₀) of the two-mode response. The stability of the resulting solutions is then assessed by finding the eigenvalues associated with the Jacobian of Equation (17) evaluated at the steady-state roots. To first order, the steady-state solution can be written as

where

and $R=4μ22+(σ2−2σ1)2$⁠.

We use the resulting analytical solution, i.e., Equation (18) in conjunction with Equation (19) to study the displacement frequency response curves for two different height-to-width ratios, namely, h/L = 0.2 and h/L = 0.3. Figure 4 depicts variation of the wave amplitude measured at approximately 1 cm from the wall with the excitation frequency for a base acceleration of 0.3 m/s².

For h/L = 0.2, the frequency-response curve exhibits two distinct peaks occurring at approximately 1.7 Hz and 2.25 Hz. The second, higher magnitude peak, occurs near the fundamental frequency of the first mode and is attributed to the primary resonance behavior. On the other hand, the first peak, associated with the lower magnitude, occurs near half the second modal frequency and is due to the two-to-one internal resonance between the first and second modes. Such internal resonance occurs because, as shown in Fig. 3(a), the second mode is nearly twice the first mode when h/L = 0.2. This results in nonlinear energy exchange between the interacting modes.

As shown earlier in Fig. 3(b), in the absence of the second mode contribution, the sloshing wave associated with the first mode is such that the surface has zero velocity at its midpoint. However, as depicted in Fig. 5, the velocity streamlines and surface profile obtained near 1.7 Hz illustrate that, due to the contribution of the second mode, there is a large vertical velocity component at the midpoint of the surface when the internal resonance is activated.

Nevertheless, as shown in Fig. 6, when the excitation frequency is shifted towards 2.57 Hz which is slightly above the first modal frequency, the velocity streamlines, for the most part, illustrate symmetric motions around the midpoint of the surface. This is typical of the first mode response and indicates that the contribution of the second mode diminishes as the excitation frequency deviates from half the second modal frequency.

For h/L = 0.3, the effect of the nonlinear energy pump becomes negligible since, as shown earlier in Fig. 3, the second modal frequency deviates significantly from twice the fundamental frequency. As a result, as shown in Fig. 4, the frequency response only includes one dominant peak near 2.4 Hz.

To obtain the voltage induced in one coil, V_o, Faraday’s law is applied. For a rectangular tank, with coils wound in the vertical direction and magnetic field applied along the $x̄$-axis, the voltage output per unit width of the harvester can be written

where dA is an element on the moving surface $Σ(t̄)$⁠, and B is the vector of magnetic flux density. Assuming small surface waves and that the voltage is induced due to the bulk motion of the fluid and not due to individual dipole rotation, the previous equation can be written as

where 2b is the depth of the container.

The magnetic flux density, B, can be further related to the applied field via $B(x̄)=μ0(1+χm)H̄(x̄)$⁠, which upon substitution into Equation (20) and carrying out the integration yields

To determine the average output voltage generated in the total number of coils, N, we multiply Equation (22) by the number of coils, N, and average the results over the width of the container, L. This yields

where $αc=−2NbhLμo(1+χm)$⁠.

To determine the current induced in the coil, we apply Kirchoff’s Law and obtain, see Fig. 2,

where $ī$ is the induced current, R_l is the load resistance, R_c and L_c are, respectively, the resistance and inductance of the collecting coil.

To investigate the validity of the theoretical model and analytical solution, the experimental setup depicted in Fig. 7 is constructed. A cubic container carrying the ferrofluid with each side measuring 10 cm is placed inside a pick-up coil. The coil is wound around a ferrite core and the whole setup is mounted on an electrodynamic shaker. The external magnetic field is applied using permanent magnets with maximum magnetic field intensity of 92 mT. The spatial distribution of the magnetic field is varied by changing the location of the external permanent magnets with respect to the container walls. The harvested voltage is measured across a resistive load connected in series with the pick-up coil. The physical properties of the ferrofluid and the harvester are listed in Table I.

The height of the surface wave is measured at the corner of the container using a laser vibrometer. Since ferrofluids absorb the laser beam, a thin layer of superhydrophobic coating is sprayed over the ferrofluid surface. To reflect the laser beam, an additional thin layer of a high gloss yellow paint is then sprayed over the hydrophobic coating.

Two cases are analyzed, namely, h/L = 0.2, and h/L = 0.3. In both cases, the steady-state amplitude of the surface wave at approximately 1 cm from the wall is measured. The magnitude of the output voltage across an optimal load of 254 Ω is also recorded under a base acceleration of 0.3 m/s² and a frequency bandwidth around the first modal frequency. In the experiment, the frequency of motion of the shaker head is varied from a value well below resonance to a value well above resonance and vise versa. At each frequency step, the steady-state is allowed to develop. The magnetic field $H̄(x̄)$ is applied using static permanent magnets with opposite polarities placed at a small distance from the container walls. The resulting field across the container is measured using a Gaussmeter and recorded then fitted into a quadratic curve as shown in Fig. 8.

First, the displacement frequency response curves are obtained experimentally and compared to the analytical solution of the equations of motion as shown in Fig. 9. A good qualitative agreement between the model and the experiment is observed. For h/L = 0.2, both the theory and experiment predict the presence of two distinct peaks in the response. However, the theoretical model over-predicts the location of the larger peak associated with the first mode. This could be due to an error in the experimental measurement of the height of the fluid column due to surface tension effects.

Similarly, when h/L = 0.3, both the theory and experiment agree qualitatively in predicting the presence of only one distinct peak in the frequency response. However, the theoretical results again over-predict the response magnitude over most of the frequency range considered.

Next, we compare the voltage response curves in both cases as depicted in Fig. 10. Good qualitative agreement between the theoretical and experimental results is observed. In both cases, the theoretical results under predict the experimental findings near the first modal frequency and slightly over-predict them near half the second modal frequency. When h/L = 0.2 larger voltage output is produced as compared to the case involving h/L = 0.3. Moreover, quite interestingly, even though the displacement response curve did not exhibit a clear peak near half the second modal frequency when h/L = 0.3, both the theoretical and experimental voltage responses exhibit clear response peaks near this frequency. This signifies that the high-frequency small-magnitude motions occurring near twice the excitation frequency contribute significantly to the output voltage.

To further validate the model under a different magnetic field, the experiment is repeated for h/L = 0.2 but using a different applied field, $H̄$⁠. The magnets are inverted such that the two magnets have similar polarity which resulted in the magnetic field profile shown in Fig. 11(a). The model results presented in Fig. 11(b) again show good qualitative agreement with the experimental findings. Both the theory and experiment predict significant amplification in the output voltage resulting from changing the spatial distribution of the magnetic field. The theoretical results, however, overestimate the nonlinearity resulting in more bending of the primary resonance peak.

This article developed a theoretical nonlinear model which governs the electro-magneto-hydrodynamics of a ferrofluid based energy harvester. An approximate analytical solution of the model is obtained using the method of multiple scales for a case involving a two-to-one internal resonance between the first two sloshing modes. The results of the theoretical model are compared to experimental findings for several design parameters. The comparison revealed very good qualitative agreement between the model and experiment, and also indicated some quantitative deviations. Such deviations could have resulted from the different assumptions invoked on the analytical model. First, it was assumed that the angle of contact between surface line and the container is ninety degrees and that this point is free to slip on the surface. Actual conditions may deviate from this assumption. Second, it is assumed that the voltage is generated due to the bulk motion of the fluid. However, voltage can also be generated due to the individual rotation of the magnetic dipoles. Third, only two-mode nonlinear interactions were considered. However, when inspecting Fig. 3, it becomes evident that other nonlinear interactions are possible in the vicinity of the height-to-width ratios considered in the experiments. Finally, the model neglected the backward coupling resulting from the electric damping. This can actually explain why the wave heights obtained analytically over predict the experimental data.

This material is based upon work supported by the National Science Foundation under CMMI Grant No. 1335049. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author and do not necessarily reflect the views of the National Science Foundation.