Obsidian/Zotero/001_artiklid/shanmugadasonsagercoefficientsliquidmetal2025.md

---

category: literaturenote citekey: shanmugadasonsagercoefficientsliquidmetal2025 title: Onsager Coefficients for Liquid Metal Flow in a Conduit under a Magnetic Field authors: "Shanmugadas, Sindu E.; Bau, Haim H." year: 2025 date: 2025-00-00 2025 doi: 10.1063/5.0276573 url: "https://arxiv.org/abs/2506.10132v1" zotero_key: KEME3G5Q zotero_storage: 2EKW7QU2 collections: transpordi_planeerimine folder: 001_artiklid firstAuthor: "Shanmugadas, Sindu E."

status: converted

Copyright ((2025) Shanmugadas and Bau. This article is distributed under a Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International (CC BY-NC-ND) License. The following article has been accepted by Physics of Fluids. After it is published, it will be found at (DOI: 10.1063/5.0276573).

Onsager Coefficients for Liquid Metal Flow in a Conduit under a Magnetic Field

Sindu E. Shanmugadas and Haim H. Bau*

Dept. Mechanical Engineering and Applied Mechanics, University of Pennsylvania, Philadelphia, PA 1904-6315. USA

*Corresponding author: bau@seas.upenn.edu

Abstract:

We analyze the flow of room and near room-temperature liquid metals in shallow, long rectangular conduits with two insulating and two perfectly conducting walls under a uniform magnetic field perpendicular to the flow direction and the insulating surfaces, focusing on moderate Hartmann numbers. A pressure gradient and Lorentz body forces may drive or oppose the flow. We derive explicit expressions for the Onsager coefficients that relate the flow rate and electric current on the one hand to the potential difference across electrodes and the pressure gradient on the other hand. We further demonstrate that these coefficients satisfy Onsager-Casimir reciprocity. These simplified expressions provide a convenient framework for analyzing, optimizing, and controlling magnetohydrodynamic (MHD) machines operating with liquid metals in applications such as power conversion, energy harvesting, pumping, actuation, valving, breaking, and sensing without moving components.

Nomenclature

(Variables without specified units are dimensionless)

Latin Symbols

• a Conduit's half height (m).
• b Conduit's half width (m).
• c Numerical coefficient.
• B0 Uniform magnetic field in the y-direction (T).
• h Normalized induced current stream function.
• G Pressure gradient (= − ) (Pa m-1 ).
• Hz Current stream function (A m-1 ).
• J Electric current density (A m-2 ).
• KEE Onsager coefficient electrical conductance (Ω-1 m-1 ).
• KEH Cross Onsager coefficient, the electrical current generated by a 1 Pa·m⁻¹ pressure gradient (2 −1).
• KHE Cross Onsager coefficient, the flow rate generated by a 1V potential difference (2 −1).
• KHH Onsager coefficient, hydrodynamic conductance (6−1−1).
• L Half conduit length (m).
• p Pressure (Pa)
• P Body force scale ( 0 2 ⁺ ) (N m-3 )
• Q Flow rate (m3 s -1 )
• r Resistance ratio
• Rex External resistance (Ω m)
• Rin Internal resistance (Ω m)
• u Normalized axial velocity (= 2 )

• vz Axial velocity (m s-1 )

• ̅ Axial average velocity =(4ab)-1 Q (m s-1 )

• x, y, z Cartesian coordinates (m)

Subscripts

• c Core (outer) region
• i Boundary layer (inner) region
• m maximum

Greek Symbols

• Numerical coefficients.
• Conversion efficiency
• Length scale for the lateral boundary layer
• Potential (V)
• ̃ Normalized potential (equation 7)
• Dynamic viscosity (Pa s)
• Kinematic viscosity (m2 s -1 )
• Electric conductivity (S m-1 )
• Normalized x coordinate (=a-1 x)
• Normalized y coordinate (=b-1 y)
• Stretched boundary layer coordinate in the x-direction
• Potential difference across electrodes (V)

Dimensionless Groups

• A Aspect ratio (a/b)
• Ha Hartmann Number *(=*0 σ μ )
• Re Reynolds Number ̄
• S Ratio of pressure gradient and Lorentz body force 0∆ 2

Introduction:

In recent years, there has been a growing interest in room-temperature and near-room-temperature liquid metals, particularly gallium-based eutectic alloys, across various fields, including flexible, stretchable, reconfigurable, and self-healing electronics; soft robotics; energy harvesting and conversion; on-demand thermal control; wearable sensors; and microfluidics [1]. Gallium-based alloys are attractive due to their high electrical and thermal conductivities, low toxicity, and biocompatibility.

Earlier, microfluidics researchers leveraged magnetic fields and electric currents to induce Lorentz body forces in aqueous solutions, driving liquid flow [2, 3, 4], stirring liquids [5, 6, 7], and controlling fluidic networks [8] without any moving members or mechanical components. However, aqueous solutions have relatively low electric conductivity and are susceptible to electrochemical reactions at electrode surfaces, significantly limiting the strength of the Lorentz body forces that can be produced. In contrast, liquid metals, such as gallium-based alloys, offer a significant advantage due to their electrical conductivity being orders of magnitude higher (~106 S m-1 ) compared to aqueous solutions (~10 S m-1 ), enabling the generation of substantial Lorentz body forces even under moderate magnetic fields, such as available from small rare-earth, permanent magnets.

The use of eutectic gallium alloys comes with several challenges. A primary concern is that gallium and its alloys readily form a surface oxide layer when exposed to oxygen. This oxide layer reduces surface energy, hindering the flow and wetting behavior of the liquid metal and adversely affecting contact electrical resistance. Strategies to mitigate the formation and effects of this oxide layer include operating in controlled environments, such as in inert (oxygen-free) atmospheres; using acidic conditions to dissolve the oxide layer [9]; using reducing electrochemistry to reduce gallium oxide; and applying coatings to the liquid metal surface to prevent gallium's exposure to oxygen [10, 11, 12]. Additionally, gallium is corrosive to certain metals, such as aluminum and copper. Gallium infiltrates the grain boundaries of these metals, leading to embrittlement and structural degradation over time. Protective coatings, such as barrier layers, can be applied to susceptible metals, preventing direct contact with gallium [10, 11, 12].

Despite these challenges, applications of room-temperature and near-room-temperature liquid metals MHD are expanding. Although liquid metals have 10-fold lower electrical conductivity than solid copper and cannot compete with solid copper in conventional electromagnetic machines, they are compatible with flexible, reconfigurable, compact devices. Applications include, among other things, harvesting energy from ocean waves [13], low-temperature heat sources [14], motion [15, 16], and body locomotion [17]; ondemand thermal control [18, 19]; actuation, sensing, and assisting in fabricating soft robotics and wearable devices.

Many minute systems that utilize liquid metal MHD comprise shallow, long conduits with rectangular cross-sections (− ≤ ≤ , − ≤ ≤ ) with aspect ratio: = ^&lt;1; electrodes at ⁼ ±; and uniform magnetic field (0) in the y-direction (Fig. 1). Permanent, rare-earth magnets often provide the magnetic field. Flows in these devices are usually laminar (Re << 103 ). Since the magnetic Reynolds number and the current flux are small, distortions to the externally applied magnetic field 0 are insignificant [20]. The Hartmann number = 0 σ μ , gauging the relative importance of induced Lorentz and viscous forces, is O(10). In the above, and are, respectively, the liquid metal electric conductivity and dynamic viscosity.

Except for the asymptotic limits of small [21] and large [22] Hartmann numbers, the MHD equations must be solved numerically. The analysis, design, optimization, and control of magnetohydrodynamic systems, such as generators and energy harvesters, pumps, actuators, breaks, and sensors, would benefit from simple relationships between control inputs, such as potential across electrodes and flow assisting/retarding pressure gradient on the one hand, and outputs, such as flow rate and current, on the other hand. Our paper aims to provide such relationships.

Our paper is organized as follows. Section 2 presents the mathematical model for fully developed MHD flow in a long conduit with a rectangular cross-section and solves the MHD equations by finite elements. Section 3 carries asymptotic analysis in the limit of aspect ratio: → 0. The result of this analysis enables us to obtain relationships between flow rate and current on the one hand and potential difference across the electrodes and pressure gradient on the other hand in the form of Onsager relations that satisfy Onsager-Casimir reciprocity. Section 4 takes advantage of these simple relations to analyze the performance of various MHD machines.

Fig. 1: Our long flow conduit − < < with the uniform, rectangular cross-section − ≤ ≤ , − ≤ ≤ . The uniform magnetic field 0 is in the positive -direction. The surfaces = ± are perfect conductors (electrodes). The surfaces = ± are insulating.

2. Mathematical Model

Consider a long duct (− < < ) with the rectangular cross-section − ≤ ≤ , − ≤ ≤ ; (Fig. 1). The coordinate aligns with the conduit's axis. The conduit's surfaces = ± are perfectly conducting electrodes, each at a uniform potential and in perfect contact with the liquid metal. The surfaces = ± are perfectly insulating (blocking electrical current). The duct is filled with a conductive liquid with uniform electric conductivity and viscosity . The conduit is subjected to a uniform magnetic field B0 in the -direction, unaffected by induced magnetic fields. Since the magnetic Reynolds number and the current in the duct are small, so is the magnetic induction. Thus, the externally applied magnetic field remains unaffected. The steady flow in the conduit is assumed to be nearly isothermal due

to the liquid metal's high thermal conductivity and the efficient heat exchange with the duct walls, facilitated by the duct's shallow geometry. We modeled the MHD flow using the time-independent Navier-Stokes equation with non-slip boundary conditions at all solid surfaces. ≫ (, ). At a distance O(a Re) from the inlet, the flow is fully-developed, such that ^{= 0}. When the flow is fullydeveloped, the momentum equation in the -direction reduces to:

$$0 = -\frac{dp}{dz} + \mu \left( \frac{\partial^2 v_z}{\partial x^2} + \frac{\partial^2 v_z}{\partial y^2} \right) + J_x B_0, \tag{1}$$

where p is the pressure, the x-direction component of the current density, and (, ) the axial velocity. The pressure gradient = − is constant.

The electric current density accounts for both Ohmic and induced currents.

$$J_x = \sigma \left( -\frac{\partial \phi}{\partial x} - v_z B_0 \right), \quad J_y = \sigma \left( -\frac{\partial \phi}{\partial y} \right),$$ (2)

where (, ) is the electric potential. Charge conservation requires the current density J to be divergence-free ( ∙ = 0). We use bold letters to describe vectors and regular font for scalar quantities. Following Hunt-Stewartson (hereafter, HS) [22], we introduce the current stream function (, )

$$J_x = \frac{\partial H_z}{\partial y}, \qquad J_y = -\frac{\partial H_z}{\partial x},$$ (3)

which satisfies charge conservation. With the introduction of , equations (1) and (2) reduce to two coupled, linear pdes with vz and Hz as the dependent variables.

$$0 = -\frac{dp}{dz} + \mu \left( \frac{\partial^2 v_z}{\partial x^2} + \frac{\partial^2 v_z}{\partial y^2} \right) + B_0 \frac{\partial H_z}{\partial y}. \tag{4}$$

$$0 = \frac{1}{\sigma} \left( \frac{\partial^2 H_z}{\partial x^2} + \frac{\partial^2 H_z}{\partial y^2} \right) + B_0 \frac{\partial v_z}{\partial y}. \tag{5}$$

At the insulating surfaces, (, ) = 1, (, −) = 2, and ∆ = 1 − 2 is the net current flow between the opposing electrodes per conduit's unit length. ∆ is known a priori only when operating under current control.

We use b as the length scale in the *-*direction ( = ) and a as the length scale in the -direction ( = ). = 0 2 ⁺ (N m-3 ) scales the body force. 0Δ 2 is the cross-section average Lorentz body force and

= − is the pressure-induced body force. To assess their relative importance, we define the dimensionless group: = 0∆ 2 .

Next, following HS, we define the dimensionless velocity (, ) and the alternation in the current stream function resulting from the flow ℎ(ξ, η). ⁼ 2 (, ) and ^{= 1 2} (1 ⁺ 2) ^{+ 2 +} 2√ ℎ(, ). Since only differences in matter, we set 1 ^{= 1} 2 Δ and 2 ^{= − 1} 2 Δ, where Δ = 1 − 2. The resulting dimensionless equations are:

$$\begin{pmatrix} A^{2} \frac{\partial^{2}}{\partial \xi^{2}} + \frac{\partial^{2}}{\partial \eta^{2}} & H{a} \frac{\partial}{\partial \eta} \ H{a} \frac{\partial}{\partial \eta} & A^{2} \frac{\partial^{2}}{\partial \xi^{2}} + \frac{\partial^{2}}{\partial \eta^{2}} \end{pmatrix} \begin{pmatrix} u \ h \end{pmatrix} = \begin{pmatrix} -1 \ 0 \end{pmatrix} \tag{6}$$

with the boundary conditions (±1, ) = (, ±1) = ℎ(±1,) ⁼ ℎ(, ±1) ^{= 0}. In the above, ⁼ 0√ is the Hartman number and is the conduit cross-section's aspect ratio.

It is apparent from equation (6) and the boundary conditions that the velocity is an even function of both and : (−, ) = (, ) and (, −) = (, ), while h is odd in and even in : ℎ(, −) = −ℎ(, ), ℎ(−, ) = ℎ(, ), and ℎ(, 0) = 0.

The velocity () satisfying the Dirichlet boundary condition, while (h) satisfies the Neumann boundary condition at = ±1 complicates the analytical solution of equations (6). For small and moderate , equations (6) can be readily solved numerically, e.g., with finite elements (Supplement S1). When the Hartmann number is large, numerical solutions are challenged by thin Hartmann boundary layers of thickness −1 next to = ±1 and Shercliff boundary layers of thickness −1/2 next to = ±1. A fine mesh is required in the boundary layers for appropriate spatial resolution.

The electric potential distribution

$$\phi(x,y) = \frac{1}{2} \left( \phi(-1,\eta) + \phi(1,\eta) \right) + \frac{abP}{\sqrt{\mu\sigma}} \tilde{\varphi}(\xi,\eta). \tag{7}$$

To determine ̃, we integrate the dimensionless 1st order equation (2):

$$\frac{\partial \tilde{\varphi}}{\partial \xi} = -\frac{1}{H_a} \left( \frac{S}{1+S} \right) - \left( \frac{\partial h}{\partial \eta} + H_a u \right). \tag{8}$$

Without loss of generality, we set ̃(0,0) = 0. ̃ is odd in : ̃(−, ) = −̃(, ) and even in : ̃(, −) = ̃(, ). We use asymptotic methods to obtain explicit, approximate relations among the flow rate, pressure gradient, current, and potential difference across the electrodes.

3. Asymptotic Solution for Small Aspect Ratio *A*.

Previously, HS carried out an asymptotic analysis of MHD flow in conduits with rectangular crosssections in the asymptotic limit of → ∞. We are, however, interested in systems wherein the magnetic fields are provided by permanent magnets (typically, 0 < 0.5 T), and conduits' dimensions are on the mm scale. In other words, we focus on moderate values of the Hartmann number (e.g., O(10)). Furthermore, microfluidic systems are often shallow ( ≪ 1). Hence, we examine → 0.

Outer (core) solution: We hypothesize that in the core (outer region), away from the electrodes ( = ±1), both 2 ² and 2ℎ ² are O(1). Thus, in the limit of ^{→ 0} (shallow conduit), we drop all terms with the pre-factor A. Equation (6) reduces to

$$\begin{pmatrix} \frac{\partial^2}{\partial \eta^2} & H_a \frac{\partial}{\partial \eta} \ H_a \frac{\partial}{\partial \eta} & \frac{\partial^2}{\partial \eta^2} \end{pmatrix} \binom{u}{h} = \binom{-1}{0}. \tag{9}$$

The solution of equation (9) is:

$$h_c(\eta) = H_a^{-1} \hat{h}_c(\eta) = -H_a^{-1} \left( \eta - \frac{Sinh(H_a \eta)}{Sinh(H_a)} \right)$$ (10)

and

$$u_c(\eta) = H_a^{-1} \, \hat{u}_c(\eta) = H_a^{-1} Coth(H_a) \left( 1 - \frac{Cosh(H_a \eta)}{Cosh(H_a)} \right). \tag{11}$$

In the above, we defined the scaled variables, ℎ̂ and ̂, for later use.

When ≫ 1, includes the Hartman boundary layers of thickness O( −1) next to = ±1. Away from = ±1, the core flow is uniform ≈ −1. In the limit of small , the velocity profile assumes a parabolic shape ^{1 2} (1 ⁻ 2), consistent with flow driven by uniform body force between two parallel plates. The vertical component of the current density in the core ̃ , = 0 for all Hartmann numbers. The dimensionless, horizontal component of the current density:

$$\tilde{J}{x} = S + H{a} \frac{\partial h}{\partial \eta} (S+1) \tag{12}$$

We use 0−1 as the current density scale instead of the more obvious −1Δ/2 to accommodate opencircuit operation when Δ = 0. In the core,

$$\tilde{J}_{x,c} = (S+1)H_a \frac{Cosh(H_a\eta)}{Sinh(H_a)} - 1$$ (13)

Integration of equation (8) yields the potential distribution in the core.

$$\tilde{\varphi}_c = \left( -\frac{Cosh(H_a)}{Sinh(H_a)} + \frac{1}{H_a} \left( \frac{G}{P} \right) \right) \xi = \left( -Coth(H_a) + \frac{1}{H_a(1+S)} \right) \xi \quad (-1+\delta < \xi < 1-\delta). \quad (14)$$

Witness that ̃ is independent of . represents the thickness of the lateral boundary layers.

As is common in singular perturbations, the core solution does not satisfy the boundary conditions at = ±1. Next to ⁼ ±1, our earlier hypothesis that 2 ² and 2ℎ ² are O(1) cannot be true. We anticipate boundary layers next to = ±1.

Inner region: To resolve the system's behavior next to = −1, we introduce the stretched coordinate () such that = −1 + , (0 ≤ ζ ≤ ∞, δ → 0). Dominant balance (equation 6) suggests, ~ −1/2 . When Ha~O(1), we have a Stokes layer of thickness O(A), wherein the top and bottom surfaces screen the non-slip condition at = −1. When ≫ 1, we have a Shercliff boundary layer of thickness ( −1/2 ).

Rescaling the equations (6) in terms of the stretched coordinate and introducing the scaled (hatted) dependent variables ℎ = −1 ℎ̂ and = −1 ̂, we obtain the boundary layer equations (0 ≤ ζ ≤ ∞, −1 ≤ η ≤ 1):

$$\frac{\partial^2 \hat{u}_i}{\partial \zeta^2} + H_a^{-1} \frac{\partial^2 \hat{u}_i}{\partial \eta^2} + \frac{\partial \hat{h}_i}{\partial \eta} = -1 \tag{15}$$

$$\frac{\partial^2 \hat{h}_i}{\partial \zeta^2} + H_a^{-1} \frac{\partial^2 \hat{h}_i}{\partial \eta^2} + \frac{\partial \hat{u}_i}{\partial \eta} = 0$$ (16)

with the boundary conditions:

$$\hat{u}_i = \hat{u}_i(\zeta, \pm 1) = \hat{u}_i(\infty, \eta) - \hat{u}_c(\eta) \text{ and } \frac{\partial \hat{h}_i(0, \eta)}{\partial \zeta} = \hat{h}_i(\zeta, \pm 1) = \hat{h}_i(\infty, \eta) - \hat{h}_c(\eta) = 0.$$

In the above, the subscripts i and c denote, respectively, inner (boundary layer) and (outer) core dependent variables. The conditions at → ∞ are 1st order asymptotic matches. The solution of the boundary layer equations is independent of the aspect ratio A. For large Ha, to the 1st-order approximation, we could have dropped the ( −1) terms. However, we retain ( −1) terms because we wish to obtain approximations valid for moderate values of the Hartmann number. Serendipitously, ( −1) terms retention provides dissipation that stabilizes the numerical scheme (Supplement S1).

Once the inner solutions are available, we determine the deviation of the boundary layer potential from the linear core potential Δ̃ = ̃ − ̃ by integrating the equation:

$$\frac{\partial \Delta \tilde{\varphi}_i}{\partial \zeta} = \frac{A}{\sqrt{H_a}} \left( \frac{\partial h_c}{\partial \eta} + H_a u_c - \frac{\partial h_i}{\partial \eta} - H_a u_i \right) \tag{17}$$

with ̃(∞, ) = 0, we have

$$\Delta \tilde{\varphi}_i(\zeta) = \frac{A}{Ha^{\frac{1}{2}}} \int{\infty}^{\zeta} \left( \left( H_a^{-1} \frac{\partial \hat{h}_c}{\partial \eta} + \hat{u}_c \right) - \left( H_a^{-1} \frac{\partial \hat{h}_i}{\partial \eta} + \hat{u}_i \right) \right) dl, \tag{18}$$

Since the integrand approaches zero rapidly as increases, it suffices to carry out the integration from a finite > , such that any increase in ^L beyond L,min does not significantly alter ̃(). Typically, ~10 suffices (Supplement S1). The difference at the left wall between the electrode potential extrapolated from the core expression and the actual electrode potential is of interest.

$$\Delta \tilde{\varphi}_{i,L} = \frac{A}{Ha^{\frac{1}{2}}} \Delta \widetilde{\varphi}^*{i,L} = \frac{A}{Ha^{\frac{1}{2}}} \int{\zeta_L}^0 \left( \left( H_a^{-1} \frac{\partial \hat{h}_c}{\partial \eta} + \hat{u}_c \right) - \left( H_a^{-1} \frac{\partial \hat{h}_i}{\partial \eta} + \hat{u}_i \right) \right) dl. \tag{19}$$

̃, is independent of . The result for the left boundary layer ( = −1) also applies to the right boundary layer ( = 1), with ̃, = −̃,.

We numerically solved the boundary layer equations (15-16) (Supplement S1) for the velocity and current stream function.

Fig. 2: The ratio of the computed boundary layer velocity and the core velocity (,) () as a function of the stretched boundary layer coordinate when (A) Ha=5 (purple triangles) and =0.1 (hollow triangle), =0.9 (solid triangle), (B) Ha=10 (pink), and =0.1 (star), =0.9 (7 point star), and (C) Ha=30 (blue), and =0.1 (square), =0.9 (square with cross). (1 − −) is depicted with a solid (green) line.

The approximate velocity in the left lateral boundary layer

$$u(\xi, \eta) \sim u_c(\eta) \left( 1 - e^{-\sqrt{\frac{H_a}{A}}(1+\xi)} \right) \quad (-1 \le \xi \le 1)$$ (20)

agrees well with the numerically computed velocities from the whole field equations (6) and boundary layer equations (15-16) (Fig. 2). The composite expression for the entire flow domain

$$u_{comp}(\xi, \eta) \sim u_c(\eta) \left( 1 - 2e^{-A^{-1}H_a^{1/2}} Cosh\left(A^{-1}H_a^{1/2}\xi\right) \right) (-1 \le \xi \le 1)$$ (21)

satisfies the boundary conditions at ( = ±1) only approximately (with an exponentially small error). Using equation (20), we obtain the approximation for the flow rate () that accounts for the velocity deficits in the lateral boundary layers.

Fig. 3: The relative error (%) between the numerically computed flow rate and the approximate flow rate (equation 22) as a function of the Hartmann number Ha and the aspect ratio A.

The approximate flow rate Qa agrees within 5% with the numerical solution of equations (6) when Ha>5 and 0<A<1. The error decreases as Ha increases and/or A decreases (Fig. 3). This error does not decrease monotonically as the Hartmann number (Ha) increases because our approximate velocity (equation 20) slightly underestimates (overestimates) the actual velocity when < 20 ( > 20) and > 0.9.

Fig. 4: The difference (∆̃, ^∗ ) between the extrapolated core potential next to the left electrode ( = −1) and the actual numerically determined electrode potential as a function of the Hartmann number. The pink asterisks, purple open diamonds, and the solid line correspond, respectively, to ∆̃, ^∗ obtained by numerical integration of the boundary layer equations (13-14 and 18), full-domain numerical simulations (with the core-potential subtracted), and the empirical correlation ∆̃, ^∗ = ^β + , where = 1, = −1, = −0.96 (² = 0.97). The inset on the right depicts the potential in the left boundary layer (solid line) and the extrapolated core potential (dashed line) as functions of when = 5, = 0.05, and = 1.5.

The potential in the lateral boundary layer adjacent to the left electrode (ξ=−1) deviates from its core linear dependence on ξ, dipping below the extrapolated core potential (Fig. 4, inset). The deviation from the extrapolated core potential correlates with the Hartmann number. ∆̃, ^∗ = ^β + , where ~1, ~ − 1, ~ − 0.96 (² = 0.97) (Fig. 4, solid line). The approximate difference () between the left and right electrodes' potential:

$$\psi = \Phi_L - \Phi_R = \frac{2 a b P}{\sqrt{\mu \sigma}} \left( Coth(H_a) + \frac{1}{H_a} \left( -\frac{1}{1+S} + \frac{A}{\sqrt{Ha}} \left( \alpha H_a^{\beta} + c \right) \right) \right). \tag{23}$$

4. Onsager Relations

MHD systems, such as the one described here, can serve multiple functions, including electrical power generators (converting mechanical energy into electrical energy), actuators (converting electrical energy into mechanical motion, such as in pumps), as well as valves, brakes, and sensors. When designing, optimizing, controlling, and operating MHD-based systems, it is convenient to have simple relations among the system's inputs, such as the potential difference between the left and right electrodes and the pressure-gradient , and the system's outputs, such as the current H and the flow rate Q.

$$\binom{\Delta H}{Q} = \begin{pmatrix} K{EE} & K{EH} \ K{HE} & K{HH} \end{pmatrix} \binom{\psi}{G}.$$ (24)

Kij are the Onsager conductivities, where the first subscript denotes the force (E for electrical and H for hydrodynamic), and the second subscript represents the flow (E for electrical current and H for fluid flow). KEE is the electric conductance – the inverse of the resistance. When > 0 and > 0, the flowinduced current < 0 is in the opposite direction to the Ohmic current, and therefore < 0. Based on section 3, we have obtained the Onsager coefficients (Supplement S2):

$$K_{EE} = \frac{A\sigma}{\sqrt{H_a} \left( \sqrt{H_a} \text{Coth}[H_a] + A \left( c + \alpha H_a^{\beta} \right) \right)} \qquad \left( \frac{1}{\Omega m} \right)$$ (25)

$$K_{EH} = -2a^2 \sqrt{\frac{\sigma}{\mu}} \frac{\left(H_a \text{Coth}[H_a] + A\sqrt{H_a}\left(c + \alpha H_a^{\beta}\right) - 1\right)}{H_a^{3/2}\left(\sqrt{H_a} \text{Coth}[H_a] + A\left(c + \alpha H_a^{\beta}\right)\right)} \qquad \left(\frac{m^2 A}{N}\right)$$ (26)

$$K_{HE} = 2a^2 \sqrt{\frac{\sigma}{\mu}} \frac{\left(\sqrt{H_a} - A\left(1 - e^{-\frac{\sqrt{H_a}}{A}}\right)\right) (H_a \operatorname{Coth}[H_a] - 1)}{H_a^2 \left(A\left(c + \alpha H_a^{\beta}\right) + \operatorname{Coth}[H_a]\sqrt{H_a}\right)} \qquad \left(\frac{m^2 A}{N}\right)$$ (27)

$$K_{HH} = \frac{4a^4 \left(\sqrt{H_a} - A\left(1 - e^{-\frac{\sqrt{H_a}}{A}}\right)\right) (H_a Coth[H_a] - 1)}{A \mu H_a^3 \left(A\left(c + \alpha H_a^{\beta}\right) + \sqrt{H_a} Coth[H_a]\right)} \qquad \left(\frac{m^6}{N s}\right)$$ (28)

In the limit of Ha>>1, these Onsager coefficients (25-28) reduce to

$$K_{EE}^{\infty} = \frac{A\sigma}{Ha}, \qquad K{HE}^{\infty} = -K_{EH}^{\infty} = \frac{2a^2}{Ha} \sqrt{\frac{\sigma}{\mu}}, \quad K{HH}^{\infty} = \frac{4a^3b}{\mu H_a^2}$$ (29)

In the limit of A<<1, our Onsager coefficients reduce to

$$K{EE}^{0} = \frac{A\sigma}{H{a}Coth[H{a}]}, \qquad K{HE}^{0} = -K{EH}^{0} = 2\alpha^{2} \sqrt{\frac{\sigma}{\mu} \frac{(H{a}Coth[H{a}] - 1)}{H{a}^{2}Coth[H_{a}]}},$$

$$K{HH}^{0} = \frac{4\alpha^{4}(H{a}Coth[H{a}] - 1)}{\mu A H{a}^{3}Coth[H_{a}]}$$ (30)

In both limits (29) and (30), the coefficients satisfy the Onsager-Casimir reciprocity condition [23, page 90], = −. This anti-symmetry arises because the magnetic field breaks time-reversal symmetry, rendering the classical, symmetric Onsager reciprocity invalid. The Onsager coefficients (26) and (27) satisfy the Onsager-Casimir reciprocity only approximately, likely because of the empirical approximations we used for the velocity and the potential distributions in the lateral boundary layers near ( = ±1). However, the difference between || and || is less than 4% when > 5 and < 1. This difference decreases rapidly as increases and/or decreases (Fig. 5).

Fig 5: The relative difference (%) between the off-diagonal Onsager coefficients 2 × 10² | + − | as a function of and A.

Equations (24) apply to various MHD machines, including electrical power generators that convert mechanical power ( > 0, > 0) into electrical power ( < 0, > 0); machines (e.g., pumps and actuators) that convert electrical power (>0, >0) into mechanical power ( < 0, > 0) as well as hybrid systems combining these functions.

The net power dissipation is ΔΨ + = Ψ² + ² ≥ 0. The off-diagonal terms cancel on account of Onsager-Casimir reciprocity. The term Ψ² represents ohmic losses, and the term ² viscous losses. When the Lorentz body force is equal in magnitude and opposed in direction to the pressure gradient Δ0 2 ^{+ = 0}, there is no flow ^{= 0} at steady-state. This leads to ^{+ 2} − 2 0 = 0.

5. Results and Discussion

Velocity and Electrical Fields

We solved equations (6) for the velocity (, ) and the current stream function ℎ(, ) using finite elements in Mathematica (Supplement S1). Given symmetry, it suffices to solve the equations for the quarter domain {−1 < < 0, 0 < < 1}. The computational domain is divided into the coarse grid region + < and the boundary region, outside the coarse region, to facilitate using a fine mesh in the vicinity of the boundaries without incurring high computational costs. The whole domain numerical results were compared and agreed within 3% with the core (outer) results of Section 3 when √ ≤ 0.28 and −1 + 3 √ ≤ ≤ 1 − 3 √ .

Fig. 6 depicts a sample of our numerical results for the current stream function, velocity, and potential fields. Readers with access to Mathematica can readily solve the MHD equations for other values of and using the supplemental Mathematica Notebook (Supplement S1).

Fig. 6: (a) The current stream function h() as a function of and (=10, A=0.2). (b) The axial velocity u( ) as a function of and (=10, A=0.2). (c) The curves above and below the axis depict, respectively, the velocity u(, 0.5) and the current stream function ℎ(, 0.5) as functions of when =5 (red triangle), 10 (purple star), and 30 (green diamond). A=0.2. The solid and dashed lines show the corresponding core velocity and current stream function. (d) The curves above and below the axis depict, respectively, the velocity u(0, ) and the current stream function h(0, ) as functions of when A=0.2 and =5 (red triangle), 10 (purple star), and 30 (green diamond). The solid and dashed lines show the corresponding core velocity and current stream function. (e) The potential (, ) as a function of and (Ha=10, A=0.2, S=2/3). (f) The potential (, 0.5) (symbol) and () (solid line) as functions of when A=0.2, = 2/3, and = 5 (red triangle), 10 (purple star), and 30 (green diamond).

With the exclusion of the O(Ha -1/2) boundary layers adjacent to ( = ±1), the current stream function is independent of , indicating that in most of the domain, the current flux density (J) is in the -direction (Fig. 6a). Currents in the -direction exist only next to the electrodes ( = ±1).

At moderate and high , the velocity field away from the solid boundaries resembles plug-flow with a nearly uniform velocity −1 ( 0Δ 2 ⁻ ) (Fig. 6b). The induced current produces damping body forces proportional to the velocity's magnitude (Hartmann break), leveling velocity variations away from the non-slip boundaries. As the Hartmann number increases, so does the -component of the current density next to the insulating surfaces = ±1, resulting in thin boundary layers. Although the physical mechanisms and length scales differ, high Hartmann number flow resembles electroosmotic flow in a conduit with a uniform zeta potential at its surface.

Both (, ) and ℎ(, ) attain their core values rapidly as the distance from the electrodes increases (Fig. 6c), even at moderate Hartmann numbers. The numerically computed (, ) and ℎ(, ) (symbols, Fig. 6d) favorably agree with the asymptotic core expressions (, ) and ℎ(, ).

The core region's electrical potential (Fig. 6e) is independent of and varies linearly with (Fig. 6f, solid line). However, the electrodes' potential cannot be obtained by linearly extrapolating the core potential. Witness the nonlinear dependence of on next to the electrode (Figs. 4 inset (solid line), Figs. 6e and 6f ( symbols)).

The Onsager coefficients provide a simple means to determine the performance characteristics of various MHD machines, as illustrated next. Below, we applied Onsager-Casimir reciprocity in our analysis of various MHD machines.

Conversion of Mechanical Power into Electrical Power (Electrical Generators)

Consider pressure-driven flow = (− ) ^&gt; ⁰ in a conduit (Fig. 1). The electrodes are connected to an external load with resistance ′ (Ω m) per conduit's unit length. The generator's internal resistance is ′ = ( ). The liquid metal's motion in the magnetic field induces an electromotive force (emf) >0. = − ′ Δ. varies from zero (short circuit, ′ =0) to open circuit potential ( ′ → ∞).

$$\Delta H = -\frac{\kappa{HE}}{1 + \kappa{HH} R'{ex}} G \text{ and } \Psi = \frac{\kappa{HE} R'{ex}}{1 + \kappa{EE} R'_{ex}} G$$ (31)

When the circuit is open ( ′ ^{→ ∞}), Δ ^{= 0}, and ⁼ . The current can be expressed as a function of the potential difference = ( − ). Maximum current intensity takes place when the electrodes are shortened ( ′ =0, =0), Δ = −. The corresponding S value is = − 2 0.

The flow rate ^{= + =} ( ^{+ 2} ′ 1+ ′ ) is lowest when the electrodes are shortcircuited (=Rex=0, S=Sm). The induced current generates Lorentz body force opposing G throughout the conduit's cross-section (Fig. 7a). As the emf increases (when the external load increases), it produces an Ohmic current opposing the induced current next to the insulating surfaces = ±1, resulting in flowassisting Lorentz force (Fig. 7b, = 0.5 ) and an increased flow rate. This effect is maximized under open circuit conditions (Fig. 7c, = 0) when no current leaves the conduit. The electrical current forms a clockwise eddy when > 0 and a counterclockwise eddy when < 0 (Fig. 8).

Fig. 7: Contour lines of the normalized current stream function ̂ = 0 2 ^{= +} ( + 1)ℎ(, ) when = 5 and = 0.2. = (a, ′ ^{= 0}), ^{= 2} (b), and ^{= 0} (c). The bottom graphs provide a magnified view of the region (−1 < < −0.98; 0.2 < < 0.5) adjacent to the electrode.

Fig. 8: The horizontal component of the normalized current (0, ) as a function of η (0 ≤ ≤ 1) at various Hartmann and numbers: = 5 and = 0 (dashed red line); = 10 and = (purple star), = /2 (six point star), and = 0 (square with cross); and = 100, = 0.2, and = 0 (green line).

The mechanical power consumed per unit length of the conduit is , and the electrical power produced per unit length of the conduit is (−∆). The MHD generator conversion efficiency

$$\varepsilong = -\frac{\Delta H \Psi}{GQ} = \frac{K{HE}^2 R'{ex}}{(1 + K{EE} R'{ex}) \left( K{HE}^2 R'{ex} + K{HH} (1 + K{EE} R'{ex}) \right)}$$ (32)

is independent of the driving pressure gradient () but it depends on the external load. As the resistance ratio = ′ / ′ increases from zero, the generator's efficiency increases, attains a maximum, and then decreases (Fig. 9). As the Hartmann number increases, the optimal resistance ratio and the maximum efficiency , increase and the efficiency's sensitivity to deviations from the optimal r decreases, resulting in a broader peak. In the limit of → ∞,

$$\varepsilong^{\infty} = \frac{\frac{R'{ex}}{R'{in}}}{1 + \frac{R'{ex}}{R'_{in}}}.$$ (33)

Fig 9: MHD generator's efficiency (g) as a function of the external resistance and internal resistance ratio when : = 5, 20, 500 and = 0.05 (Solid line), 0.5 (dashed line), 1 (dotted line).

The maximum efficiency , is achieved when

$$r = \frac{R'{ex}}{R'{in}} = A \,\sigma \frac{\sqrt{K{HH}}}{\sqrt{K{EE}(K{HE}^2 + K{EE}K_{HH})}}$$ (34)

The efficiency at the optimal load

$$\varepsilon{g,m} = \frac{\kappa{HE}^2 + 2K{EE}K{HH} - 2\sqrt{K{HH}}\sqrt{K{EE}(K{HE}^2 + K{EE}K{HH})}}{K{HE}^2}$$ (35)

increases as the Hartmann number increases (Fig. 10). When → ∞, , → 1.

The max efficiency occurs when ( ′ ′ ) ≫1 ~ √. In the range of aspect ratios considered. The maximum efficiency does not vary much with variations in .

Fig 10: MHD generator's efficiency (,) at optimal load as a function of when = 0.05 (solid, purple line), 0.1 (horizontal dashed, blue line), 0.5 (vertical dashed, yellow line), and 1 (dash-dot, red line).

Conversion of Electrical Power into Mechanical Power (Actuator, Pump)

Next, we consider the case when the MHD machine is connected to an external power supply, and current is injected into the liquid metal through the electrodes. The resulting Lorentz body force can propel the liquid metal against an adverse pressure gradient ( < 0) or adverse force. Here, is the pressure gradient induced by the external load. This MHD system can circulate liquid metal for thermal control, displace immiscible fluids, or displace a piston or elastic members for actuation, as well as provide flow control without mechanical valves in fluidic networks [8]. In the absence of flow, the stagnation pressure gradient = − −1 balances the Lorentz forces generated by the injected current.

The pump's mechanical efficiency ^{= −} ∆ varies with the adverse pressure head (G)

$$\varepsilon_p = \frac{G}{G_s} \left( 1 - \frac{G}{Gs} \right) \left( \frac{K{EE} K{HH}}{K{HE}^2} + \frac{G}{G_s} \right)^{-1}$$ (35)

As / increases, the efficiency () increases, attains a maximum, and decreases (Fig. 11). attains its maximum when

$$\frac{G_{op}}{Gs} = \frac{-K{EE}K{HH} + \sqrt{K{EE}K{HH}(K{HE}^2 + K{EE}K{HH})}}{K_{HE}^2}$$ (36)

The optimal conversion efficiency is:

$$\varepsilon{p,m} = \left(\frac{G{op}}{Gs}\right)^2 \frac{K{HE}^2}{K{EE}K{HH}} \tag{37}$$

As the Hartmann number increases, so does the maximum efficiency , (Fig. 11). When → ∞, , → 1 (Fig. 12).

Fig. 11: MHD pump efficiency () as a function of . = 5 (red), 20 (purple), and 500 (green). = 0.05 (solid line), 0.5 (dashed line), and 1 (dotted line).

Fig 12: Optimal pump efficiency (,) as a function of the Hartmann number . = 0.05 (solid, purple line), 0.1 (horizontal dashed, blue line), 0.5 (vertical yellow, solid line), and 1 (orange, dash-dot line).

MHD sensors

Measurement of electrical signals, such as the potential difference across or total current passing through the electrodes, informs the flow rate and the pressure drop in the conduit. Consider, for example, open circuit conditions (∆ = 0). = Ψ and ⁼ (1 ^{+ 2} ) .

MHD Valve and Brake

The induced current in the moving liquid metal generates a Lorentz body force that opposes the flow, providing a (Hartmann) brake. The Lorentz body force can be modulated to further slow the flow by adjusting the magnetic field intensity and/or supplying current from a power supply ( < 0 when > 0). Thus, the MHD machine provides a means to control the flow rate without a mechanical valve [8].

When || is sufficiently large to stop motion, we have a MHD brake. MHD brakes provide several advantages over conventional friction-based brakes. Instead of dissipating the braking energy as heat, the MHD brake provides the opportunity to recover and harvest most of the braking energy. The MHD brake provides contactless, soft braking, avoiding the wear and performance degradation of friction-based brakes.

Conclusions

Recently, interest has grown in applying room- and near-room-temperature liquid metals to miniature systems for power harvesting, actuation, flow control, and sensing. These systems typically operate with permanent magnets at moderate Hartmann numbers. Analyzing such systems often requires repeated numerical solutions of the MHD equations during design, optimization, and control processes. This paper reformulates the relationships between control inputs, such as the electrodes' potential difference () and the pressure gradient (), and outputs, such as flow rate () and current (), in the form of Onsager relations. These simplified algebraic relations facilitate the control and optimization of various MHD machines for pumping, flow control, and actuation when the pressure gradient and electrodes' potentials are independently controlled (Fig. 13 left) and for electrical power generation and sensing when an external load controls electrode potential (Fig. 13 right).

Fig 13: Characterization of the operating regimes of the MHD machine as a function of the pressure gradient and the potential difference across the electrodes (left). ̇ = (Wm-1 ) and ̇ = Δ (Wm-1 ) are, respectively, the mechanical and electrical powers per unit length. is the flow rate (m3 s -1 ),

taken as positive when in the positive z-direction. is taken as positive when in the positive x direction. When operating as a generator (right), the external load resistance ′ (Ω) controls the electrode potential difference .

We explicitly demonstrated that the cross Onsager coefficients satisfy the anti-symmetric Onsager– Casimir reciprocity = −. This anti-symmetry arises because the magnetic field breaks timereversal symmetry, and the classical, symmetric Onsager's reciprocal relations are no longer valid.

To determine the Onsager coefficients, we employed multiple methods, including asymptotic analysis in the limits of small aspect ratio ( → 0), small Hartmann number ( → 0), and large Hartmann number ( → ∞)), as well as full domain, finite element simulations. All approaches yielded consistent results. This favorable agreement serves as a mutual verification of our findings.

The Onsager conductivities can be determined from experimental data, enabling validation of this paper's results. The flow rate (), current (), pressure gradient (), and potential difference across the electrodes ( ) are all measurable quantities. For example, a simple experimental setup could resemble the MHD machines sketched in Fig. 13 (right). In this setup, a duct with a small rectangular crosssectional area, like the one in Fig. 1, connects two relatively large reservoirs containing liquid metal at different levels. The difference in liquid levels controls the pressure gradient. By varying the pressure gradient () and the potential difference ( ) and measuring the flow rate () and the electrical current () repeatedly, one can extract the Onsager conductivities reliably by minimizing an objective (cost) function such as ∑ ((ℓ ⁻ ℓ ) 2 + (ℓ ⁻ ℓ ) 2 ^ℓ ).

Supplementary Material

• S1. Finite Element Calculations of the Velocity and Electrical Fields. A Mathematica Notebook containing the scripts used to solve and visualize the velocity, potential, and current distributions.
• S2. Onsager Coefficients and Performance Characteristics. A Mathematica Notebook presenting Onsager Coefficients and performance calculations of various magnetohydrodynamic (MHD devices.

Funding statement

SES acknowledges support from the USA National Science Foundation Graduate Research Fellowship (DGE-2236662). Any opinions, findings, conclusions, or recommendations expressed in this paper are those of the authors and do not necessarily reflect the views of the National Science Foundation.

References

• [1]L. Zhu, B. Wang, S. Handschuh-Wang, and X. Zhou, "Liquid Metal–Based Soft Microfluidics," Small, vol. 16, no. 9, p. 1903841, 2020, doi: 10.1002/smll.201903841.
• [2]J. Zhong, M. Yi, and H. H. Bau, "Magneto hydrodynamic (MHD) pump fabricated with ceramic tapes," Sens. Actuators Phys., vol. 96, no. 1, pp. 59–66, Jan. 2002, doi: 10.1016/S0924- 4247(01)00764-6.
• [3]A. Lemoff and A. Lee, "An AC magnetohydrodynamic micropump," Sens. Actuators B Chem., vol. 63, pp. 178–185, May 2000, doi: 10.1016/S0925-4005(00)00355-5.
• [4]P. U. Arumugam et al., "Characterization and Pumping: Redox Magnetohydrodynamics in a Microfluidic Channel," J. Electrochem. Soc., vol. 153, no. 12, p. E185, Oct. 2006, doi: 10.1149/1.2352040.
• [5]H. H. Bau, J. Zhong, and M. Yi, "A minute magneto hydro dynamic (MHD) mixer," Sens. Actuators B Chem., vol. 79, no. 2, pp. 207–215, Oct. 2001, doi: 10.1016/S0925-4005(01)00851-6.
• [6]M. Yi, S. Qian, and H. H. Bau, "A magnetohydrodynamic chaotic stirrer," J. Fluid Mech., vol. 468, pp. 153–177, Oct. 2002, doi: 10.1017/S0022112002001635.
• [7]J. West et al., "Application of magnetohydrodynamic actuation to continuous flow chemistry," Lab. Chip, vol. 2, no. 4, pp. 224–230, Nov. 2002, doi: 10.1039/B206756K.
• [8]H. H. Bau, J. Zhu, S. Qian, and Y. Xiang, "A magneto-hydrodynamically controlled fluidic network," Sens. Actuators B Chem., vol. 88, no. 2, pp. 205–216, Jan. 2003, doi: 10.1016/S0925-4005(02)00325- 8.
• [9]S. Holcomb et al., "Oxide-Free Actuation of Gallium Liquid Metal Alloys Enabled by Novel Acidified Siloxane Oils," Langmuir, vol. 32, no. 48, pp. 12656–12663, Dec. 2016, doi: 10.1021/acs.langmuir.6b03501.
• [10]Y. Ding, M. Zeng, and L. Fu, "Surface Chemistry of Gallium-Based Liquid Metals," Matter, vol. 3, no. 5, pp. 1477–1506, Nov. 2020, doi: 10.1016/j.matt.2020.08.012.
• [11]Y. Wang and Y. Xie, "Interfacial interaction-induced super-wettability of gallium-based liquid metals: a review," J. Mater. Chem. A, vol. 12, no. 13, pp. 7396–7417, Mar. 2024, doi: 10.1039/D3TA07297E.
• [12]J.-H. Kim, S. Kim, M. D. Dickey, J.-H. So, and H.-J. Koo, "Interface of gallium-based liquid metals: oxide skin, wetting, and applications," Nanoscale Horiz., vol. 9, no. 7, pp. 1099–1119, Jun. 2024, doi: 10.1039/D4NH00067F.

• [13]J. C. Domínguez-Lozoya, S. Cuevas, D. R. Domínguez, R. Ávalos-Zúñiga, and E. Ramos, "Laboratory Characterization of a Liquid Metal MHD Generator for Ocean Wave Energy Conversion," Sustainability, vol. 13, no. 9, Art. no. 9, Jan. 2021, doi: 10.3390/su13094641.

• [14]H. Teimouri and A. Behzadmehr, "Numerical investigation of a thermosyphon MHD electrical power generator," Energy Convers. Manag., vol. 187, pp. 378–397, May 2019, doi: 10.1016/j.enconman.2019.02.085.

• [15]K. Panchadar, D. West, J. A. Taylor, and T. Krupenkin, "Mechanical energy harvesting using a liquid metal vortex magnetohydrodynamic generator," Appl. Phys. Lett., vol. 114, no. 9, p. 093901, Mar. 2019, doi: 10.1063/1.5078384.

• [16]D. West, J. A. Taylor, and T. Krupenkin, "Alternating current liquid metal vortex magnetohydrodynamic generator," Energy Convers. Manag., vol. 223, p. 113223, Nov. 2020, doi: 10.1016/j.enconman.2020.113223.

• [17]D. Dai, J. Liu, and Y. Zhou, "Harvesting biomechanical energy in the walking by shoe based on liquid metal magnetohydrodynamics," Front. Energy, vol. 6, no. 2, pp. 112–121, Jun. 2012, doi: 10.1007/s11708-012-0186-x.

• [18]L. Xiang, S. Yang, Q. Wang, and J. Wu, "Design and testing of a direct current electromagnetic pump for liquid metal," Therm. Sci. Eng. Prog., vol. 50, p. 102560, May 2024, doi: 10.1016/j.tsep.2024.102560.

• [19]Y. Yerasimou, V. Pickert, B. Ji, and X. Song, "Liquid Metal Magnetohydrodynamic Pump for Junction Temperature Control of Power Modules," IEEE Trans. Power Electron., vol. 33, no. 12, pp. 10583–10593, Dec. 2018, doi: 10.1109/TPEL.2018.2806622.

• [20]P. A. Davidson, An Introduction to Magnetohydrodynamics. in Cambridge Texts in Applied Mathematics. Cambridge: Cambridge University Press, 2001. doi: 10.1017/CBO9780511626333.

• [21]M. Qin and H. H. Bau, "When MHD-based microfluidics is equivalent to pressure-driven flow," Microfluid. Nanofluidics, vol. 10, no. 2, pp. 287–300, Feb. 2011, doi: 10.1007/s10404-010-0668-2.

• [22]J. C. R. Hunt and K. Stewartson, "Magnetohydrodynamic flow in rectangular ducts. II," J. Fluid Mech., vol. 23, no. 3, pp. 563–581, Nov. 1965, doi: 10.1017/S0022112065001544.

• [23]R. (Rolf) Haase, Thermodynamics of irreversible processes. New York : Dover, 1990. Accessed: Feb. 04, 2025. [Online]. Available: http://archive.org/details/thermodynamicsof0000haas

This article (and its supplemental files) may be downloaded for personal use only. Any other use requires prior permission of the author and AIP Publishing.

This article was accepted for publication in Physics of Fluids (DOI: 10.1063/5.0276573).

Supplement S1: Finite Element Simulation of MHD Flow in a Long Conduit with a Rectangular Cross-Section

Supplement to: Sindu Shanmugadas & Haim H. Bau*, 2025, "Onsager Coefficients for Liquid Metal Flow in a Conduit under a Magnetic Field,"

Dept. Mechanical Engineering and Applied Mechanics, University of Pennsylvania, Philadelphia, PA 1904-6315 *Corresponding author: bau@seas.upenn.edu

Mathematica bersion 14.2

We solve the two coupled, dimensionless pdes for the steady, fully-developed axial velocity u(ξ,η) and the current stream function h(ξ,η).

$$\begin{pmatrix} A^{2} \frac{\partial^{2}}{\partial \xi^{2}} + \frac{\partial^{2}}{\partial \eta^{2}} & H{a} \frac{\partial}{\partial \eta} \ H{a} \frac{\partial}{\partial \eta} & A^{2} \frac{\partial^{2}}{\partial \xi^{2}} + \frac{\partial^{2}}{\partial \eta^{2}} \end{pmatrix} \begin{pmatrix} u \ h \end{pmatrix} = \begin{pmatrix} -1 \ 0 \end{pmatrix}$$

$$u(\pm 1, \eta) = u(\xi, \pm 1) = \frac{\partial h(\pm 1, \eta)}{\partial \xi} = h(\xi, \pm 1) = 0.$$

Hence to forth, we refer, respectively, to the normalized coordinates -1<η<1 and -1<ξ<1 as "y" and "x" and the dimensionless axial velocity as v(x,y).

1. Finite Element Solution for a Quarter Domain

In[ ]:= Needs["DifferentialEquationsInterpolatingFunctionAnatomy"]

Exploiting symmetries, we solve the MHD equations for the quarter domain {0∂u(0,y) ^{∂x = ∂}u(x,0) ^∂y =0. h(x,y) is odd with respect to y=0 and even with respect to x=0. h(x,0)= ^∂h(0,y) ^∂x =0. The computational domain is divided into a core (outer) region xn + yn < 0.9 with relatively coarse mesh and an inner region xn + yn > 0.9 with a finer grid.

In[1]:= L1 = ImplicitRegion[-1 ≤ x ≤ 0 && 0 ≤ y ≤ 1, {x, y}];
     refinementRegion := ImplicitRegion[x^2 + y^2 > 0.9, {{x, -1, 0}, {y, 0, 1}}];
     mrfl = With[{rmf = RegionMember[refinementRegion]},
         Function[{vertices, area}, Block[{x, y}, {x, y} = Mean[vertices];
           If[rmf[{x, y}], area > 0.000005, area > 0.0025]]]];

Below, we specify the coupled PDEs and boundary conditions.

op1[A_, Ha_] :=
        Inactive[Div][{{A^2 , 0}, {0, 1}}.Inactive[Grad][v[x, y], {x, y}], {x, y}] + Ha D[h[x, y], y] + 1;
      op2[A_, Ha_] :=
        Inactive[Div][{{A^2 , 0}, {0, 1}}.Inactive[Grad][h[x, y], {x, y}], {x, y}] + Ha D[v[x, y], y];
      BC$Dirichlethq = {v[x, 1]  0, v[-1, y]  0, h[x, 0]  0, h[x, 1]  0};
In[8]:= Ha = 10; A = 0.2;
      {nh, nv} =
       NDSolveValue[{Activate[op2[A, Ha]]  0, Activate[op1[A, Ha]]  0, BC$Dirichlethq}, {h, v},
        Element[{x, y}, L1], Method  {"FiniteElement", "MeshOptions"  {MeshRefinementFunction  mrfl}}]
Out[9]= InterpolatingFunction Domain: {{-1., 0.}, {0., 1.}}
                                        Output: scalar
                                  Data not saved. Save now
                                                            ,
       InterpolatingFunction Domain: {{-1., 0.}, {0., 1.}}
                                        Output: scalar
                                  Data not saved. Save now

Apparently, listing the second PDE first in NDSolve leads to a better conditioned stiffness matrix than otherwise.

2. Data on the FE grid

In[10]:= Show[nv["ElementMesh"]["Wireframe"]] (*Show the finite element mesh *)

In[11]:= nv["ElementMesh"] (* number of elements *)
      NDSolve`FEM`ElementMesh[{{-1., 0.}, {0., 1.}}, {NDSolve`FEM`TriangleElement[<95 256>]}]
In[12]:= nv["ElementMesh"]["MeshOrder"] (* verify that Mathematica uses quadratic elements*)
      2

3. Core (outer) Solution :

In[13]:= vc[Ha_, y_] := 1 / Ha * Cosh[Ha] / Sinh[Ha] (1 - Cosh[Ha * y] / Cosh[Ha])
      hc[Ha_, y_] := -1 / Ha (y - Sinh[Ha * y] / Sinh[Ha])
      jxc[Ha_, y_, S_] := (S + 1) Ha Cosh[Ha y] / Sinh[Ha] - 1 (* current component in the x-direction*)
      fc[Ha_, x_, S_] := (-Coth[Ha] + 1 / Ha / (1 + S)) x (*potential*)

$$\begin{aligned} & & & & & & & & & & & & & & & & & &$$

4. The velocity v(x,y) and the current stream function h(x,y)

In[18]:= gv = Plot3D[ nv[x, y], {x, -1, 0}, {y, 0, 1} , LabelStyle Directive[Medium], AxesLabel {"ξ", "η", "v"}, PlotTheme "Classic"]

In[19]:= gh = Plot3D[Ha nh[x, y], {x, -1, 0}, {y, 0, 1}, AxesLabel {"ξ", "η", "Ha h"}, LabelStyle Directive[Medium], PlotTheme "Classic"]

In[20]:= Plot{Ha nv[0, y], Ha nh[0, y]}, {y, 0, 1}, AxesLabel  {"η"},
       LabelStyle  Directive[Medium], AxesOrigin  {0, 0}, PlotStyle  Thickness[0.01],
       Epilog  Table[{Red, Circle[{y, Ha vc[Ha, y]}, 0.02]}, {y, 0, 1, 0.1}],
         Table[{Green, Circle[{y, Ha hc[Ha, y]}, 0.02]}, {y, 0, 1, 0.1}], Text"Ha v(0,η)",
          {0.3, 0.9}, Background  White, Text"Ha h(0,η)", {0.6, -0.5}, Background  White,
       AspectRatio  Automatic, ScalingFunctions  None, PlotStyle  "Classic"
     (*FE solution (solid line); analytical core solution (symbols) *)
       1.0
                    Ha v(0,η)

In[21]:= Plot[{Ha nv[x, 0], Ha nh[x, 0.3]}, {x, -1, 0}, AxesLabel  {"ξ"},
      LabelStyle  Directive[Medium], AxesOrigin  {0, 0}, PlotStyle  Thickness[0.01],
      Epilog  {Table[{Red, Circle[{x, Ha vc[Ha, 0]}, 0.02]}, {x, -1, 0, 0.1}],
         Table[{Green, Circle[{x, Ha hc[Ha, 0.3]}, 0.02]}, {x, -1, 0, 0.1}], Text["Ha v(ξ,0)",
          {-0.4, 0.9}, Background  White], Text["Ha h(ξ,0.3)", {-0.5, -0.2}, Background  White]},
      AspectRatio  Automatic, ScalingFunctions  None] (*FE solution (solid line);
     analytical core solution (symbols) *)

5. Potential

1 + 1 d p

∂h

• Ha u = -

1

S

= - 1

∂φ

$$\begin{split} & \frac{1}{\partial \xi} = -\frac{1}{H_a} \left( 1 + \frac{1}{P} \frac{1}{dz} \right) - \left( \frac{1}{\partial \eta} + H_a u \right) = -\frac{1}{H_a} \left( \frac{1}{1+S} \right) - \left( \frac{1}{\partial \eta} + H_a u \right) \ & S = \frac{B0 \Delta H}{2 \, a \, G} \ & \ln[s] := \, \mathsf{S} = 2 \, / \, \mathsf{3}; \ & \ln[s] := \, \mathsf{dphidx}[\mathsf{x}-, \mathsf{y}-, \mathsf{S}-] := - \, (\mathsf{Derivative}[\theta, 1] \, [\mathsf{nh}] \, [\mathsf{x}+, \mathsf{y}] + 10 \star \mathsf{nv}[\mathsf{x}+, \mathsf{y}]) \, - \, 1 \, / \, 10 \, (\mathsf{S} \, / \, (1+\mathsf{S}) \, ) \ & \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad$$

• Ha u

∂h

NIntegrate: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

nvphi = Plot3D[phi10[x, y], {x, -1, -0.8}, {y, 0, 1}, PlotTheme "Classic", AxesLabel {"ξ", "η", "φ"}, LabelStyle Directive[Medium]]

NIntegrate: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

Plot comparing core phi (dashed lines) to FEA solutions for Ha=10 (purple)

NIntegrate: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

6. Boundary Layer Problem

∂2u

-1 ^∂2u

To compute the system's behavior in the left boundary layer next to x=-1, we introduce the stretched coordinate ξ -1 + A Ha ζ.

$$\frac{\partial^2 \hat{u}_i}{\partial \varsigma^2} + H_a^{-1} \frac{\partial^2 \hat{u}_i}{\partial \eta^2} + \frac{\partial \hat{u}_i}{\partial \eta} = -1$$

$$\frac{\partial^2 \hat{u}_i}{\partial \varsigma^2} + H_a^{-1} \frac{\partial^2 \hat{u}_i}{\partial \eta^2} + \frac{\partial \hat{u}_i}{\partial \eta} = 0$$

$$\hat{u}_i(0, \eta) = \hat{u}_i(\zeta, 1) = \frac{\partial \hat{u}_i(\zeta, 0)}{\partial \eta} = \hat{u}_i(\infty, \eta) - H_a \hat{u}_c(\eta) = \hat{h}_i(\zeta, 0) = \hat{h}_i(\zeta, 1) = \frac{\partial \hat{u}_i(0, \eta)}{\partial \zeta} = \hat{h}_i(\infty, \eta) - H_a \hat{h}_c(\eta) = 0$$

$$In[]:= \text{lbl} = 12 \text{ ( the BL thickness *); Ha} = 5.;$$

$$\text{BLdomain} = \text{Rectangle}[{0, 0}, {1b1, 1}]; \text{ (* define domain and equations*)}$$

$$\text{refinementRegionBL} := \text{ImplicitRegion}[y > 0.8, {{xx, 0, 1b1}, {y, 0.8, 1}}];$$

$$\text{mrfb} = \text{With}[{\text{rmf} = \text{RegionMember}[\text{refinementRegionBL}]},$$

$$\text{Function}[{\text{vertices, area}}, \text{Block}[{xx, y}, {xx, y} = \text{Mean}[\text{vertices}];$$

$$\text{If}[\text{rmf}[{xx, y}], \text{ area} > 0.00005, \text{ area} > 0.0025]]]];$$

$$\text{be1}[\text{Ha}_] := \text{D}[\text{vb}[xx, y], {xx, 2}] + \text{Ha}^*(-1) \text{D}[\text{vb}[xx, y], {y, 2}] + \text{D}[\text{hb}[xx, y], y] + 1$$

$$\text{be2}[\text{Ha}_] := \text{D}[\text{hb}[xx, y], {xx, 2}] + \text{Ha}^*(-1) \text{D}[\text{hb}[xx, y], {y, 2}] + \text{D}[\text{vb}[xx, y], y]$$

$$\text{bbBC$$Dirichlet}[\text{Ha}_] := {\text{vb}[xx, 1] := 0, \text{vb}[0, y] := 0,$$

$$\text{vb}[\text{lb1}, y] := \text{Ha*vc}[\text{Ha}, y], \text{hb}[xx, 1] := 0, \text{hb}[xx, 0] := 0, \text{hb}[\text{lb1}, y] := \text{Ha*hc}[\text{Ha}, y];$$

$$In[*]:= \text{Show}[\text{RegionPlot}[\text{BLdomain}], \text{ RegionPlot}[\text{refinementRegionBL}], \text{ AspectRatio} \to \text{Automatic}]$$

$$ut[*]:= \text{Show}[\text{RegionPlot}[\text{BLdomain}], \text{ RegionPlot}[\text{refinementRegionBL}], \text{ AspectRatio} \to \text{Automatic}]$$

In[ ]:= {bv, bh} = NDSolveValue[
        {be1[Ha]  0, be2[Ha]  0, bbBC$Dirichlet[Ha]}, {vb, hb}, Element[{xx, y}, BLdomain],
        Method  {"FiniteElement", "MeshOptions"  {MeshRefinementFunction  mrfb}}]
      InterpolatingFunction Domain: {{0., 12.}, {0., 1.}}
                                       Output: scalar
                                 Data not saved. Save now
                                                          ,
       InterpolatingFunction Domain: {{0., 12.}, {0., 1.}}
                                       Output: scalar
                                 Data not saved. Save now
In[ ]:= Plot[{bv[xx, 0.] / Ha / vc[Ha, 0.], bh[xx, 0.5] / Ha / hc[Ha, 0.5]},
       {xx, 0, lbl}, PlotRange  {{0, lbl}, {0, 1.1}}, PlotStyle  Thickness[0.01]]
      0 2 4 6 8 10 12 0.0
      0.4
      0.6
      0.8
      1.0
When Ha=5, the velocity and stream function achieve their asymptotic values within 1% at xx~5. 
In[ ]:= Plot[{bv[1, y], bh[1, y]}, {y, 0, 1}, PlotStyle  Thickness[0.01]]
       0.4

0.2 0.4 0.6 0.8 1.0

Witness the smooth dependence on y.

-0.2

7. Data on the FE grid in the Boundary Layer

In[ ]:= Show[bv["ElementMesh"]["Wireframe"]] (*Show the finite element mesh *)
In[ ]:= bv["ElementMesh"] (* number of elements *)
     NDSolve`FEM`ElementMesh[{{0., 8.}, {0., 1.}}, {NDSolve`FEM`TriangleElement[<57 483>]}]

In[ ]:= bv["ElementMesh"]"MeshOrder"

8. How far is far enough?

2

Since we solve the BL equations numerically, we must define a finite computational domain. We hypothesize that the inner (BL) variables approach their asymptotic core values once lbl exceeds a certain value. The results of the previous section suggest that when Ha>5, the boundary layer variables reach their asymptote value at xx~5. Hence, in most of our calculations,when Ha>5, we use lbl=8.

9. The Perils of Neglecting (Ha -1 ) terms in the BL equations

For low and moderate Hartmann numbers, the terms proportional to Ha -1 cannot be ignored. Although in the limit Ha--> infinity, the .Ha -1 terms could be formally neglected, they provide dissipation to dampen spurious numerical (Gibbs) oscillations. Below, we compare results obtained from integrating the BL equations in the presence and the absence of Ha -1 terms.

In the following BL equations the Ha -1 terms are absent. The solutions still depend on the Ha number through the far field matching conditions.

Without the Ha -1 term, the boundary layer equations become:

In[ ]:= be1h[Ha_] := D[vb[xx, y], {xx, 2}] + D[hb[xx, y], y] + 1
      be2h[Ha_] := D[hb[xx, y], {xx, 2}] + D[vb[xx, y], y]
In[ ]:= {bvh, bhh} =
       NDSolveValue[{be1h[Ha]  0, {be2h[Ha]  0}, bbBC$Dirichlet[Ha]}, {vb, hb}, Element[{xx, y},
          BLdomain],
         Method  {"FiniteElement", "MeshOptions"  {MeshRefinementFunction  mrfb}}]
      InterpolatingFunction Domain: {{0., 12.}, {0., 1.}}
                                         Output: scalar
                                   Data not saved. Save now
                                                             ,
       InterpolatingFunction Domain: {{0., 12.}, {0., 1.}}
                                         Output: scalar
                                   Data not saved. Save now

In[ ]:= Plot[{bvh[0.5, y], bhh[0.5, y]}, {y, 0.9, 1}, PlotStyle Thickness[0.01]]

In the absence of the dissipate terms Ha -1 ^∂2u ∂η² and Ha -1 ^∂2h ∂η² , the solutions exhibit numerical, non-physical oscillations.

11. How well does exponential behavior vC[Ha, y] 1 - e-xx approximates the numerical BL solution for the velocity

The relative difference between the exponential approximation vc(y) (1 - e-x x) and the boundary layer finite element solution as a function of xx when Ha=5, 10, and 30. The error is smaller than 3%.

12. The deviation of the electrode potential Δ φ * from extrapolated core value

$$\Delta \, \tilde{\varphi}_{i,L} = \frac{A}{Ha^{\frac{1}{2}}} \, \Delta \, \tilde{\varphi}^*{i,L} = \frac{A}{Ha^{\frac{1}{2}}} \, \int{\zeta_L}^0 \left( \left( H_a^{-1} \, \frac{\partial \hat{h}_c}{\partial \eta} + \hat{u}_c \right) - \left( H_a^{-1} \, \frac{\partial \hat{h}_i}{\partial \eta} + \hat{u}_i \right) \right) d \, l \, .$$

When ζ^L > lbl, the BL (inner solution denoted with subscript i) matches the core solution (denoted with subscript c), The integrand equals to zero. hc Ha hc and uc Ha uc

In[ ]:= D[hc[Ha, nu], nu] + Ha * vc[Ha, nu] // FullSimplify - 1 Ha ⁺ Coth[Ha] We carry out our calculations without the prefactor A Ha and find ^{Δ φ}* i,L at the electrode's surface In[ ]:= dphidx[xx_, y_, Ha_] := - ¹ Ha ⁺ Coth[Ha]-(Derivative[0, 1][bh][xx, y] / Ha ⁺ bv[xx, y]) In[ ]:= deltaf[l_, y_, Ha_] := NIntegrate[dphidx[xx, y, Ha], {xx, lbl, l}, PrecisionGoal 6, MaxRecursion 50] deltaf0, 0, 10

13. Whole region FE-computed and Approximate Flow Rates

Integrate [2 * vc Ha, y, {y, -1, 1}, {x, -1, 0}] // Simplify (*approximate flow rate*)

Out[*] = $$\frac{4 e^{-\frac{\sqrt{100}}{A}} \left(A + e^{\frac{\sqrt{100}}{A}} \left(-A + \sqrt{Ha}\right)\right) \left(-1 + \text{Ha Coth [Ha]}\right)}{\text{Ha}^{5/2}}$$

$$In[*] := \text{Qapprox [A, Ha]} := \frac{4 \left(1 - \frac{A}{\sqrt{Ha}} \left(1 - e^{-\frac{\sqrt{100}}{A}}\right)\right) \left(-1 + \text{Ha Coth [Ha]}\right)}{\text{Ha}^2}$$

$$In[*] := \text{Qapprox [0.2, 10]}$$

$$Out[] = 0.337232$$ FE whole domain solution for the flow rate. Below, Ha=10 and A= 0.2 $$In[] := \text{NIntegrate 4 * nv[x, y], {y, 0, 1}, {x, -1, 0}}$$

$$Out[*] = 0.338776$$

14. Current Density

-0.83899

The current stream function $$H_z=\frac{H_1+H_2}{2}+\frac{\Delta\,H}{2}\,\eta+a^2\left(\frac{B_0\,\Delta\,H}{2\,a}+G\right)\sqrt{\frac{\sigma}{\mu}}\,h(\xi,\eta)$$ Let $H_1=\frac{\Delta\,H}{2}$ , and $H_2=-\frac{\Delta\,H}{2}$
$$H_z=\frac{\Delta\,H}{2}\,\eta+a^2\left(\frac{B0\,\Delta\,H}{2\,a}+G\right)\,\sqrt{\frac{\sigma}{\mu}}\,h\left(\xi\,,\,\eta\right) = \mathrm{S}\,\mathrm{Let}\,H-\mathrm{a}\,h\left[\mathrm{Let}\,,\,\mathrm{Let}\,\right] \,(\mathrm{S}+1)$$

The electric current is maximized when the device is short circuited (φ 0). Based on our Onsager relations, the max current is: Δ Hmax KE H G -KH E G.

KH E (Supplemental S2) :

$$ln[-]:= \text{ KHE [A_, Ha_]}:= \left(2 \text{ a^2 2} \left(1 + \frac{\text{A}\left(-1 + \text{e}^{-\frac{\sqrt{160}}{\text{A}}}\right)}{\sqrt{\text{Ha}}}\right) \sqrt{\sigma} \left(-1 + \text{Ha Coth [Ha]}\right)\right) / \left(\text{Ha}^{3/2} \sqrt{\mu} \left(\text{A}\left(\text{C} + \text{Ha}^{\beta} \alpha\right) + \sqrt{\text{Ha Coth [Ha]}}\right)\right)$$

For concreteness, we use thermochemical properties of liquid Gallium and specific conduit dimensions. [Sindu: please add units. You did not specify B0]

$$\mu := 2.4 * 10^{\circ} (-3)$$ $\sigma := +2.3 * 10^{\circ} 6$ $a := .001$ $\alpha := 1$ $\beta := -1$ $c := -0.956$ $B0 := Ha/a * Sqrt\mu/\sigma$

At maximum current, $$Sm = -\frac{K{HE}*GB0}{2~a~G}$$ , $${\rm Sm\,[A,\,Ha_]}:= -{\rm KHE\,[A,\,Ha]*B0\,/\,(2~a)}$$

The x-component of current density Jx S + Ha ∂h ∂η (S + 1)

$$J_y = -H_a \, \frac{\partial h}{\partial \xi} \, (S+1)$$

When Ha=10 and A=0.2:

$$ln[*]:= Hz[x{,y{,s}}] := S * y + S 10 * nh[x, y] + 10 * nh[x, y]$$

$ln[*]:= Jx[x{,y{,s}}] := S + 10 * (S + 1) Derivative[0, 1] [nh][x, y]$
$Jy[x{,y{,s}}] := -10 * (S + 1) Derivative[1, 0] [nh][x, y]$

In[ ]:= contourplot = ContourPlot[Hz[x, y, Sm[0.2, 10]], {x, -1, 0}, {y, 0, 1}, PlotLegends  Automatic,
         ContourLabels  True, FrameLabel  {"ξ", "η"}, BaseStyle  {FontSize  21, FontColor  White},
         FrameStyle  Directive[21], ColorFunction  (ColorData["Rainbow"][Rescale[#, {-0.8, 0}]] &)
     , PlotRange  All, ColorFunctionScaling  False];
     vectorfieldsmax = VectorPlot[{Jx[x, y, Sm[0.2, 10]], Jy[x, y, Sm[0.2, 10]]}, {x, -1, 0},
        {y, 0, 1}, VectorPoints  Coarse, VectorSizes  1.3, VectorStyle  Arrowheads[0.05],
        VectorColorFunction  "BrightBands", PlotLegends  Automatic, FrameLabel  {"ξ", "η"}]
     Show[{contourplot, vectorfieldsmax}]
          -1.0 -0.8 -0.6 -0.4 -0.2 0.0
        0.0
        0.2
        0.4
        1.0
              -1
                                                       -0.8
                              -0.6
                                                   -0.4
             -0.2
          0.2
          0.4
          0.6
          0.8
          1.0
       η-0.8
                                                                  -0.6
                                                                  -0.4
                                                                  -0.2

1.0

-1.0 -0.8 -0.6 -0.4 -0.2 0.0

ξ

0.0

Copyright ((2025) Shanmugadas and Bau. This article is distributed under a Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International (CC BY-NC-ND) License. The following article has been accepted by Physics of Fluids. After it is published, it will be found at (DOI: 10.1063/5.0276573).

Supplement S2: Finite Element Simulation of MHD Flow in a Long Conduit with a Rectangular Cross-Section

Supplement to: Sindu Shanmugadas & Haim H. Bau*, 2025, "Onsager Coefficients for Liquid Metal Flow in a Conduit under a Magnetic Field,"

Dept. Mechanical Engineering and Applied Mechanics, University of Pennsylvania, Philadelphia, PA 1904-6315 *Corresponding author: bau@seas.upenn.edu

Onsager Coefficients & Efficiency Calculations

Onsager coefficients are provided in a format suitable for calculations with Mathematica. These coefficients are then used to estimate performance characteristics of various MHD machines.

Ψ - the difference between the (computed) left electrode potential and the electrode potential extrapolated from the core solution.

$$\Delta\,\tilde{\Phi_L} = \frac{A}{\sqrt{H_a}} ! \left(! \alpha H_a{}^\beta + c \right)\,$$

• empirical correlation

$$\tilde{\Phi}{c}\left(\xi\right) = \left(-C \, o \, t \, h \, (H \, a) - \frac{1}{H{a}} \left(\frac{1}{P} \, \frac{d \, p}{d \, z}\right)\right) \xi - \text{core potential as a function of position}$$

$$\Psi = \frac{2 a b P}{\sqrt{\sigma \mu}} \left( \frac{A}{\sqrt{H_a}} \left( \alpha H_a{}^\beta + c \right) + C o t h (H a) + \frac{1}{H_a} \left( \frac{1}{P} \frac{d p}{d z} \right) \right) - \text{the potential difference between the left (x=-1) and right (x=1) electrodes}$$

D H Δ H = H¹ - H2 - total current between the electrodes per conduit's unit length.

The pressure gradient $$G = \frac{-d p}{d z} = -dpdz$$ .

For convenience, in Mathematica expressions, we replaced ξand η with "x" and "y," respectively.

A. Onsager Coefficients

A1. Current (ΔH) Expression

We start with the Ψ expression:

$$\label{eq:problem} \begin{split} &\text{In[110]:=} \ &\text{P[dpdz] := -dpdz + DH * B0 / (2 a)} \ &\text{Psi[Ha, dpdz, A, DH_] :=} \ & \left(2*a*bP[dpdz] / \left(\sqrt{\sigma\mu}\right)\right) \left(\left(A / \sqrt{Ha}\right) (\alpha*Ha^\beta+c) + \text{Coth[Ha]} + 1 / \text{Ha} \left(dpdz / P[dpdz]\right)\right) \ &\text{Psi[Ha, dpdz, A, DH]} \ &\text{Out[112]:=} \end{split}$$

$$\frac{2\,\text{a}\,\text{b}\,\left(\frac{\text{B0\,DH}}{2\,\text{a}}-\text{dpdz}\right)\,\left(\frac{\text{dpdz}}{\left(\frac{\text{B0\,DH}}{2\,\text{a}}-\text{dpdz}\right)\,\text{Ha}}+\frac{\text{A}\left(-\text{0.956}+\frac{1}{\text{Nb}}\right)}{\sqrt{\text{Ha}}}+\text{Coth}\left[\text{Ha}\right]\right)}{\sqrt{\mu\,\sigma}}$$

The above expression for Psi is simplified and the coefficients of G and Ψ are identified;

$$\Psi = a{11} \left( -\frac{dP}{dz} \right) + a{12} DH$$ In[113]:=

a11[A_, Ha_] := $$\frac{2 \text{ a b} \left(-\frac{1}{\text{Ha}} + \frac{\text{A} \left(\text{c-Ha}^{\beta} \cdot \text{a}\right)}{\sqrt{\text{Ha}}} + \text{Coth} \left[\text{Ha}\right]\right)}{\sqrt{\mu \, \sigma}}$$

a12[A_, Ha_] := $$\frac{\text{b BO } \left(\text{A } \left(\text{c} + \text{Ha}^{\beta} \, \alpha\right) + \sqrt{\text{Ha}} \, \text{Coth} \left[\text{Ha}\right]\right)}{\sqrt{\text{Ha}} \, \sqrt{\mu \, \sigma}}$$

In[115]:=

Psi[Ha, dpdz, A, DH]-(DH*a12[A, Ha] + (-dpdz)*a11[A, Ha]) // FullSimplify (*verify that the coefficient form is equal to original expression*)

The expressions are rearranged to isolate DH.

$$DH = \frac{1}{a{12}} \psi + \frac{-a{11}}{a{12}} \left( -\frac{dP}{dz} \right)$$ $$DH = K{EE} \psi + K_{EH} \left( -\frac{dP}{dz} \right)$$

0.

In[116]:=

KEE[A_, Ha_] := $$\frac{A \sigma}{\sqrt{\text{Ha} (A (c + \text{Ha}^{\beta} \alpha) + \sqrt{\text{Ha} Coth[Ha]})}}$$

KEH [A_, Ha_] := $$-\left(\left(2\,a^2\,\sqrt{\sigma/\mu}\,\left(-1+A\,c\,\sqrt{Ha}\,+A\,Ha^{\frac{1}{2}+\beta}\,\alpha+Ha\,Coth\,[Ha]\right)\right)\right/\,\left(Ha^{\,\alpha}\,(3/2)\,\left(A\,\left(c+Ha^{\,\beta}\,\alpha\right)\,+\,\sqrt{Ha}\,Coth\,[Ha]\right)\right)\right)$$

In[118]:=

DH[Ha_, A_, G_, ψ_] := KEE[A, Ha]*ψ + KEH[A, Ha]*(G) (*DH expressed with Onsager Coefficients*)

A2. Flow Rate (Q) expression

The approximate velocity field that accounts for the velocity deficit in the boundary layer next to x=-1:

$$u(\xi, \eta) \approx u_c(\eta) \left(1 - e^{\left(-\frac{\sqrt{H\alpha}}{A}(1+\xi)\right)}\right)$$

In[119]:=

uc[y_] := 1/Ha*CothHa;

Integrate[2 ucy, {x, -1, 0}, {y, -1, 1}] // FullSimplify (* integral of the velocity field over the dimensionless cross-section *)

$$\frac{4\left(1+\frac{A\left(-1+e^{\frac{\sqrt{G_a}}{A}}\right)}{\sqrt{Ha}}\right)\,\left(-1+\text{Ha Coth}\left[\text{Ha}\right]\right)}{\text{Ha}^2}$$

$$4\left(1+\frac{A\left(-1+e^{-\frac{\sqrt{16}}{A}}\right)}{\sqrt{\text{Ha}}}\right)\left(-1+\text{Ha Coth[Ha]}\right)$$

$$Q[\text{Ha, dpdz, A, DH}]:=\frac{}{\text{Ha}^{2}}\text{a}^{3}\text{b}*\text{P[dpdz]}/\mu$$

The above expression for the dimensional flow-rate is rewritten as

$$Q = a{21} \Delta H + a{22} \left( -\frac{dp}{dz} \right)$$ . Below, we replace $\Delta H$ with the expression from section A1.

a21[A_, Ha_] := $$\frac{1}{\text{Ha}^2 \, \mu}$$ 2 a² b B0 $\left[1 + \frac{\text{A}\left(-1 + \text{e}^{-\frac{\sqrt{\text{Ha}}}{\Lambda}}\right)}{\sqrt{\text{Ha}}}\right]$ (-1 + Ha Coth[Ha])

a22 [A_, Ha_] := $$\frac{1}{\text{Ha}^2 \, \mu} \, 4 \, \text{a}^3 \, \text{b} \left( 1 + \frac{\text{A} \left( -1 + \text{e}^{-\frac{\sqrt{\text{Ha}}}{\text{A}}} \right)}{\sqrt{\text{Ha}}} \right) \, (-1 + \text{Ha Coth [Ha]})$$

In[124]:=

Q[Ha, dpdz, A, DH]-(DH*a21[A, Ha] + (-dpdz)*a22[A, Ha]) // FullSimplify (*verify that the coefficient form is equal to original expression*)

Out[124]=

0

$$Q = K{HE} \Psi + K{HH} \left( -\frac{d p}{d z} \right)$$

The simplified coefficients are below:

In[125]:=

KHE[A_, Ha_] :=

$$\left(2 \text{ a}^2 \left(1 + \frac{\text{A}\left(-1 + \text{e}^{-\frac{\sqrt{16}}{\Lambda}}\right)}{\sqrt{\text{Ha}}}\right) \sqrt{\sigma} \left(-1 + \text{Ha Coth[Ha]}\right)\right) / \left(\text{Ha}^{3/2} \sqrt{\mu} \left(\text{A}\left(\text{c} + \text{Ha}^\beta \alpha\right) + \sqrt{\text{Ha Coth[Ha]}}\right)\right)$$

$$\text{KHH}\left[\text{A}_\text{, Ha}_\right] := \left(\text{4 a}^4 \left(\text{1} + \frac{\text{A}\left(-\text{1} + \text{e}^{-\frac{\sqrt{\text{Ha}}}{\text{A}}}\right)}{\sqrt{\text{Ha}}}\right) \left(-\text{1} + \text{Ha}\left(\text{Coth}\left[\text{Ha}\right]\right)\right) \right) \left(\text{Ha}^{5/2} \text{ A } \mu \left(\text{A}\left(\text{c} + \text{Ha}^\beta \alpha\right) + \sqrt{\text{Ha}}\left(\text{Coth}\left[\text{Ha}\right]\right)\right) \right) \right)$$

Q[Ha_, A_, G_, $$\psi_{-}$$ ] :=
KHE[A, Ha] * $\psi$ + KHH[A, Ha] *G (*The flow rate Q expressed with Onsager Coefficients*)

A3. Onsager Coefficient Matrix

$$\begin{pmatrix} K{\rm EE} & K{\rm EH} \ K{\rm HE} & K{\rm HH} \end{pmatrix} \begin{pmatrix} \Psi \ G \end{pmatrix}$$ In[128]:= onsmat := {{KEE[A, Ha], KEH[A, Ha]}, {KHE[A, Ha], KHH[A, Ha]}}

In[20]:= onsmat // MatrixForm

Out[20]//MatrixForm=

$$\left( \begin{array}{c} A \, \sigma \ \hline \sqrt{\text{Ha}} \, \left( A \, \left( c + \text{Ha}^{\theta} \, \alpha \right) + \sqrt{\text{Ha}} \, \text{Coth} \left[ \text{Ha} \right] \right)} \ \ 2 \, a^{2} \left( 1 + A \, \left( \sqrt{\text{Ha}} + A \, \text{Ha}^{\frac{1}{2} + \beta} \, \alpha + \text{Ha} \, \text{Coth} \left[ \text{Ha} \right] \right)} \ \hline + A \, \left( 1 + A \, \left( \sqrt{\text{Ha}} + A \, \text{Ha}^{\frac{1}{2} + \beta} \, \alpha + \text{Ha} \, \text{Coth} \left[ \text{Ha} \right] \right)} \ \hline + A \, \left( 1 + A \, \left( \sqrt{\text{Ha}} + A \, \text{Ha}^{\frac{1}{2} + \beta} \, \alpha + \text{Ha} \, \text{Coth} \left[ \text{Ha} \right] \right)} \ \hline + A \, \left( 1 + A \, \left( \sqrt{\text{Ha}} + A \, \text{Ha}^{\frac{1}{2} + \beta} \, \alpha + \text{Ha} \, \text{Coth} \left[ \text{Ha} \right] \right)} \ \hline + A \, \left( 1 + A \, \left( \sqrt{\text{Ha}} + A \, \text{Ha}^{\frac{1}{2} + \beta} \, \alpha + \text{Ha} \, \text{Coth} \left[ \text{Ha} \right] \right)} \ \hline + A \, \left( 1 + A \, \left( \sqrt{\text{Ha}} + A \, \text{Ha}^{\frac{1}{2} + \beta} \, \alpha + \text{Ha} \, \text{Coth} \left[ \text{Ha} \right] \right)} \ \hline + A \, \left( 1 + A \, \left( \sqrt{\text{Ha}} + A \, \text{Ha}^{\frac{1}{2} + \beta} \, \alpha + \text{Ha} \, \text{Coth} \left[ \text{Ha} \right] \right)} \ \hline + A \, \left( 1 + A \, \left( \sqrt{\text{Ha}} + A \, \text{Ha}^{\frac{1}{2} + \beta} \, \alpha + \text{Ha} \, \text{Coth} \left[ \text{Ha} \right] \right)} \ \hline + A \, \left( 1 + A \, \left( \sqrt{\text{Ha}} + A \, \text{Ha}^{\frac{1}{2} + \beta} \, \alpha + \text{Ha} \, \text{Coth} \left[ \text{Ha} \right] \right)} \ \hline + A \, \left( 1 + A \, \left( \sqrt{\text{Ha}} + A \, \text{Ha}^{\frac{1}{2} + \beta} \, \alpha + \text{Ha} \, \text{Coth} \left[ \text{Ha} \right] \right)} \ \hline + A \, \left( 1 + A \, \left( \sqrt{\text{Ha}} + A \, \text{Ha}^{\frac{1}{2} + \beta} \, \alpha + \text{Ha} \, \text{Coth} \left[ \text{Ha} \right] \right)} \right)$$

A4. Onsager-Casimir Reciprocity

Since we use approximate expressions for the velocity and potential in the boundary layers next to the electrodes, we cannot demonstrate reciprocity symbolically except in limiting cases (see manuscript). Instead, we show numerically that KE H + KH E ≈ 0.

To verify that the O-C reciprocity is approximately satisfied, we compute (KEH+KHE)/(KEH-KHE) as a function of Ha and A.

In[21]:= (KEH[A, Ha] + KHE[A, Ha])/(KEH[A, Ha] - KHE[A, Ha]) // FullSimplify

$$-\frac{\left[1+\frac{\sqrt{\alpha{}}}{\sqrt{\text{Hz}}}\right] \sqrt{\sigma} \ \left(-1+\text{Ha Coth [Ha]}\right)}{\sqrt{\mu}} - \sqrt{\frac{\sigma}{\mu}} \ \left(-1+\text{A C }\sqrt{\text{Ha}} \ + \text{A Ha}^{\frac{1}{2}+\beta} \ \alpha + \text{Ha Coth [Ha]}\right)}{\sqrt{\mu}} - \frac{\left[1+\frac{\sqrt{\alpha{}}}{\sqrt{\pi_{}}}\right]}{\sqrt{\sigma} \ \left(-1+\text{Ha Coth [Ha]}\right)} + \sqrt{\frac{\sigma}{\mu}} \ \left(-1+\text{A C }\sqrt{\text{Ha}} \ + \text{A Ha}^{\frac{1}{2}+\beta} \ \alpha + \text{Ha Coth [Ha]}\right)}{\sqrt{\mu}} + \sqrt{\frac{\sigma}{\mu}} \ \left(-1+\text{A C }\sqrt{\text{Ha}} \ + \text{A Ha}^{\frac{1}{2}+\beta} \ \alpha + \text{Ha Coth [Ha]}\right)$$

The above expression simplifies to

In[129]:=

$$\begin{aligned} \text{diffons} := & - \left[ \left( A \left( 1 - e^{-\frac{\sqrt{||a|}}{A}} - Ha \left( c + Ha^{\beta} \; \alpha \right) + \left( -1 + e^{-\frac{\sqrt{||a|}}{A}} \right) Ha \; \text{Coth} \left[ Ha \right] \right) \right] / \ \ \left( \sqrt{Ha} \left( -1 + A \; c \; \sqrt{Ha} \; + A \; Ha^{\frac{1}{2} + \beta} \; \alpha + Ha \; \text{Coth} \left[ Ha \right] + \left( 1 + \frac{A \left( -1 + e^{-\frac{\sqrt{||a|}}{A}} \right)}{\sqrt{Ha}} \right) \; \left( -1 + Ha \; \text{Coth} \left[ Ha \right] \right) \right] \right) \right] \end{aligned}$$

The following contour plot shows the relative difference (%) between the off-diagonal Onsager coefficients: 2 *102 KE H ⁺ KH E KE H - KH E

In[173]:= $$\alpha := 1$$ $\beta := -1$ $c := -0.956$

ContourPlot2*10 ^ 2 Abs[diffons], {A, 0, 1}, {Ha, 5, 20}, PlotRange All, ColorFunction "DarkBands", Contours {0.1, 0.5, 1, 2}, FrameLabel "A", "Ha", BaseStyle {FontSize 15, Black}, PlotLegends Automatic, ContourLabels True

When A<0.8 and Ha>5, the O-C reciprocity is satisfied with less than 2% difference.

B Efficiency Calculations

B1. MHD Generator

Consider pressure-driven flow in a conduit (Fig. 1). The electrodes are connected to an external load with resistance Rx (Ω ^m) per conduit's unit length.

The internal resistance Ri n 1 σ A (^{Ω m}). The conductor's motion in the magnetic field induces an electromotive force (EMF) Ψ>0.

The generator efficiency is the ratio of the electrical power produced and the mechanical power consumed. ψ -R x * D H .

$$\epsilon_g = \frac{-\Delta H \Psi}{Q G} = \frac{\Delta H^2 R x}{Q G}$$

We rewrite DH and Q in terms of Rx instead of ψ. We rename the expressions for the flow rate and current Qg and DHg;

Qg[Ha_, A_, G_, Rx_] := (KHH[A, Ha]-(KHE[A, Ha]*KEH[A, Ha]*Rx)/(1 + KEE[A, Ha]*Rx)) (G)

DHg[Ha_, A_, G_, Rx_] := (KEH[A, Ha]*G)/(1 + KEE[A, Ha]*Rx)

The generator efficiency is:

In[30]:= (DHg[Ha, A, G, Rx]^ 2 Rx)/(Qg[Ha, A, G, Rx] *G) // Simplify Rx ^{σ A} -0.956 ^{+ 1} Ha ⁺ Ha Coth[Ha] 2 -1 ^{+ A} (1. - 0.956 Ha) Ha + Ha Coth[Ha] 2 1 + ^A -1 ⁺ - Ha Ha A 1. - 0.956 Ha + 1. Ha Rx σ Ha ⁺ 1. Ha Coth[Ha] 2 (-1 + Ha Coth[Ha])

$$\left[ \left( -\text{0.956} + \frac{\textbf{1}}{\text{Ha}} \right) \sqrt{\text{Ha}} + \frac{\text{Ha Coth} \left[ \text{Ha} \right]}{\text{A}} + \frac{\text{Rx } \sqrt{\mu} \sqrt{\sigma} \sqrt{\frac{\sigma}{\mu}} \left( -\text{1} + \frac{\text{A } (\text{1.-0.956 Ha})}{\sqrt{\text{Ha}}} + \text{Ha Coth} \left[ \text{Ha} \right] \right)}{\text{1} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}} \left( \text{A} \left( -\text{0.956} + \frac{\textbf{1}}{\text{Ha}} \right) + \sqrt{\text{Ha}} \text{ Coth} \left[ \text{Ha} \right] \right)}} \right) \right] \right) + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}} \left( \text{A} \left( -\text{0.956} + \frac{\textbf{1}}{\text{Ha}} \right) + \sqrt{\text{Ha}} \text{ Coth} \left[ \text{Ha} \right] \right)}}{\text{1} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{Ha}}} + \frac{\text{A Rx } \sigma}{\sqrt{\text{$$

Introducing Ri and the resistance ratio RxRi = Rex/Rin, the expression for efficiency becomes:

$$\ln[31] := -\left(\left(e^{\frac{\sqrt{Ha}}{A}} \sqrt{Ha} \ RxRi \left(-1 + A \ C \ \sqrt{Ha} + A \ Ha^{\frac{1}{2} + \beta} \ \alpha + Ha \ Coth \ [Ha]\right)^2\right) / \left(\left(-A + e^{\frac{\sqrt{Ha}}{A}} \left(A - \sqrt{Ha}\right)\right) \ (1 + RxRi) \right)$$

$$\left(-1 + Ha \ Coth \ [Ha]\right) \left(A \left(C \ \sqrt{Ha} + Ha^{\frac{1}{2} + \beta} \ \alpha\right) + RxRi + Ha \ Coth \ [Ha]\right) \right) / / Full Simplify$$

$$\frac{e^{\frac{\sqrt{NB}}{A}} \sqrt{\text{Ha}} \ \text{RxRi} \left(-1 + \frac{A \ (1, -0.956 \ \text{Ha})}{\sqrt{\text{Ha}}} + \text{Ha} \ \text{Coth} \left[ \ \text{Ha} \right] \right)^2}{\left(-A + e^{\frac{\sqrt{NB}}{A}} \left(A - \sqrt{\text{Ha}} \ \right) \right) \ \left(1 + \text{RxRi} \right) \ \left(-1 + \text{Ha} \ \text{Coth} \left[ \ \text{Ha} \right] \right) \ \left(\frac{A \ (1, -0.956 \ \text{Ha})}{\sqrt{\text{Na}}} + \text{RxRi} + \text{Ha} \ \text{Coth} \left[ \ \text{Ha} \right] \right)}{\sqrt{\text{Na}}} \right)}$$

The above is equivalent to:

$$\begin{aligned} & \mathsf{Effg}[\mathsf{Ha}{-}, \, \mathsf{A}{-}, \, \mathsf{RxRi}_{-}] \, := \, \bigg( \, \mathsf{RxRi} \, \left( -1 + \mathsf{A} \, \mathsf{C} \, \sqrt{\mathsf{Ha}} \, + \mathsf{A} \, \mathsf{Ha}^{\frac{1}{2} + \beta} \, \alpha + \mathsf{Ha} \, \mathsf{Coth} \, [\mathsf{Ha}] \, \right)^2 \bigg) \bigg/ \ & \left( \left( 1 + \frac{\mathsf{A} \left( -1 + \mathsf{e}^{-\frac{\sqrt{\mathsf{Ha}}}{\mathsf{A}}} \right)}{\sqrt{\mathsf{Ha}}} \right) \, \left( -1 + \mathsf{Ha} \, \mathsf{Coth} \, [\mathsf{Ha}] \right) \, \left( \, \mathsf{A} \, \left( \mathsf{C} \, \sqrt{\mathsf{Ha}} \, + \mathsf{Ha}^{\frac{1}{2} + \beta} \, \alpha \right) + \mathsf{RxRi} + \mathsf{Ha} \, \mathsf{Coth} \, [\mathsf{Ha}] \right) \, \left( 1 + \mathsf{RxRi} \right) \end{aligned} \right)$$

Below, the MHD generator's efficiency is depicted as a function of the resistances ratio RxRi for various Ha and A:

$$\begin{split} & \ln |38| = \text{Plot} \Big[ { \text{Effg[5, 0.05, RxRi], Effg[5, 0.5, RxRi], Effg[5, 1, RxRi], Effg[20, 0.05, RxRi], } \ & \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad$$

Next, we find the resistance ratio that maximizes MHD generator's efficiency:

In[44]:= Solve[D[Effg[Ha, A, RxRi], RxRi] 0, RxRi]

$$\left{\left{\text{RxRi}\rightarrow\text{-1. }\sqrt{\frac{\text{1. A}}{\sqrt{\text{Ha}}}}-\text{0.956 A }\sqrt{\text{Ha}}+\text{1. Ha Coth [Ha]}\right.\right},\ \left{\text{RxRi}\rightarrow\sqrt{\frac{\text{1. A}}{\sqrt{\text{Ha}}}}-\text{0.956 A }\sqrt{\text{Ha}}+\text{1. Ha Coth [Ha]}\right.\right}$$

$$In[45] = Effg[Ha, A, \sqrt{Ac \sqrt{Ha} + AHa^{\frac{1}{2}+\beta}} \alpha + Ha Coth[Ha]] // FullSimplify$$

$$\left( \left( -1 + \frac{A (1. - 0.956 \text{ Ha})}{\sqrt{\text{Ha}}} + \text{Ha Coth [Ha]} \right)^2 \sqrt{\frac{A (1. - 0.956 \text{ Ha})}{\sqrt{\text{Ha}}}} + \text{Ha Coth [Ha]} \right) \right/$$

$$\left[ \left( 1 + \frac{\mathsf{A}\left(-1 + e^{-\frac{\sqrt{\mathsf{IG}}}{\mathtt{A}}}\right)}{\sqrt{\mathsf{Ha}}} \right) \left(-1 + \mathsf{Ha}\,\mathsf{Coth}\,[\mathsf{Ha}]\right) \left( 1 + \sqrt{\frac{\mathsf{A}\left(1 \cdot - 0.956\,\mathsf{Ha}\right)}{\sqrt{\mathsf{Ha}}}} + \mathsf{Ha}\,\mathsf{Coth}\,[\mathsf{Ha}] \right) \right] \right] = \left( -\frac{\mathsf{A}\left(-1 + e^{-\frac{\sqrt{\mathsf{IG}}}{\mathtt{A}}}\right)}{\sqrt{\mathsf{Ha}}} + \frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\sqrt{\mathsf{IG}}}{\mathtt{A}}}\right)}{\mathsf{A}}}\right)}{\sqrt{\mathsf{Ha}}} \right) \right) = \left( -\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-}}\right(-1 + e^{-\frac{\mathsf{A}\left(-1 + e^{-}\right(-1 + e^{$$

$$\left( \frac{ \mathsf{A} \ ( \texttt{1.} - \texttt{0.956} \, \mathsf{Ha} ) }{ \sqrt{\mathsf{Ha}} } \ + \, \mathsf{Ha} \, \mathsf{Coth} \, [\, \mathsf{Ha} \, ] \ + \ \sqrt{ \frac{ \mathsf{A} \ ( \texttt{1.} - \texttt{0.956} \, \mathsf{Ha} ) }{ \sqrt{\mathsf{Ha}} } } \ + \, \mathsf{Ha} \, \mathsf{Coth} \, [\, \mathsf{Ha} \, ] \ \right)$$

Further simplification leads to the max generator efficiency Effgmax:

$$\ln[52] := \operatorname{Effgmax}[\operatorname{Ha}{-}, \operatorname{A}{-}] := \left( \mathbf{1} + \operatorname{Ac} \sqrt{\operatorname{Ha}} + \operatorname{AHa}^{\frac{1}{2} + \beta} \alpha + \operatorname{Ha} \operatorname{Coth}[\operatorname{Ha}] - 2 \sqrt{\operatorname{Ac} \sqrt{\operatorname{Ha}} + \operatorname{AHa}^{\frac{1}{2} + \beta} \alpha + \operatorname{Ha} \operatorname{Coth}[\operatorname{Ha}]} \right) / \left( \left( \mathbf{1} + \frac{\operatorname{A} \left( -\mathbf{1} + \mathbf{e}^{-\frac{\sqrt{\operatorname{Ha}}}{A}} \right)}{\sqrt{\operatorname{Ha}}} \right) (-\mathbf{1} + \operatorname{Ha} \operatorname{Coth}[\operatorname{Ha}]) \right)$$

The MHD generator's efficiency at optimal resistance ratio as a function of the Hartman number for various aspect ratios A.

B2. MHD Pump

Next, consider the case when current is injected into the electrodes. As the result of the Lorentz body force, the liquid metal is propelled against an adverse pressure gradient G<0. The maximum flow rate is obtained in the absence of adverse pressure (G=0)

Efficiency = $$\frac{-QG}{\Delta H \psi}$$ . The stagnation pressure gradient $Gs = \frac{-K{HE} \Psi}{K_{HH}}$ balances the Lorentz forces generated by the injected current. We rewrite the

efficiency expression in terns of G Gs

In[141]:=

dh := kee*Ψ + keh (G) q := khe*Ψ + khh*G

$$Gs := -\frac{khe \Psi}{khh}$$

$$In[60]:= (-qG) / (dh \Psi) // FullSimplify$$

$$-\frac{\mathsf{G}\;(\mathsf{G}\;\mathsf{khh}\;+\;\mathsf{khe}\;\Psi)}{\Psi\;(\mathsf{G}\;\mathsf{keh}\;+\;\mathsf{kee}\;\Psi)}$$

$$G/Gs (1-G/Gs) \left(\frac{khh (kee)}{khe^2} + G/Gs (1-(keh+khe)/(khe))\right)^{(-1)} // FullSimplify$$

(*Verify that we get the same pump efficiency expression using the G/Gs term*)

Out[145]=

$$-\frac{\mathsf{G}\;\left(\mathsf{G}\;\mathsf{khh}\;+\;\mathsf{khe}\;\Psi\right)}{\Psi\;\left(\mathsf{G}\;\mathsf{keh}\;+\;\mathsf{kee}\;\Psi\right)}$$

This pump efficiency:

$$\epsilon_p = \frac{G}{G_s} \left( 1 - \frac{G}{G_s} \right) \left( b_1 + \frac{G}{G_s} \left( 1 - b_2 \right) \right)^{-1}$$

Gst ^{a l l}

In[62]:= b1[A_, Ha_] := KEE[A, Ha]*KHH[A, Ha]/(KHE[A, Ha]^ 2) // FullSimplify

In[79]:= FullSimplify[b2[A, Ha] := (KEH[A, Ha] + KHE[A, Ha])/KHE[A, Ha]]

In[80]:= FullSimplify[b2[A, Ha], Assumptions $$\rightarrow \mu > 0$$ ] // Chop

$$\frac{\text{A (1. -1. Ha Coth [Ha])} + e^{\frac{\sqrt{168}}{A}} \; (-0.956 \, \text{A Ha} + 1. \, \text{A Ha Coth [Ha]})}{\left(-1. \, \text{A} + e^{\frac{\sqrt{168}}{A}} \; \left(1. \, \text{A} - 1. \; \sqrt{\text{Ha}} \; \right) \right) \; (-1. \, + \text{Ha Coth [Ha]})}$$

$$\label{eq:lossenergy} \begin{split} & \ln[81] := \ b2\left[A{_}, \ Ha\right] \ := 1 - \frac{-1 + A \ c \ \sqrt{Ha} \ + A \ Ha^{\frac{1}{2} + \beta} \ \alpha + Ha \ Coth \left[Ha\right]}{\left(1 + \frac{A \left(-1 + e^{-\frac{\sqrt{Ma}}{A}}\right)}{\sqrt{Ha}}\right) \ \left(-1 + Ha \ Coth \left[Ha\right]\right) \end{split}$$

In[82]:= effpggs[Ha_, A_, GGS_, Ψ_] := GGS (1 - GGS) (b1[A, Ha] + GGS (1 - b2[A, Ha]))^(-1)

In[73]:= Plot{effpggs[5, 0.05, GGS, 0.3], effpggs[5, 0.5, GGS, 0.3], effpggs[5, 1, GGS, 0.3], effpggs[20, 0.05, GGS, 0.3], effpggs[20, 0.5, GGS, 0.3], effpggs[20, 1, GGS, 0.3], effpggs[500, 0.05, GGS, 0.3], effpggs[500, 0.5, GGS, 0.3], effpggs[500, 1, GGS, 0.3]}, {GGS, 0, 1}, PlotStyle {Thick, Dashed, Dotted}, PlotTheme "ThickLines", AxesLabel " ^G ", "ϵm",

LabelStyle Directive[Medium], PlotRange All, PlotLegends Automatic

Below, we calculate the pressure gradient ratio that optimizes the pump/actuator efficiency as well as the optimal efficiency.

$$ln[93]:=$$ Solve[D[effpggs[Ha, A, GGS, $\Psi$ ], GGS] $==$ 0, GGS, Assumptions $\rightarrow \mu > 0$ ] // FullSimplify

$$\left{\left{\mathsf{GGS} \to \frac{1}{1-\sqrt{\left(\mathsf{A} \,\mathsf{C} \,\sqrt{\mathsf{Ha}} \,+\, \mathsf{A} \,\mathsf{Ha}^{\frac{1}{2}+\beta} \,\alpha \,+\, \mathsf{Ha} \,\mathsf{Coth} \,[\mathsf{Ha}]\right)}}\right},\,\,\left{\mathsf{GGS} \to \frac{1}{1+\sqrt{\left(\mathsf{A} \,\mathsf{C} \,\sqrt{\mathsf{Ha}} \,+\, \mathsf{A} \,\mathsf{Ha}^{\frac{1}{2}+\beta} \,\alpha \,+\, \mathsf{Ha} \,\mathsf{Coth} \,[\mathsf{Ha}]\right)}}\right}\right}$$

Select the positive GGS and plug into the efficiency expression:

$$ln[94]:= GGSmax := \frac{1}{1 + \sqrt{A c \sqrt{Ha} + A Ha^{\frac{1}{2} + \beta} \alpha + Ha Coth[Ha]}}$$

In[95]:= effpggs[Ha, A, GGSmax, Ψ] // FullSimplify

$$\left( e^{-\frac{\sqrt{\text{Ha}}}{A}} \left( \mathsf{A} + e^{\frac{\sqrt{\text{Ha}}}{A}} \left( -\mathsf{A} + \sqrt{\mathsf{Ha}} \right) \right) \sqrt{\sigma} \ \left( -1 + \mathsf{Ha} \, \mathsf{Coth} \, [\mathsf{Ha}] \right) \sqrt{\mathsf{A}} \, \mathsf{C} \, \sqrt{\mathsf{Ha}} + \mathsf{A} \, \mathsf{Ha}^{\frac{1}{2} + \beta} \, \alpha + \mathsf{Ha} \, \mathsf{Coth} \, [\mathsf{Ha}] \right) / \left( \sqrt{\mathsf{Ha}} \, \left( 1 + \sqrt{\mathsf{A}} \, \mathsf{C} \, \sqrt{\mathsf{Ha}} + \mathsf{A} \, \mathsf{Ha}^{\frac{1}{2} + \beta} \, \alpha + \mathsf{Ha} \, \mathsf{Coth} \, [\mathsf{Ha}] \right) \right) / \left( -1 + \mathsf{A} \, \mathsf{C} \, \sqrt{\mathsf{Ha}} + \mathsf{A} \, \mathsf{Ha}^{\frac{1}{2} + \beta} \, \alpha \right) \sqrt{\mu} \, \sqrt{\frac{\sigma}{\mu}} + \mathsf{Ha} \, \sqrt{\mu} \, \sqrt{\frac{\sigma}{\mu}} \, \mathsf{Coth} \, [\mathsf{Ha}] + \sqrt{\sigma} \, \left( 1 + \sqrt{\mathsf{A}} \, \mathsf{C} \, \sqrt{\mathsf{Ha}} + \mathsf{A} \, \mathsf{Ha}^{\frac{1}{2} + \beta} \, \alpha + \mathsf{Ha} \, \mathsf{Coth} \, [\mathsf{Ha}] \right) \right) \right)$$

In[168]:=

Effpmax[Ha_, A_] :=

$$\left(\left[1 + \frac{A\left(-1 + e^{-\frac{\sqrt{16a}}{A}}\right)}{\sqrt{\text{Ha}}}\right] \left(-1 + \text{Ha Coth [Ha]}\right)\right/ \left(1 + \sqrt{A \text{ c } \sqrt{\text{Ha}} + A \text{ Ha}^{\frac{1}{2} + \beta} \alpha + \text{Ha Coth [Ha]}}\right)^{2}$$

This is equivalent to:

$$\ln[49] := \left(1 + \frac{A\left(-1 + e^{-\frac{\sqrt{8a}}{A}}\right)}{\sqrt{\text{Ha}}}\right) \quad (-1 + \text{Ha Coth [Ha]}) \text{ GGSmax^2} // \text{FullSimplify}$$

$$\frac{\left[1 + \frac{A\left[-1 + e^{\frac{\sqrt{\alpha}}{A}}\right]}{\sqrt{\text{Ha}}}\right] \left(-1 + \text{Ha Coth [Ha]}\right)}{\left[1 + \sqrt{\frac{A\left(1. - e.956\,\text{Ha}\right)}{\sqrt{\text{Ha}}}} + \text{Ha Coth [Ha]}\right]^2}$$

In[104]:=

Plot[{Effpmax[Ha, 0.05], Effpmax[Ha, 0.1], Effpmax[Ha, 0.5], Effpmax[Ha, 0.6]}, {Ha, 1, 100}, PlotRange {0, 1}, PlotStyle "Rainbow", PlotTheme "ThickLines", LabelStyle {15, GrayLevel[0]}, PlotLegends {"A=0.05", "A-0.1", "A=0.5", "A=0.6"}]

Out[104]=

B2.3 Max Pump Efficiency at a given pressure head

We determine the electrode potential difference that maximizes the efficiency.

$$Simplify[(-Q[Ha, A, G, \Psi] *G) / (\Psi *DH[Ha, A, G, \Psi]), Assumptions \rightarrow \sigma > 0 && \mu > 0 && A > 0 && Ha > 0]$$

$$-\left[\left(2\,a^{2}\,e^{-\frac{\sqrt{6a}}{A}}\,G\left(A\left(-1+e^{\frac{\sqrt{6a}}{A}}\right)-e^{\frac{\sqrt{6a}}{A}}\,\sqrt{Ha}\right)\,\left(2\,a^{2}\,G+A\,Ha\,\sqrt{\mu\,\sigma}\,\,\Psi\right)\,\left(-1+Ha\,Coth\,[Ha]\,\right)\right]\right/$$

$$\left[A\,Ha^{3/2}\,\mu\,\Psi\left(-2\,a^{2}\,G\left(\sqrt{\frac{\sigma}{\mu}}\,-A\,\left(c+Ha^{\beta}\,\alpha\right)\,\sqrt{\frac{Ha\,\sigma}{\mu}}\right)-A\,Ha\,\sigma\,\Psi+2\,a^{2}\,G\,Ha\,\sqrt{\frac{\sigma}{\mu}}\,Coth\,[Ha]\,\right)\right]$$

Pumpef [Ha_, A_, G_, \Psi_] := $$-\left(\left(2\ a^2\ e^{-\frac{\sqrt{\mu_a}}{\Lambda}}\ G\left(A\left(-1+e^{\frac{\sqrt{\mu_a}}{\Lambda}}\right)-e^{\frac{\sqrt{\mu_a}}{\Lambda}}\ \sqrt{Ha}\right)\left(2\ a^2\ G+A\ Ha\ \sqrt{\mu\ \sigma}\ \Psi\right)\ (-1+Ha\ Coth\ [Ha])\right)\right/$$

$$\left(A\ Ha^{3/2}\ \mu\ \Psi\left(-2\ a^2\ G\left(\sqrt{\frac{\sigma}{\mu}}\ -A\left(c+Ha^\beta\ \alpha\right)\ \sqrt{\frac{Ha\ \sigma}{\mu}}\right)-A\ Ha\ \sigma\ \Psi+2\ a^2\ G\ Ha\ \sqrt{\frac{\sigma}{\mu}}\ Coth\ [Ha]\right)\right)$$

Ψs := -khh G/khe // FullSimplify

$$\Psi_s = \frac{-2 a^2 G}{A H a \sqrt{\sigma \mu}}$$

The pump efficiency can be simplified and written in terms of Ψ/Ψs , written as ΨΨs:

Pumpeff[Ha_, A_, $$\Psi\Psi$$ S_] := $\left(e^{-\frac{\sqrt{ha}}{\Lambda}}\left(A\left(-1+e^{\frac{\sqrt{ha}}{\Lambda}}\right)-e^{\frac{\sqrt{ha}}{\Lambda}}\sqrt{Ha}\right)\left(-1+\Psi\Psi$ S)\left(-1+HaCoth[Ha]\right)\right)\right(\Ha^{1/2}\PP\PS)\left(1-A\left(C+Ha^\beta\alpha\right)\right\right\Ha)

Plot[{Pumpeff[5, 0.05, $$\Psi\Psi$$ s], Pumpeff[5, 0.2, $\Psi\Psi$ s], Pumpeff[5, 0.5, $\Psi\Psi$ s], Pumpeff[20, 0.05, $\Psi\Psi$ s], Pumpeff[20, 0.05, $\Psi\Psi$ s], Pumpeff[500, 0.05, $\Psi\Psi$ s], Pumpeff[500, 0.05, $\Psi\Psi$ s], Pumpeff[500, 0.05, $\Psi\Psi$ s], PlotStyle $\rightarrow$ {Thick, Dashed, Dotted}, PlotTheme $\rightarrow$ "Classic", AxesLabel $\rightarrow$ {" $\Psi$ $\Psi$ $S$ $S$ $S$ $S$ $S$ $S$ $S$ $S$ $S$ $S$

Ψs is the electrode potential difference under stall. We determine the optimal Ψ/Ψs:

This is just the inverse of G/GS max, since G/Gs= (Ψ/Ψs) -1

In[146]:=

G / Gs

Out[146]:=

$$-\frac{G \, khh}{khe \, \Psi}$$

In[149]:=

$$(\Psi / \Psi s) ^ (-1)$$

Out[149]:

$$-\frac{G \, khh}{khe \, \Psi}$$

In[ ]:= Solve[D[Pumpeff[Ha, A, ΨΨs, α, β, c], ΨΨs] 0, ΨΨs] // FullSimplify

$$\left{\left{\Psi\Psi S\to\mathbf{1}-\sqrt{A\;c\;\sqrt{\text{Ha}}\;+A\;\text{Ha}^{\frac{1}{2}+\beta}\;\alpha\;+\;\text{Ha}\;\text{Coth}\left[\text{Ha}\right]}\;\right},\;\left{\Psi\Psi S\to\mathbf{1}+\sqrt{A\;c\;\sqrt{\text{Ha}}\;+A\;\text{Ha}^{\frac{1}{2}+\beta}\;\alpha\;+\;\text{Ha}\;\text{Coth}\left[\text{Ha}\right]}\;\right}\right}$$

Select the positive quantity:

ΨΨsm:=1 + $$\sqrt{\text{Ac}\sqrt{\text{Ha}}}$$

• A Ha $^{\frac{1}{2}+\beta}$ α + Ha Coth[Ha];
  Pumpeff[Ha, A, ΨΨsm] // FullSimplify

Out[156]=

$$-\left[\left(\mathbf{1.} \ e^{-\frac{\sqrt{\text{Ma}}}{A}} \left(-\mathbf{1.} \ A + e^{\frac{\sqrt{\text{Ma}}}{A}} \left(\mathbf{1.} \ A - \mathbf{1.} \ \sqrt{\text{Ha}}\right)\right) \ \left(-\mathbf{1.} + \text{Ha} \ \text{Coth} \ [\text{Ha}]\right) \ \sqrt{\frac{A \ (\mathbf{1.} - 0.956 \ \text{Ha})}{\sqrt{\text{Ha}}}} + \text{Ha} \ \text{Coth} \ [\text{Ha}]\right]\right] / \left(\left(\mathbf{1.} + \sqrt{\frac{A \ (\mathbf{1.} - 0.956 \ \text{Ha})}{\sqrt{\text{Ha}}}} + \text{Ha} \ \text{Coth} \ [\text{Ha}]\right)\right) \ \left(A \ (\mathbf{1.} - 0.956 \ \text{Ha}) + \sqrt{\text{Ha}} \left(\text{Ha} \ \text{Coth} \ [\text{Ha}] + \sqrt{\frac{A \ (\mathbf{1.} - 0.956 \ \text{Ha})}{\sqrt{\text{Ha}}}} + \text{Ha} \ \text{Coth} \ [\text{Ha}]}\right)\right)\right)\right)$$

This simplifies to:

$$\begin{aligned} & \text{Pumpeffm}\left[\text{Ha}{_},\,\text{A}{_}\right] \,:=\, \left(\text{e}^{-\frac{\sqrt{\text{Ha}}}{\Lambda}}\,\left(\text{A} + \text{e}^{\frac{\sqrt{\text{Ha}}}{\Lambda}}\,\left(-\text{A} + \sqrt{\text{Ha}}\,\right)\right)\,\left(-\text{1} + \text{Ha}\,\text{Coth}\left[\text{Ha}\right]\right)\,\right) / \ & \left(\sqrt{\text{Ha}}\,\left(\text{1} + \text{A}\,\text{c}\,\,\sqrt{\text{Ha}}\, + \text{A}\,\text{Ha}^{\frac{1}{2} + \beta}\,\alpha + \text{Ha}\,\text{Coth}\left[\text{Ha}\right]\, + 2\,\,\sqrt{\text{A}\,\text{c}\,\,\sqrt{\text{Ha}}\, + \text{A}\,\text{Ha}^{\frac{1}{2} + \beta}\,\alpha + \text{Ha}\,\text{Coth}\left[\text{Ha}\right]}\,\right)\right) \end{aligned}$$

Plot{Pumpeffm[Ha, 0.05], Pumpeffm[Ha, 0.1], Pumpeffm[Ha, 0.5], Pumpeffm[Ha, 0.6]}, {Ha, 1, 100}, PlotRange {0, 1}, PlotStyle "Rainbow", PlotTheme "ThickLines", PlotLegends {"A=0.05", "A=0.1", "A=0.5", "A=0.6"}, AxesLabel "Ha", "ϵ^{g Ψ},m", LabelStyle Directive[ Medium]

Comparing the obtained max efficiencies between the pump maximized via G/Gs or Ψ/Ψs:

In[170]:=

Clear $$[\alpha, \beta, c]$$

Pumpeffm[Ha, A] - Effpmax[Ha, A] // FullSimplify

0

B3. MHD sensor calculations

Measurement of electrical signals, such as the potential across the electrodes Ψ and/or the current Δ H inform on the flow rate (Q) and pressure drop in the conduit G. Consider, for example, open circuit conditions (Δ H =0). From the current equation we get G to be:

Gsensor := -KEE [A, Ha] * $$\psi$$ / KEH [A, Ha]

This simplifies to:

In[176]:=

Gsensor := $$\frac{\text{A Ha } \sqrt{\mu \, \sigma \, \psi}}{2 \, \text{a}^2 \, \left(-1 + \text{A C } \sqrt{\text{Ha}} + \text{A Ha}^{\frac{1}{2} + \beta} \, \alpha + \text{Ha Coth} \, [\text{Ha}]\right)}$$

In[ ]:= Q[Ha, A, Gsensor, ψ] // FullSimplify

$$\begin{split} &\left(2~a^2~e^{-\frac{\sqrt{\eta_a}}{A}}~\left(A+e^{\frac{\sqrt{\eta_a}}{A}}~\left(-A+\sqrt{Ha}~\right)~\right)~\psi~\left(-1+Ha~Coth~[Ha]~\right) \right.\ &\left.\left(\left(-1+A~c~\sqrt{Ha}~+A~Ha^{\frac{1}{2}+\beta}~\alpha\right)~\sqrt{\mu}~\sqrt{\sigma}~+\sqrt{\mu}~\sigma~+Ha~\sqrt{\mu}~\sqrt{\sigma}~Coth~[Ha]~\right)~\right)\right/\ &\left.\left(Ha^2~\mu~\left(A~\left(c+Ha^\beta~\alpha\right)~+\sqrt{Ha}~Coth~[Ha]~\right)~\left(-1+A~c~\sqrt{Ha}~+A~Ha^{\frac{1}{2}+\beta}~\alpha~+Ha~Coth~[Ha]~\right)~\right) \end{split}$$

This is equivalent to the following expression:

$$In[-]:=$$ KHEA, Ha $\psi$ // FullSimplify (*verify that answer is the same*)

$$2 \ a^2 \left( 1 + \frac{A \left( -1 + e^{\frac{\sqrt{m}}{A}} \right)}{\sqrt{\text{Ha}}} \right) \sqrt{\sigma} \ \psi \ \left( -1 + \text{Ha Coth [Ha]} \right) \ \left( 1 + \frac{\sqrt{\mu} \ \sqrt{\frac{\pi}{\mu}}}{\sqrt{\sigma} \left( -1 + \text{Ac } \sqrt{\text{Ha}} + \text{A Ha}^{\frac{1}{2}, \theta} \alpha + \text{Ha Coth [Ha]} \right)} \right) \ + \left( 1 + \frac{A \left( -1 + e^{\frac{\sqrt{m}}{A}} \right)}{\sqrt{\sigma} \left( -1 + \text{Ac } \sqrt{\text{Ha}} + \text{A Ha}^{\frac{1}{2}, \theta} \alpha + \text{Ha Coth [Ha]} \right)} \right)$$

This simplifies to:

In[178]:=

$$\operatorname{Qs}\left[\operatorname{Ha}{-},\,\operatorname{A}{-},\,\alpha{-},\,\beta{-},\,\operatorname{c}_{-}\right] := \frac{2\,\operatorname{a}^{2}\left(1+\frac{\operatorname{A}\left(-1+\operatorname{e}^{-\frac{\sqrt{\ln\alpha}}{\Delta}}\right)}{\sqrt{\operatorname{Ha}}}\right)\,\sqrt{\sigma}\,\,\psi\,\left(-1+\operatorname{Ha}\,\operatorname{Coth}\left[\operatorname{Ha}\right]\right)}{\operatorname{Ha}\,\,\sqrt{\mu}\,\left(-1+\operatorname{A}\,\operatorname{c}\,\,\sqrt{\operatorname{Ha}}\,+\operatorname{A}\,\operatorname{Ha}^{\frac{1}{2}+\beta}\,\alpha+\operatorname{Ha}\,\operatorname{Coth}\left[\operatorname{Ha}\right]\right)}$$

In[179]:=

Limit[Qs[Ha, A, 1, -1, -0.956], Ha $$\rightarrow$$ Infinity]

Out[179]=

0. if $$\left( \mathbf{a^2} \mid \frac{\mathbf{1}}{\sqrt{\mu}} \mid \sqrt{\sigma} \mid \psi \right) \in \mathbb{R} & A > 0$$

B3.1 MHD Sensor Applying the Onsager Casimir Reciprocity

If we use Onsager Relations instead: KEH=-KHE

$$ln[s]:=$$ KHEA, Ha $\psi$ // FullSimplify Out[s]=

$$\left( \text{2 a}^2 \; \text{e}^{-\frac{\sqrt{\text{Ma}}}{A}} \; \sqrt{\sigma} \; \psi \left( \text{A (-1 + Ha Coth [Ha])} \; + \; \text{e}^{\frac{\sqrt{\text{Ma}}}{A}} \; \left( \text{A} + \; \left( -\text{A Ha} + \text{Ha}^{3/2} \right) \; \text{Coth [Ha]} \right) \right) \right) \right) \left( \text{Ha}^2 \; \sqrt{\mu} \; \left( \text{A (c + Ha}^\beta \; \alpha) \; + \; \sqrt{\text{Ha}} \; \text{Coth [Ha]} \right) \right)$$

In[180]:=

ORQS [Ha_, A_, $$\alpha$$ , $\beta$ , $c_$ ] := $$\left(2 \text{ a}^2 \text{ e}^{-\frac{\sqrt{\mu}a}{\Lambda}} \sqrt{\sigma} \psi \left(\text{A (-1 + Ha Coth [Ha])} + \text{e}^{\frac{\sqrt{\mu}a}{\Lambda}} \left(\text{A + (-A Ha + Ha^{3/2}) Coth [Ha])}\right)\right) \right)$$

$$\left(\text{Ha}^2 \sqrt{\mu} \left(\text{A (c + Ha}^\beta \alpha) + \sqrt{\text{Ha} Coth [Ha]}\right)\right)$$

In[ ]:= Limit[ORQS[Ha, A, 1, -1, -0.956], Ha Infinity]

$$\boxed{ \textbf{0. if } \left( \textbf{a}^2 \mid \frac{\textbf{1}}{\sqrt{\mu}} \mid \sqrt{\sigma} \mid \psi \right) \in \mathbb{R} & \texttt{A} > \textbf{0} }$$

B4. MHD Brake

Here we define a brake where the fluid's motion is due to a pressure gradient is in the positive z-direction. We apply a negative potential such that a force in the negative z-direction is created to slow down the fluid. As long as Q ≥ 0 , we remain in the break region. Once the flow rate becomes negative, we are in a pump region once again.

$$H = K{EE} \psi + K{EH} G$$

$$Q = K{HE} \Psi + K{HH} G$$

When Q ≥ 0, we have:

$$\psi \le -\frac{K{H\,H}}{K{H\,E}}$$

In[181]:=

$$\psi$$ bmax[A_, Ha_, G_] := -KHH[A, Ha] / KHE[A, Ha] *G // FullSimplify

We can plug that into the H expression:

In[182]:=

Hbrake[A_, Ha_, G_] := (KEH[A, Ha] - KEE[A, Ha]*KHH[A, Ha]/KHE[A, Ha]) G // FullSimplify

From Physics we know that at Q=0 the forces cancel out such that:

$$\frac{H B_0}{2 a} + G = 0$$

$$\frac{H b r a k e}{G} = \frac{-2 a}{B_0}$$

In[*]:= Simplify[Hbrake[A, Ha, G] / G, Assumptions $$\rightarrow \mu > 0$$ ]

$$-\frac{2 a^2 \sqrt{\frac{\sigma}{\mu}}}{Ha}$$

Plugging in $$H_a = a B_0 \sqrt{\frac{\sigma}{\mu}}$$ , this becomes: $\frac{-2 a}{B_0}$ as desired