2.1 Electrical circuit

The electrical circuit consists of a voltage source (LiPo batteries) of internal resistance r_i, maintaining an electric potential U₀ between two rectangular electrodes of size l_x × l_y × l_z separated by a distance H in the y-direction. These two electrodes generate an electric current in the fluid of conductivity σ (the fluid being salt water in our case). In presence of a magnetic field b and assuming that the fluid is flowing at a velocity u_d inside the thruster, the current density can be described using the local Ohm’s law (e.g. [35] or [36] on MHD propulsion) (1) where j is the current density and E the electric field. Assuming uniform fields, the rectangular geometry allows the integration of Eq (1), leading to (2) using |E| = U/H, where U is the electrical potential existing in the fluid. This leads to the global Ohm’s law (3) where r = H/(σl_xl_z) is the fluid electrical resistance in absence of fluid flow, B the volume-averaged magnetic field strength, and k the volume-averaged sinus of the angle between u_d and b. Note that U may differ from the potential U₀ imposed by a voltage source because of the fluid electrolysis chemical reactions, this will be presented in section 2.2.

The evolution of the fluid conductivity σ as a function of salt concentration C can be described using the Kohlrausch law, also called Debye-Hückel-Onsager equation (see [37] for a theoretical derivation of the law). This law reads (4) which gives σ = 0 for C = 0 and C₀ = (a₀/b₀)². Eq (4) also gives a maximum conductivity at the concentration (5) The parameters a₀ and b₀ can be estimated with theoretical expressions [37], but, to be as close as possible from actual measurements, we choose here to fit tabulated values of the literature for a NaCl electrolyte [38]. At 20°C, this gives a₀ ≈ 2071 ⋅ 10⁻⁴ S.m².kg⁻¹ and b₀ ≈ 98.32 ⋅ 10⁻⁴ S.m^7/2.kg^−3/2, leading to the maximum conductivity σ_max ≈ 13.6 S.m⁻¹ for C_max ≈ 197 kg.m⁻³. Note that C₀ ≈ 444 kg.m⁻³ is far larger than the solubility of NaCl in water at 20°C, and cannot thus be reached.

Experiments in salt water show that two supplementary effects have to be taken into account in Eq (3). The fringing effects of the electrodes and the fluid electrolysis. When a voltage difference is applied between two conductive objects (the electrodes here) the generated electric field extends over a distance larger than the electrode itself. This is called the fringing effect, it means that the electric current actually flows on a length larger than l_x, leading to a smaller fluid electrical resistance r* < r. In our configuration, the fringing effect can be estimated analytically using equation (10) of [39], giving , with (6) where ζ = l_y/H, and l_y is the electrode thickness. The actual fluid resistance (without fluid flow) is thus rather .

2.2 Fluid electrolysis

Imposing a voltage difference between two electrodes in an electrolyte (e.g. salt water) leads to fluid electrolysis, i.e. non-spontaneous electrochemical reactions driven by an electric current. In a NaCl solution, we can consider the water reduction to hydroxide and hydrogen gas at the cathode (see [40, 41] for details) (7) and oxidation of chloride to chlorine (8) at the anode. The overall electrolysis of aqueous NaCl results in hydrogen and chlorine gas formation and can be written as (9) Note that the standard potential of Na⁺ reduction is E⁰ = −2.71 V whereas the one for reduction of water is E⁰ = −1.23 V which means that, in aqueous solutions, water reduction will prevail. Chlorine gas Cl₂ will rapidly be dissolved in water (giving ions, see e.g. [42]) meaning that hydrogen gas H₂ would mainly be observed in this case. Depending on the electrode material, reactions implying the electrodes can also occur. For instance, aluminum electrodes can be damaged by oxidation [33] (10) with a standard potential of −1.66 V. Note that stainless electrodes can be used, but copper electrode should be avoided because an oxyde or chloride film is formed on the copper anode limiting the current [33]

Having the lowest standard potential, the water oxidation prevails onto other reactions, and the electrolysis starts as soon as the voltage is larger than 1.23 V. However, any voltage larger than 2.71 V will drive the three reactions. Note that the production of hydroxide OH⁻(aq) at the cathode during the electrolysis will drastically increase the basicity of the solution.

Over-potentials at the electrodes, due to electrolysis driven chemical reactions, can be estimated with the Butler-Volmer equation [43]. Two limiting cases of this equation appear in low and high over-potential regions. The high limit gives the so-called Tafel equation where the over-potential δU is given by (11) where A₀ is the so-called Tafel slope, j is the current density, and j₀ is the so-called exchange current density [43]. Using j ≈ I/(l_xl_z), Eq (11) can be written under its usual form (12) where E₀ = −A₀ ln j₀ and A₀ are two constants related to the reactions involved during the electrolysis. As discussed previously, to start the electrolysis of NaCl electrolytes, the voltage has to be larger than 1.23 V, and one can thus expect E₀ ≈ 1.23 V [33]. The Tafel slope associated with the water oxidation is 0.3 according to [14], and we thus expect A₀ = 0.3 (e.g. [14]) as a typical value. This value also agrees with the detailed over-potential measurements of [44], who note that the Tafel slope A₀ decreases as the conductivity increases. Thus, for a given current density, lower over-potentials are expected with higher conductivity solution.

Finally, the total electrical potential U₀ is the sum of the over-potential δU and the electrical potential U in the fluid. Using Eqs (3) and (12), and taking fringing effect into account, it leads to (13) where and is given by Eq (6). Note that Eq (13) can be solved for the current, leading to (14) with LambertW the Lambert W function and the voltage V = U₀ − E₀ − ku_dBH, where V > 0 for a propulsive thruster.

For very large currents, the rate of electrolysis driven reactions is large, and it may be questioned if this could reduce the local ions concentrations, which would limit the current. This phenomenon is actually quite common in electrolysis, leading to a limiting current density, reached when the electrolysis have consummated all the reactant present in the thin diffusion layer in contact with the electrodes. It can also be noticed that, in presence of a magnetic field as our case, the ions transporting the electric current are deviated by the Lorentz force. The use of a scalar conductivity σ, hidden in the scalar total resistance R in Eq (13), which is considered as independent of the magnetic field, can then be a priori questioned. Finally, a strong flow may influence the electrolysis reactions, which also questions the validity of the over-potential terms in Eq (13). These three important questions are discussed in S1 Appendix, where it is shown that Eq (13) remains valid in the usual ranges of parameters, i.e. for our small scale ship models as well as large scale MHD ships.

2.3 Measurements of the electrical properties

To test the electrical properties of the thruster, a conductometric cell made of two aluminum plates has been built. The electrode size is 2 × 1.7 cm, with a thickness of 1 mm separated one to another by a distance of 7 mm. Using a DC current generator, the current I has been measured for 10 voltages U₀ in the range 0.5 − 10 V, and for various concentrations of salt in the range 5 − 35 kg.m⁻³. For a given concentration, we have checked that our data can successfully be fitted with a simplified version of Eq (13), where the magnetic term has been removed, i.e. with a function of the form U₀ = E₀ + A₀ ln I + RI. These fits give the values of E₀, A₀ and R shown in Fig 3. The data are in excellent agreement with the expected theoretical values. Moreover, Fig 3 shows also that, taking into account fringing effect on the electrodes considerably improve the predictions of the fluid’s electrical properties with salt concentration. One can thus conclude that the theory predicts correctly the electrical behavior of the thruster in the ranges considered here, without any adjustable parameter.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Evolution of E₀, A₀ and R with the NaCl concentration C.

The solid lines correspond to expected values (section 2), i.e. A₀ = 0.3 V (lowest one), E₀ = 1.23 V (intermediate one) and R = r* (the uppermost one). The dashed line shows the resistance r, i.e. the resistance without taking fringing effect into account.

https://doi.org/10.1371/journal.pone.0178599.g003

Since fringing effects can be predicted by Eqs (6) and (13) can be used to obtain the evolution of the conductivity σ from the values of U₀ and I. However the term A₀ ln I makes analytical solution of Eq (13) rather complicated. A simplification can be done using asymptotic behaviors of Eq (13). In the limit of zero current, the flow velocity inside the thruster is zero and thus Eq (13) reduces to U₀ = E₀, on the opposite, for large currents, RI ≫ A₀ ln I and Eq (13) can be approximated by U₀ ≃ RI + ku_dBH. Combining these two asymptotic expressions leads to a simplified version of Eq (13) (15) valid in both limits of small and large currents.

Given the difficulty to measure (U₀, I) on a moving ship, a static ship has been first considered. Indeed, since the term ku_dBH is negligible in Eq (15), these values should not change appreciably for a moving ship. Data are fitted with Eq (15), where ku_dBH has been neglected, to provide R. Then, using Table 1 and Eq (6), the conductivity σ can be estimated for each salt concentration. The results, shown in Fig 4, are in good agreement with the Kohlrausch law Eq (4), obtained from the tabulated values of the literature. All the different electrical measurements are in good agreement with the expected values from the theory, confirming the validity of Eq (13).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Evolution of the conductivity with the NaCl concentration.

Our data (circles) are in excellent agreement with the Kohlrausch law (solid line), given by Eq (4), using (a₀ ≈ 2071 ⋅ 10⁻⁴ S.m².kg⁻¹, b₀ ≈ 98.32 ⋅ 10⁻⁴ S.m^7/2.kg^−3/2). Note that the NaCl solubility in water is 360 kg.m⁻³ at 25°C.

https://doi.org/10.1371/journal.pone.0178599.g004

3.1 Magnetic circuit

The magnetic circuit consists of two neodymium N40 cuboidal magnets, of size L_x × L_y × L_z, separated by a distance W + L_z in the z-direction (such that W is the distance between their surface) and magnetized along z with a residual magnetic field density B_r. The generated magnetic field is then channelized by placing the magnets on a high magnetic permeability (ferromagnetic) U shaped piece of iron which will be called magnetic bridge in the rest of the paper. The magnetic bridge is to magnetic circuit what electrical wires are to electric one; it allows to close the magnetic field. The magnetic field generated by these two cuboidal magnets can be evaluated using the analytical expression of [45], based on [46]. Note that the results fully agree with [47] provided that the magnetization field B_re_z is added within the magnet. The obtained results allow the calculation of the volume-averaged magnetic field B required in section 2. The use of a magnetic bridge is helpful for two main reasons. First, from an experimental point of view, it forces the magnetic field to flow through the bridge, reducing the generation of unwanted forces between the magnetic circuit and the ferromagnetic materials in the surrounding environment. Second, it slightly increases the magnetic flux density by reducing the magnetic reluctance of the circuit.

This small increase of magnetic flux due to the magnetic bridge can be estimated with simple arguments. Given the problem symmetry, we only consider for this estimation a single magnet with a magnetic bridge. The magnetic flux circulating through the magnet (of reluctance ) can flow through different media following two different main paths. It can either flow in the magnetic bridge (of reluctance ) then through the gap between the iron bridge and the magnet (of reluctance ), or flow only in the surrounding medium (of reluctance ), which can be seen as the salt water reluctance in the absence of iron bridge. This ‘choice’ of two possible ways for the magnetic flux leads to the association in parallel of and for the total water reluctance, in series with , which gives (16) Since the magnetic bridge has a high magnetic permeability, is a priori negligible. In the limit of very small gap, , whereas in the limit of infinite gap, we recover . In the first case, the magnet is ‘shunted’, and generates the magnetic flux density B_r, whereas the field of an isolated magnet is recovered in the second case. The magnetic flux is thus increased by the bridge of a factor given by (17) Designating f = B_r/B₀ the ratio between the mean magnetic flux density generated by the shunted magnet (i.e. B_r) and the mean magnetic field generated by the isolated magnet (i.e. B₀), the magnetic flux conservation gives . Thus, noting , we obtain , and the magnetic flux amplification factor reduces to (18) where the factor f depends on the magnet shape. However, magnets are commonly designed with the Evershed criterion, which prescribes f = 2 for magnets with a linear demagnetization curve of the form , where is the magnetizing field. Most modern neodymium magnets present a demagnetization curve of this kind, and their recoil lines are thus very close to their demagnetization curve: these modern magnets operate thus along their demagnetization curve. Using the Evershed criterion, which maximizes the energy product BH in the magnet by finding the rectangle of maximum area fitting below the demagnetization curve , one can obtain that is maximum for B = B_r/2, i.e. for f = 2. Indeed, using the equation y = A₁x + A₀ for our linear demagnetization curve, the product xy = x(A₁x + A₀) is maximum for (x, y) = (A₀/(2A₁), A₀/2), thus for B = B_r/2 if A₀ = B_r (a circular demagnetization curve of equation would maximize xy for , which is close to B_r/2.).

Using the formula of [45], f can be calculated for different magnet dimensions (Fig 5). To channel the flow, the MHD thruster dimension in the flow direction Ox is expected to be large compared to the two others in the normal directions (Oy and Oz, the latter being the magnet magnetization direction). Choosing arbitrarily a thruster 3 times longer (Lx) than wide (Lz) or high (Ly), we only consider magnets with L_y/L_x ≤ 1/3 and L_z/L_x ≤ 1/3, i.e. magnets in the rectangle delimited by the dashed lines in Fig 5. In this rectangle, the Evershed criterion shows that L_z/L_x should be as large as possible and L_y/L_x as small as possible. We thus choose magnets with quite large L_z/L_x, but also with a quite large L_y/L_x because of the magnets availability. The magnets used in this study are represented by the red star in Fig 5. Note that a rough but simple analytical estimate of f can be obtained for our magnets by estimating B₀ with equation (1) in S4 Appendix, i.e. with the axial field of a solenoid of radius a and length L_z. At the surface of the magnet (i.e. for z = L_z/2), this provides . With the apparent radius a ≈ min(L_x/2, L_y/2), we obtain f ≈ 3.1, which is quite close to the actual value f ≈ 3.5.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Contour lines of f as a function of the cuboidal magnet dimensions.

The dashed black lines show L_y/L_x = 1/3 and L_z/L_x = 1/3, and the magnets used in this work are represented by the star (Oz being the magnet magnetization direction and Ox the mean flow velocity in the thruster).

https://doi.org/10.1371/journal.pone.0178599.g005

As w vanishes, the magnetic bridge increases the magnetic flux by a certain factor, given by Eq (18), between 1 (isolated magnet, i.e no magnetic bridge) and 2 (shunted magnet, i.e no air gap between the magnetic bridge and the magnet). It is important to note that the magnetic bridge has usually the same length L_x as the magnet, but a smaller thickness δ, i.e. a smaller cross section δL_x, which increases the magnetic field within the bridge (magnetic flux conservation). However, beyond a certain value B_sat for the magnetic field, the ferromagnetic material constituting the magnetic bridge saturates and the magnetic field leaks then in the surrounding environment. To avoid that, the magnetic bridge thickness has to be large enough, i.e. (19) with B_sat ∼ 1 − 2 T for most ferromagnetic alloys. In practice, using iron and neodymium magnets leads easily to saturation, and the gain is then reduced to a smaller gain because of magnetic leaks. Relaxing the hypothesis of negligible assumed above, the magnetic flux conservation gives . Using Eq (16), we obtain that the gain is actually given by Eq (18) provided that w is replaced by (20)

3.2 Simulations of the magnetic properties

Using the commercial software COMSOL (based on the finite elements method), we solved the magnetic flux conservation for a single magnet, two magnets, and for the total magnetic circuit (two magnets and the magnetic bridge). This allows us to check our different formula related to the magnetic circuit, especially the expression of the field generated by a cuboidal magnet (given by [45]). Our simulations also show that the magnetic bridge gives a gain of 60% in absence of magnetic saturation, in good agreement with the value given by Eq (18).

Then, using a Hirst GM05 Gaussmeter, we measured a magnetic field of ≈300 mT at the distance W/2 of the surface of the magnet’s north pole, in the middle of the thruster. With a single magnet, and in absence of magnetic bridge, COMSOL simulations result in a field of B_z = 124.26 mT, in excellent agreement with the field B_z = 124.04 mT given by the formulas of [45]. One can deduce the experimental magnetic bridge amplification factor , in good agreement with our theoretical gain of 22%.

3.3 Measurements of the magnetic properties

To exploit further the static ship measurements used to determine the electrical properties of the thruster, we also attempted to measure the ship traction force F with dynamometers (bollard pull tests). However, most of our measurements of F were largely disturbed by radiated/reflected waves, residual currents, etc. Reliable results have only been obtained for the largest values of C because F becomes then significant compared to perturbations. These results are shown in Fig 6.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Evolution of the traction force with the current intensity.

The theory predicts F = IBH, leading to a collapse of the data on a single line (BH being constant here). Note that the fit of the data gives a non-zero intercept, probably due to friction in our dynamometers.

https://doi.org/10.1371/journal.pone.0178599.g006

Since the theoretical force generated by the thruster corresponds simply to the Lorentz force, generated by electromagnetic fields, F = IBH, all the measured traction forces are expected to collapse on a straight line of slope BH. A linear fit of the data shows that the intercepts are actually non-zero and negative, indicating an internal friction of 0.04 − 0.05N in our dynamometer. By contrast, the slopes are quite close, in the range 3.2 − 3.6 mN.A⁻¹, in satisfying agreement with the theoretical slope BH ≈ 4.2 mN.A⁻¹ (error of 20%).

4.1 Fluid velocities and thrust

The MHD ship is propelled by a thruster which combine the electrical and magnetic circuits described in sections 2 and 3, respectively. In the thruster, the electric current along Oy and the magnetic field along Oz generate a Lorentz force in the x-direction. This force generates a flow which propels the thruster, and thus the ship, in the opposite direction, at the velocity −u_∞e_x. In the frame moving with the thruster at −u_∞e_x, the mean flow velocity in the thruster is noted u_d. A schematic representation of the thruster and the different fluid velocities and cross section areas used for this study is shown in Fig 7. Note that, contrary to the representation in Fig 7, the present theoretical framework is developed for any kind of velocity profile inside the thruster. Nevertheless, as detailed below and in S2 Appendix, a uniform profile is a correct approximation for this configuration.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Sketch of the thruster and the different velocities taken into account in this study.

Four vortices are displayed to represent the singular head loss at the entrance and exit of the thruster.

https://doi.org/10.1371/journal.pone.0178599.g007

The cross section area in the thruster is constant, equal to S_d = Hl_z. Since the fluid is assumed to be incompressible, volume conservation implies that the mean velocity u_d in this section is also constant. Far from the thruster, the mean fluid velocity is u_∞. Using volume conservation, the section S_∞ of fluid drawn in the thruster is thus given by (21) where λ is the velocity ratio (or the section areas ratio).

At steady state, the momentum conservation for the whole system (water and ship), reduces to the balance between the drag F_D of the ship and the mean thrust F of the thruster. The drag F_D is given by (22) where ρ is the fluid density, S_w the ship cross-section used for the ship drag, and C_d(u_∞) the total drag coefficient of the ship. To estimate C_d(u_∞), it is customary to write it as the sum of the skin friction, form drag and wave making drag coefficients (23) Note that is related to the force created by the friction between the fluid and the surfaces over which it is flowing. This skin friction drag force is usually estimated by considering the drag force generated by a fluid flowing over one side of a flat plate (parallel to the flow), which leads for instance to the Blasius law [48], (24) where S_wet is the total surface area of the plate in contact with the fluid, and Re_L = u_∞L/ν the Reynolds number based on the length L of the plate. Note that formula Eq (24) is typically valid for Re_l < 5 ⋅ 10⁵. Except for the internal wet surface of our propulsive thruster where the skin friction is taken into account through major head loss, our ship can be considered as a sum of flat plates of length L_i, leading to the total skin friction drag coefficient in Eq (22), (25) with , where θ_i is the angle between the plate and Ox, and L_i (resp. δ_i) is the length of each plate along Ox (resp. in the direction perpendicular to Ox). In our experimental setup, we have (in mm^3/2) (26) where is the immersion depth of the thruster top (in mm in this formula).

The coefficient , related to the formation of a wake, depends on the exact shape of the hull and is typically of order 1. Focusing on our experimental ship, one can confirm that . Indeed, for an immersed circular disk at Re > 10³, i.e. for the cylindrical support of the thruster, and for immersed 2D wedges with a half-vertex angle of 27° ([49]), which corresponds to the (plumb) stem angle of our hulls.

The third coefficient , related to gravity waves generation, is usually more difficult to estimate [50] and depends on the Froude number , where h is the water depth. In shallow water, it has been shown [51] that corresponding to our experimental regime Fr < 1.

The mean thrust F is given by [36] (27) such that the final balance equation F = F_D is (28) where is the mass flux through the thruster (mass conservation). Note that the velocity u_∞ is not the inlet velocity, as ambiguously called by [35], because the inlet pressure is not a priori equal to the ambient pressure. As pointed out by [14] and [31], this point is the reason why the momentum balance has been sometimes erroneous in the literature anterior to the 90’s (see also [30, 52]).

In Eq (28), β_d (resp. β_∞) is the momentum coefficient, or momentum correction factor, of the flow in the thruster (resp. of the upstream flow). For a given flow profile u through a section S, this coefficient is defined by , where is the mean flow. This correction factor is used to take into account non-uniform velocity profile when integrating the momentum equation. Typical values are β = 4/3 (resp. β = 6/5) for the Poiseuille flow in a cylindrical (resp. 2D plane) duct, and β = 1 for a uniform velocity profile. However, the magnetic field will modify slightly the velocity profile in the thruster, leading to different values of β. To quantify this effect, one can consider the usual Hartmann flow between two planes. It shows that the magnetic field makes the flow closer to a uniform velocity profile, and thus β is closer to 1 (see S2 Appendix for calculation details). In any case, β remains close to 1, and we thus consider β_d ≈ β_∞ = 1, i.e. uniform flows, when comparing with our experimental results.

Substituting Eq (21) into Eq (28) gives (29) which has to be solved numerically for arbitrary C_d(λu_d). For the analytical calculations to be tractable, a constant C_d can be assumed in Eq (29), leading to (30) which relates the ship velocity u_∞ and the thruster outflow mean velocity u_d (in the frame of the thruster). In the limit , Eq (30) shows that . On the other hand in the limit where , λ tends to λ = 0 (as λ ≈ β_d/β_∞ according to the next order).

In order to close the system, another equation relating u_∞ and u_d is needed. Using the steady Navier-Stokes equation (31) where ν is the kinematic viscosity, u the velocity field, and p the pressure, an average generalized Bernoulli equation can be derived when averaging over the thruster section S_d. This equation balances the flow kinetic energy gain through the thruster with the averaged total work of the Lorentz force and the total head loss χ. It can be expressed as (e.g.[36]) (32) Where is the volume-averaged sinus of the angle between j and b, and α_d (resp. α_∞) the energy coefficient, or kinetic energy correction factor, of the flow in the thruster (resp. of the upstream flow). For a given flow profile u through a section S, this coefficient is defined by , where is the mean flow velocity. As for the momentum correction factor β, α is used to take into account non-uniform velocity profiles when integrating the momentum equation. In the configuration studied here, α remains close to one (see S2 Appendix for details) and will be approximated to unity for the rest of the study.

The total head loss χ in the thruster consists of a linear head loss Λ₁ corresponding to the viscous friction of the fluid on the walls of the thruster, and two singular head losses, Λ₂ and Λ₃, due to the thruster entrance and exit, respectively. We thus have (33) The linear head loss, due to the relative velocity u_d of the fluid with respect to the wall, is given by the Darcy Weisbach equation i.e. (34) where f_D is the Darcy friction factor, and D_h the hydraulic diameter of the thruster. For a pipe, D_h is simply the internal diameter, but for non-circular ducts, D_h is rather given by the estimate D_h = 4A/P, with A the duct cross section, and P the duct perimeter of cross section. For a rectangular duct, an even more accurate estimation is given by the Huebscher formula [53], (35) which differs from equation D_h = 4A/P by 10%.

Various expressions exist in the literature to estimate the Darcy friction factor f_D, depending on the Reynolds number Re = u_dD_h/ν. For Re < 2300, f_D is well estimated by the Hagen-Poiseuille law giving f_D = 64/Re, whereas the Blasius estimate f_D = 0.3164/Re^1/4 is a good approximation for 4000 < Re < 10⁵ in smooth pipes. Many other formulas exist, especially to take the pipe roughness into account. In any case, it is important to notice that the fluid velocity is required to calculate f_D, which introduces a supplementary non-linearity in the system. Note also that the dissipation, and thus f_D, is modified by the magnetic field. However, the Hartmann number H_a is around 1 in our experiment (see S2 Appendix for details), and the dependency of f_D with the magnetic field can thus be neglected at leading order (see [54] for quantitative estimates of this effect, e.g. their Fig 17).

Tabulated excess head coefficients or singular head loss coefficients, for pipe entrances and exits, are used to estimate the singular head loss Λ₂ + Λ₃. These values have been obtained for a fluid at rest entering a pipe, or a pipe outflow in a tank of fluid at rest. To use these values, we have to consider the mean fluid velocity in the inertial frame of reference, where the fluid surrounding the thruster is at rest. In the inertial frame of reference, the mean flow velocity in the thruster is u_d − u_∞. Given that the singular head loss coefficient for an inward projecting/re-entrant (protruding pipe in a tank) is 0.78, and the one for a pipe exit is 1, we thus have (36) where ξ = 1.78 is the sum of the singular head loss coefficients. Note that the presence of other singular head losses would simply modify Eq (36). Here again, ξ is a priori modified by the presence of a magnetic field, a effect which is neglected here because of the moderate values of the Hartmann number H_a reached in our experiment (see S2 Appendix for details).

One can wonder if the head loss associated to the entrance and exit of the thruster could have been reduced using a different geometry in these zones. When a fluid exits a pipe into a much larger body of the same fluid, the velocity is reduced to zero and all of the kinetic energy is dissipated, thus the losses in the system are one velocity head, regardless of the exit geometry. However, the entrance loss coefficient can be made very small by using an appropriate rounded entrance geometry, which gives the lower bound ξ ≥ 1. Note that the ship drag may be significantly increased by a different entrance geometry, reducing the overall interest of such a modification. One may wonder if the electrolysis may modify the flow or the ship drag via the bubbles generation. According to [42], the mean bubble diameter is 1 − 100 μm and the volume gas fraction is 10⁻⁴ − 10⁻³ for seawater electrolysis in conditions close to our experimental setup. According to [55], we do not expect any influence of the bubbles on the drag, thus on the flow, which is confirmed by the very good agreement between our theory and our experimental results (see section 4.4).

4.2 Summary: Thruster governing equations

As shown in the previous sections, the thruster dynamics can be described as follow: the imposed voltage generates an electric current, given by Eq (13), which generates a flow governed by Eq (32), inducing an opposite electric current by a feedback term in Eq (13). Here, the magnetic field generated by this flow induced electric current is thus neglected. The validity of this approximation can be estimated using the magnetic Reynolds number R_m = udH/ν_m, with ν_m = (σμ)⁻¹ the magnetic diffusivity and μ the fluid magnetic permeability. When R_m ≪ 1, the induced field is negligible compared to the imposed one, and the unknown current I and velocity u_d are then given by Eqs (13), (29) and (32). Thus, the three unknowns λ, u_d and I are governed by (37) (38) (39) where (40) (41) allowing the calculation of u_∞ and S_∞ using Eq (21). For a static thruster, u_∞ = 0, so λ = 0, and the problem is then only governed by the two Eqs (38) and (39). Note also that the inlet velocity (and pressure) is not involved in these equations, but can be calculated a posteriori using volume conservation and Bernoulli equation.

Following [32], it is of interest to estimate how the thruster electrical efficiency varies. Where is the mechanical power imparted on the fluid and the electrical power given by the LiPo Battery. In the literature, the load factor (42) is often introduced [1], which is the ratio between the effective voltage imposed to the fluid and the voltage induced by the flow. Using the load factor K, η can conveniently be written as (43) where K₀ = E₀/(ku_dBH) is usually neglected in the literature (E₀ ≪ U₀). Since K > 1 for a thruster, a good efficiency is reached for K ≥ 1 [36]. Considering the simpler Eq (15), K reduces to K = 1 + RI/(ku_dBH), which gives (44) i.e. (45) with the volume of water in the thruster. Under the usual assumptions of the literature (K₀ = 0, k = 1, , r_i = 0), Eq (45) reduces to (46) which shows that, for a given velocity u_d, the efficiency approaches 1 when is increased [1]. Maximizing η by deploying the highest possible magnetic field in the largest available volume is actually common to the four families of magnetohydrodynamic thruster [4, 5, 56–59]. However, Eqs (43) and (44) also give as [60] (47) showing that approaches 0 when the second term of Eq (45) is minimized. Thus, maximizing η gives a vanishing thrust. As pointed out in [60], rather than η, one should thus optimize . Neglecting K₀, the value of η which maximizes is η = 1/(2k), which corresponds to a load factor K = 2. With k = 1, one thus should expect that 50% of the electrical power is consumed in Joule heating.

Considering typical speed and size of commonly used ships, it is interesting to estimate how an MHD thruster can compete with usual propulsion method. Since its efficiency increases with B (Eq (46)), one can thus determine the typical magnetic field required to obtain an acceptable efficiency for these ships. Using the typical ship length L and velocity u_d, Eq (46) gives (48) with , , and where τ = L/u_d is the typical time corresponding to the time required for the ship to move for a distance equal to its length. Note that the field required by Eq (48) depends on the typical time scale only, and not on the ship length. Considering a ship of length L = 10 m traveling in seawater where ρ = 10³ kg.m⁻³ and σ = 5 S.m⁻¹, at a velocity u_d = 10 m.s⁻¹, leading to a typical time τ = 1 s, Eq (48) gives a magnetic field B = 10 T to maximize . The same magnetic field will be required for a small scale ship of size 10 cm traveling at u_d = 10 cm.s⁻¹ since the typical time τ = 1s remains the same. Practical magnetic field intensity which can be achieved are typically smaller than B = 10 T, therefore the current MHD thrusters cannot compete with usual propulsion method. Conversely, maximum magnet field strength being typically 0.1 − 1 T, using the same typical values in Eq (46) shows that one can expect an efficiency η ≈ 0.01 − 1%.

4.3 Analytical solutions for the thruster

Eqs (37)–(39) can be solved numerically but analytical solutions are much harder to obtain. In order to make analytical progress, the simplified electrical Eq (15) is considered, the Darcy friction factor is supposed to be f_D = 64/Re, which corresponds to a laminar flow in the thruster (Re < 2300), and C_d is assumed to be constant. This latter hypothesis allows to solve Eq (37), leading to the solution Eq (30) for λ. Note also that this latter hypothesis is not needed to obtain analytical solution for the static thruster since Eq (37) is irrelevant in this particular case. As shown in S3 Appendix, analytical solutions of the system of Eqs (15)–(39) can then be obtained under these assumptions.

Based on order of magnitude arguments, results on the different solutions can be obtained. For instance, K being dimensionless implies that a typical field B_typ = (U₀ − E₀)/(ku_dH) exists for the thruster. Assuming a balance between U₀ − E₀ and the other terms of Eq (15) gives I ∼ (U₀ − E₀)/R and then , allowing the evaluation of the typical flow velocity in the thruster u_d. Since K > 1 gives a thruster and K < 1 an electricity generator, it is thus expected that the solutions change for K ∼ 1, i.e. when B ∼ B_typ, with .

For instance, when B ≪ B_typ, the solutions can be reduced to (for f_D = 0) (49) (50) showing that u_d increases with B. This limit corresponds to the limit where the voltage ku_dBH induced by the flow is negligible compared to U₀ − E₀. The analytical solutions also show that u_d decreases with B for large B ≫ B_typ, which shows the existence of an optimal field B_opt for u_d (which is thus bounded when B is varied). This optimum is actually obtained when K = 2, with the optimal field B_opt = 2^−1/3 B_typ ≈ B_typ, which could be expected. This optimum is the same as the one obtained in section 4.2 for , and the efficiency η = 1/(2k) is recovered for E₀ = 0 (see equation (17) in S3 Appendix).

As shown in S3 Appendix, the analytical solutions also predict that the efficiency is maximized for a certain voltage U₀. Looking for B and U₀ which simultaneously maximizes u_d and η, respectively, gives U₀ = 3E₀ and , with an efficiency of η = 1/(3k).

4.4 Measurements of the ship velocity

Time evolution of the ship velocity u_∞ is derived from video recording its displacement for the three batteries and six different salt concentrations, varying from the average seawater salinity (∼35 g.L⁻¹) to the Dead Sea salinity (∼300 g.L⁻¹). Fig 8 shows the displacement of the boat as a function of time for three different experiments. Increasing the salt concentration and battery voltage increases the ship velocity, but, for all configurations, a terminal constant velocity is obtained before reaching the end of the tank.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. Time evolution of the ship displacement for three different experimental configurations, varying salt concentration and battery power (symbols).

Black-dashed lines represent linear fits over the last 6 points of each experiment. The terminal velocity of the boat is defined as the slope of the black-dashed line.

https://doi.org/10.1371/journal.pone.0178599.g008

In Table 2, the measured current and terminal ship velocity are given for all the eighteen experimental runs. Note that the measured current is systematically small compared to the maximum current I_m sustainable by the battery (given in Table 1), insuring that the current is not limited by the battery but by the electrical circuit resistance. Table 2 also shows that, for a given battery, u_∞ increases with the concentration of NaCl, but seems then to decrease for concentration C > 210 kg.m⁻³ (the presence of this maximum is expected, see section 2.3). Note also that the involved battery power is quite large, up to 1 kW, for a maximum terminal velocity of 30 cm.s⁻¹. One can wonder if an internal inductive thruster would have given better performances. Considering an inductive internal thruster of typical size 10 cm, powered by a 200A superconducting magnet rotating at 50 Hz and immersed in a similar fluid (ρ = 1100 kg.m⁻³, σ = 10 S.m⁻¹), the typical velocity (resp. magnetic field) is expected to be ≈10 cm.s⁻¹ (resp. ∼0.1 − 1 T), i.e. comparable to our slowest measurements [61].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Experimental measurements.

https://doi.org/10.1371/journal.pone.0178599.t002

Focusing on the terminal velocity, it is useful to collapse the results on a single curve making the comparison between the different configurations easier. To do so, it should be first noticed that, modifying the battery obviously changes the electric current I, but also the total mass m = m_ship + m_LiPo, thus the immersion depth of the thruster. Hence, the ship velocity u_∞ has to be expressed as a function of I, m and C. The regular head loss being negligible in our experimental setup, Eq (39) gives , leading to (51) For all different ship cross sections S_w considered here, 2β_dS_wC_d/S_d ≈ 30 for our ship configuration, which is large compared to . Then, assuming a constant C_d, Eq (30) gives λ ≈ [2β_dS_d/(S_wC_d)]^1/2, which gives a typical value of λ ≈ 0.35. This allows to simplify Eq (51) into (52) considering the limit λ ≪ 1, and , β_d = α_d = 1, ξ ≈ 2, C_d ≈ 1 (see Table 1). The thruster solid cross section has been neglected in Eq (52) compared to the one of the ship, leading to S_w ≈ γm/ρ, where the constant only depends on the hull geometry (see also Eq (55)).

Eq (52) shows that, at first order, the ship velocity does not depend on the fluid density ρ or the thruster cross section S_d. It also shows that, at first order, the variable I/m should allow a collapse of the data.

Experimental results for each concentration, for the three batteries, are presented in Fig 9 using the variable I/m. The terminal velocity has been obtained from the six last positions recorded by the camera (Fig 8). The approximated Eq (52) is plotted (solid line), as well as the exact theory, plotted for the three batteries (dashed lines), using Eqs (37)–(39), the values of Table 1, and (53) with m = m_ship + m_LiPo and being the immersion depth of the top of the thruster. As expected, the theoretical velocities collapse quite well for the three batteries. One can also notice a certain inflection point on the theoretical curves, corresponding to a concentration beyond which the velocity decreases when C is increased. Without any adjustable parameters, the theory predictions of the maximum velocity reached by the ship are in good agreement with the experiments.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 9. Evolution of the ship velocity as a function of I/m. Experimental results for the three batteries (for each concentration) are represented by symbols.

The three dotted lines, which collapse rather well, correspond to Eqs (37)–(39) for the three batteries. The solid line is given by Eq (52).

https://doi.org/10.1371/journal.pone.0178599.g009

5.1 Expected orders of magnitude

To design the MHD ship, it is useful to study the typical orders of magnitude we can expect for a small scale experimental setup. In order to predict these orders of magnitude a typical length of 10 cm is chosen for the ship and Lithium-Polymer (LiPo) batteries as power source. Neodymium magnets are considered to generate the magnetic field. These magnets are the strongest type of permanent magnets commercially available and a good measure of their strength is given by their grades, defined as their maximum energy product. These grades usually range between N35 and N52, and can be related to B_r using the empirical equation (54) where B_r is given in T, and where X is the grade number (e.g. 40 for a magnet of grade N40). A priori, the highest the grade, the more powerful the ship is, however, a high grade magnet is more inclined to physically break. Moreover, across the whole range of grades, B_r only varies between 1.17 T and 1.48 T. A good compromise between magnet strength and magnet solidity is the grade N40, one of the most common grade and the one used in this study. According to Eq (54), the residual magnetic field density obtained with these magnets is B_r ≈ 1.26 T. Then, using the Evershed criterion (see section 3) a typical field of half this value can be expected (f = 2), thus B ≈ 0.63 T. Using seawater as working fluid, typical values for density, ρ ≈ 10³ kg.m⁻³, minimum voltage needed for electrolysis E₀ ≈ 1 V, conductivity σ ≈ 10 S.m⁻¹, drag coefficient C_d ≈ 1 can be chosen to obtain the orders of magnitude.

These typical values lead to a fluid electrical resistance r = H/(σl_xl_z) = 1 Ω, which is large compared to the typical internal resistance r_i ≈ 0.01 Ω of a LiPo battery. For simplification, fringing effects are neglected so R ≈ r = 1 Ω, and, assuming uniform flows in the thruster, λ ≈ 0.7 which gives . Note that the thruster suction area is S_∞/S_d = 1.4 times larger than S_d, which constrains the experimental tank size. Using these parameters, Eqs (49) and (50) give I = 9A and u_d ≈ 43 cm.s⁻¹, i.e. u_∞ ≈ 30 cm.s⁻¹. This gives a load factor of K ≈ 330 and a ship efficiency of η = u_dBH/U₀ ≈ 0.3%. These expected orders of magnitude are actually confirmed by our experimental results (see section 1). The Reynolds number obtained using this values is Re ≈ 10⁴, which is in the range of validity for the Blasius estimate considered here (see section 4). This also shows that the induced electric field (ku_dBH ∼ 0.03 V with k ∼ 1) is negligible (RI ≈ 9 V), and thus that the simplifying hypothesis k = 0 is fully relevant here. We can therefor safely use Eqs (49) and (50).

In Fig 10, axes have been chosen to gather a lot of information for an arbitrary thruster of size H, operating in a flow u_d. First, if follows from the definition of the load factor K (Eq (42)) that the thick solid tilted line K = 1 separates the two possible behaviors of the thruster, i.e. a flow or a voltage generator. Actually, the solutions for a given K are represented by dotted lines parallel to the solid line K = 1 (here, K = 2 and K = 100 are shown). The thick dashed lines represent how the solutions evolve when B is varied, whereas the star corresponds to our typical MHD small scale ship model. It shows that increasing B increases the velocity, as expected, until the maximum velocity reached for K = 2 (see section 4.2 or S3 Appendix). Beyond K = 2, increasing B decreases the velocity, even if the thruster efficiency continues to increase towards one (the solid curve tends along the line K = 1). Fig 10 also shows that, in any cases, the induced magnetic field can be neglected for usual values of (U₀ − E₀)/B. Note finally that our typical boat remains below the line K = 100, even with very optimistic values of B, which confirms that the induced voltage can be neglected. However, this induced voltage cannot be neglected if we consider a ship of typical size 10 m using a magnetic field of B = 4 T (see solid curve of Fig 10). Such a ship can actually represent the Yamato 1, a ship built in the early 1990s by Mitsubishi Heavy Industries [17]. This prototype was able to reach speeds of 12 km.h⁻¹ with a thruster cross section size of ∼1 m [18, 62]. These values agree with the maximum u_dH ≈ 4 m².s⁻¹ shown in Fig 10, which gives u_d ≈ 14 km.h⁻¹ for H = 1 m.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 10. Thruster working regimes (seawater, k = 1).

Below (resp. above) the thick solid tilted line K = 1, the thruster behaves as a flow (resp. a voltage) generator. The dotted lines parallel to the thick line K = 1 show solutions for constant values of K. Our typical thruster (the star) evolves along the thick dashed curve when B is changed in the range B = 10⁻² − 1 T. The other thick dashed curve is similar, for a thruster size of 10 m and a magnetic field from B = 10⁻² T to B = 4 T.

https://doi.org/10.1371/journal.pone.0178599.g010

5.2 Magnets distance maximizing the velocity

Considering the experimental ship described in Table 1, one can wonder whether, the distance between the magnets, W has been well chosen. Varying W will naturally change B, but also the ship mass and cross section S_w. Noting m_LiPo the battery mass, S_w can be calculated using the thruster cross section S_th and the immersion depth of the ship floats (S_w = S_th + S_hull, with the hull section area S_hull), (55) where m = m_ship + m_LiPo is the ship total mass, S_thL_x is the volume of water displaced by the thruster only, and δm = 2.2(W − 0.024) is the mass difference (in kg) for the ship when W (in meters) varies.

Using l_z = W − ϵ, with ϵ = 6 mm the total thickness of the insulating material, Eqs (38) and (39) can be solved using values in Table 1 and Eq (18) to take into account the W dependency on the magnetic bridge effect, the mass of the ship, and S_th. Using S_th = S₀ + ϵ₂l_z, where S₀ = 960 mm² and ϵ₂ = 8 mm. The results, plotted in Fig 11, show that an optimal W exists.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 11. Ship velocity as a function of W.

For the intermediate battery 4s (the dashed line being the value for our experiments), with σ ≈ 8.7 S.m⁻¹, B_r ≈ 1.26 T (skin friction modifications neglected).

https://doi.org/10.1371/journal.pone.0178599.g011

To optimize the thruster on various parameters simultaneously, a simple estimate of this optimal W is required. However, calculating B as a function of W requires to average in space the exact magnetic field of a cuboidal magnet (given by [45]) which leads to complex calculations. A simpler expression of the magnetic field generated by the magnets can be obtained assuming that the magnetic field generated by a single magnet is given by the axial magnetic field of a cylindrical magnet of radius a and length L_z. Then, the mean field B between two cylindrical magnets can be estimated by (see S4 Appendix for details) (56) Calculations with the exact magnetic field of a cuboidal magnet show that Eq (56) correctly captures the evolution of B with W and hence is quite useful for rapid estimations. For instance, in [25], B is rather fitted by a quadratic polynomial in W, leading the authors to conclude erroneously that an optimum for BW exists. Actually, the exact calculations of B with the formulas of [45] confirm that BW evolves as B_rL_zW/(a + W), i.e. does not present any optimum.

Considering Eq (50), with R ≈ H/(σLW), and (assuming S_d ∝ W, and neglecting regular head losses), gives a constant u_d for small W, but varies as W^−1/2 for large W. Equating the two asymptotic expressions shows that this change of behavior appears around W ≈ a for r_i ≈ 0, which is the optimal W. Since λ ∼ W^1/2 for small W and λ ∼ 1 for large W, u_∞ ∼ W^1/2 for small W, and u_∞ varies as W^−1/2 for large W, exhibiting an optimum around W_opt ≈ a.

5.3 Interelectrode distance maximizing the velocity

As for the distance separating the magnets, considering the values given in Table 1, one can also wonder if the distance between electrodes H has been well chosen. Varying H will naturally change the current I, but also the thruster cross section S_d. The evolution of u_∞ as a function of the separation distance between the electrodes H can be calculated solving Eqs (38) and (39) with the values in Table 1. Using S_w = 3268 mm², corresponding to the case with 4s battery, the evolution of the ship velocity as a function of the distance between the electrodes can be calculated. Fig 12 shows this evolution and an optimum distance of H ≈ 2 cm is obtained, maximizing u_∞.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 12. Ship velocity as a function of H.

For the intermediate battery 4s (the dashed line being the design chosen for our experiments), with σ ≈ 8.7 S.m⁻¹, B_r ≈ 1.26 T (skin friction modifications neglected).

https://doi.org/10.1371/journal.pone.0178599.g012

As for the magnetic field, to optimize the thruster on various parameters simultaneously, an estimate of this optimal H is required. Considering the Eq (50), with R ≈ H/(σLW), and , assuming S_d ∝ H, and neglecting regular head losses, the flow velocity inside the thruster evolves as u_d ∼ H^−1/2 for small H. This unphysical behavior comes from the fact that the regular head losses have been neglected here but are actually not negligible at small H. Using the Huebscher law Eq (35) and the Hagen-Poiseuille law f_D = 64/Re, valid for asymptotically small H, it can be shown that u_d rather varies as H^1/4. If H is larger than a critical value (around the magnet size L_y), the electric current flows in a zone where the magnetic field has an opposite direction, leading to an opposite Lorentz force direction. In this case, the mean field B naturally decreases rapidly. Since λ ∼ W^1/2 for small W and λ ∼ 1 for large H, u_∞ ∼ W^3/4 for small H, and decreases rapidly beyond typically L_y. To be sure to avoid the rapid decrease zone, we have chosen a value of H slightly smaller than L_y, with H = 1.4 cm < L_y = 2.4 cm, close to the optimum H ≈ 2 cm represented in Fig 12.

5.4 Salt concentration maximizing velocities

Since C modifies σ but also ρ, one can wonder if the concentration C_max, which maximizes the fluid conductivity σ (see Eq (5)), is different from the optimal concentration maximizing the velocities u_d or u_∞. Noting ρ₀ the fluid density for C = 0, the density reads as ρ = ρ₀ + C, and thus , where Γ is independent of C. Considering Eq (50) in the particular case of vanishing internal resistance (r_i = 0), u_d is maximum for the concentration (57) which gives Eq (5) in the limit of large ρ₀. Using the values given in section 2 for a₀, b₀, and ρ₀ = 10³ kg.m⁻³, Eq (57) gives C_opt ≈ 176 kg.m⁻³, which is different from C_max ≈ 197 kg.m⁻³. Solving numerically the full Eqs (37)–(39) shows that u_d is indeed maximum around C ≈ C_opt.

The fluid velocity u_d inside the thruster does not dependent on the thruster immersion depth but the ship velocity u_∞ is directly affected by the fluid drag force. The drag force being proportional to the ship cross-section S_w, which varies with ρ, and thus C, the optimum concentration for u_d might not be the one maximizing the ship velocity u_∞. In the parameters range considered in this work, this actually balances the C dependency of , as shown by the estimate Eq (52) of u_∞ where the C dependency only comes from I, i.e. from σ. One can thus finally expect that u_∞ is maximum around C_max ≈ 197 kg.m⁻³. Solving numerically Eqs (37)–(39) gives a maximum u_∞ for C ≈ 191 kg.m⁻³, in good agreement with our measurements (see Fig 9). Note that, if the ship considered in this work was a submarine, S_w would be constant, and the concentration maximizing u_∞ would rather be C_opt ≈ 176 kg.m⁻³.