ELECTROHYDRODYNAMIC ACCELERATION OF CHARGING PROCESS IN A LATENT HEAT THERMAL ENERGY STORAGE MODULE

Description

CROSS-REFERENCES TO RELATED APPLICATIONS

This application claims the benefit of U.S. Provisional Patent Application No. 63/671,914, filed Jul. 16, 2024, the entire contents of which are hereby incorporated by reference for all purposes in its entirety.

BACKGROUND

Thermal energy storage using latent heat of phase change materials (PCMs) is a technique used in power plants, building ventilation, electronic cooling, and electric vehicle battery thermal management. For close to room temperature applications, engineers can use organic PCMs that can be eco-friendly in nature, low cost, and easily available. However, organic PCMs possess low thermal conductivity. Thus, the charging (melting) process in the energy storage systems using these PCMs is low. This can negatively affect the net energy storage density.

BRIEF SUMMARY

A latent heat thermal energy storage (LHTES) unit that employs unipolar charge injection-induced electrohydrodynamic (EHD) flow to accelerate the melting process is disclosed herein. A phase change material (PCM) based latent heat thermal energy storage (LHTES) unit can have high energy storage density, scalability, and near-constant temperature operation. The LHTES unit can be assisted with charge injection-induced electrohydrodynamic (EHD) flow to enhance the charging rate. The LHTES unit can enhance, compared with other units, the charging rate under different rates, whether weak, medium, or strong, of charge injection regimes. The EHD flow can intensify the flow velocity, can alter the flow structure, and can increase the heat transfer associated with the LHTES unit. The melting process can become more uniform and faster with the assistance of EHD flow. EHD flow at strong charge injection regime nullifies the effect of gravity and can lead to equal performance irrespective of the orientation of the LHTES unit. Shorter melting times and increased power storage capacity for the LHTES unit can be achieved in the presence of EHD flow.

In some embodiments, a latent heat thermal energy storage (LHTES) unit can include a shell, a tube, a phase change material (PCM), a first electrode, and a second electrode. The shell may be an outer housing, and the tube may extend through an interior of the shell. The PCM may be located between the shell and the tube. The first electrode may be located on the shell, and the second electrode may be located on the tube. The first electrode and the second electrode may be arranged to apply an electric potential from the tube to the shell and across the PCM, and the second electrode may be arranged to inject unipolar charge into the PCM to induce electrohydrodynamic (EHD) flow within the PCM.

In some embodiments, a first melting rate of the PCM exposed to the EHD flow may be higher than a second melting rate of the PCM not exposed to the EHD flow, and a first time to charge the LHTES unit may be lower than a second time to charge a different LHTES unit that does not include the EHD flow.

In some embodiments, the second electrode may be an emitter for unipolar charge injection, and the first electrode may be a collector electrode for unipolar charge injection.

In some embodiments, the second electrode may be configured to inject the unipolar charge into the PCM to cause approximately uniform melting of the PCM regardless of an orientation of the shell.

In some embodiments, the shell may be oriented vertically such that the tube extends parallel with respect to a direction of gravity.

In some embodiments, the shell may be oriented horizontally such that the tube extends perpendicularly with respect to a direction of gravity.

In some embodiments, the PCM may be an organic PCM that includes paraffin wax.

In some embodiments, the second electrode may be configured to adjust a charge injection into the PCM based on an intensity of the EHD flow and a performance of the LHTES unit.

In some embodiments, the EHD flow in the PCM may be configured to generate a uniform cylindrical velocity field between the shell and the tube.

In some embodiments, a method can be used for inducing melting in a latent heat thermal energy storage (LHTES) unit. The method can include applying an electric field between a first electrode of the LHTES unit and a second electrode of the LHTES unit in which the LHTES unit includes a phase change material (PCM). The method can include inducing electrohydrodynamic (EHD) flow in the LHTES unit with a unipolar charge injection. The method can include inducing, by the EHD flow, a phase change in the PCM.

In some embodiments, a method can be used for controlling a melting rate of phase change material (PCM) in a latent heat thermal energy storage (LHTES) unit. The method can include receiving input parameters comprising properties about the PCM, a target melting rate, and an orientation of the LHTES unit. The method can include determining an electric field and charge injection regime to induce EHD flow for achieving the target melting rate based on the properties about the PCM and the orientation. The method can include outputting control signals to a power supply to apply the determined electric field between a first electrode of the LHTES unit and a second electrode of the LHTES unit, wherein applying the electric field results in unipolar charge injection and EHD flow for melting the PCM.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a set of diagrams of an electro-hydrodynamics-assisted shell-and-tube latent heat thermal energy storage module, according to some embodiments.

FIG. 2 is a graph of time evolution of liquid fraction using three grids in different orientations, according to some embodiments.

FIG. 3 is an example of a 3D view of the mesh grid, according to some embodiments.

FIG. 4 is a set of plots of validation and results with respect to a numerical simulation involving the LHTES unit, according to some embodiments.

FIG. 5 is an example of a validation of numerical solver used in the simulation with respect to experimental results of unipolar charge injection induced electrohydrodynamic (EHD) flow assisted melting in a cubical cavity, according to some embodiments.

FIG. 6 is a set of plots that illustrated a growth of liquid fraction with time in an LHTES unit in different orientations, according to some embodiments.

FIG. 7 is a set of plots of kinetic energy density for the LHTES unit in a vertical orientation, according to some embodiments.

FIG. 8 is a set of plots of kinetic energy density for the LHTES unit in a horizontal orientation, according to some embodiments.

FIG. 9 is a set of plots of variation of mean temperature with respect to time for the LHTES unit, according to some embodiments.

FIG. 10 is a diagram of liquid fraction, velocity, and temperature fields for natural convection melting in a vertical orientation for the LHTES unit, according to some embodiments.

FIG. 11 is a diagram of liquid fraction, velocity, and temperature fields for natural convection melting in a horizontal orientation for the LHTES unit, according to some embodiments.

FIG. 12 is a diagram of liquid fraction, charge density, velocity, and temperature fields for applied voltage and charge injection in a vertical orientation for the LHTES unit, according to some embodiments.

FIG. 13 is a diagram of contours of velocity distribution in the vertical orientation of the LHTES unit, according to some embodiments.

FIG. 14 is a diagram of liquid fraction, charge density, velocity, and temperature fields for applied voltage and charge injection in a horizontal orientation for the LHTES unit, according to some embodiments.

FIG. 15 is a diagram of contours of velocity distribution in the horizontal orientation of the LHTES unit, according to some embodiments.

FIG. 16 is a set of plots of power stored for the LHTES unit at different charge injection regimes and at different applied voltages, according to some embodiments.

FIG. 17 is a set of plots of liquid fraction over time for the LHTES unit in different orientations, according to some embodiments.

FIG. 18 is a set of plots of maximum liquid velocity over time for the LHTES unit in different orientations, according to some embodiments.

FIG. 19 is a set of plots of mean Nusselt number over time for the LHTES unit in different orientations, according to some embodiments.

DETAILED DESCRIPTION

Certain aspects and features involve a three-dimensional shell-and-tube latent heat thermal energy storage (LHTES) module incorporated with a unipolar charge injection mechanism. The free charges can be influenced by the Coulomb force of the applied electric field and can move towards the oppositely charged electrodes. The moving charges can transfer momentum to the liquid phase change material (PCM) and can induce an electrohydrodynamic (EHD) fluid motion. The vortices of the EHD flow can increase fluid mixing and convectional heat transfer, which can result in faster melting rates compared to conventional LHTES units without any acceleration mechanisms.

Energy storage units have become an integral part of energy systems based on renewable sources, recovery of waste heat, building cooling and ventilation, battery thermal management and electronics. High volumetric efficiency, mechanical and chemical stability, and fatigue resistance are associated with LHTES units working based on melting and solidification of PCMs. Eco-friendly, organic PCMs suffer from low thermal conductivity and, hence, the energy storage rates and thermal response speed of LHTES units are limited.

Energy storage rates and thermal response speeds can be enhanced via passive methods and/or active methods. Passive methods may be simple in design and/or execution and may include (i) increasing the surface area for heat transfer and (ii) inclusion of additives to increase the thermal conductivity. Using extended fin structures and encapsulated PCMs can increase the surface area for heat transfer. Application of metallic fins can increase the net area available for heat transfer and can be effective during the solidification (discharging process) of PCM in LHTES units. Micro-encapsulation of PCMs makes each capsule behave as an individual LHTES unit. Each micro-capsule has a high surface-to-volume ratio and can exhibit high melting and solidification rates. Another passive method is to use metallic nano-additives or embedding metallic foams to enhance the effective thermal conductivity of PCMs. In addition to increasing the thermal conductivity, including nanoparticles can improve the thermal stability of PCMs and energy storage capacity for a given volume of PCM. Additionally or alternatively, nano-additives serve as nucleation points for the onset of melting. Metallic foams can aid in improving the structural stability of the PCM.

Active methods can use secondary fluid flow generated by an external mechanism to intensify the convective heat transfer. Acoustic and/or mechanical vibration-induced flows can be used to speed up the melting rate in LHTES units. The flow induced by the vibrations can intensify the flow in the liquid region, which can lead to higher melting rates. However, incorporating devices to generate acoustic or mechanical vibrations can increase the complexity of the systems. LHTES units fitted with vibration generating devices may possess poor mechanical stability and high total power consumption. Magnetic fields may be used to control the melting process, and the magnetic field's orientation may relate to a speed of the melting process. The magnetic field may induce swirl motion in the melted PCM zone and can improve the thermal response of the LHTES unit.

Electrohydrodynamic (EHD) flows generated by applying an electric field can be used as an active heat transfer enhancement method. EHD-based active heat transfer enhancement possesses several positive characteristics such as simplicity, short response time, and operation without mechanical vibration and noise. EHD flows can be employed in single or multi-phase (liquid-gas) systems to enhance heat transfer. Additionally or alternatively, the electric forces can pull semi-solid PCM in the interface region into the liquid bulk (solid extraction phenomenon). The formation of electro-convective flow cells can improve the heat transfer and melting rates. In some examples, electric field strength, charge generation mechanism, electric field orientation, and polarity of applied electric potential can be factors influencing the EHD-assisted melting process.

Eco-friendly organic PCMs exhibit good dielectric properties and possess high dielectric strength. Thus, organic PCMs are suitable for EHD applications, and the additional electric power consumption is low. EHD-based enhancement of melting of organic PCMs can be an effective and economical option. In 2-D geometry, the EHD flow vortices can increase the mixing and thermal transport in the liquid PCM, leading to shorter melting times. The EHD vortices can nullify the effect of gravity and can lead to the same melting performance at both vertical and horizontal orientations.

The 3D flow structures of the electro-convective flow can influence flow and melt interface morphology. The charging performance of a 3D LHTES unit in the presence of unipolar charge injection can be determined. A 3D numerical investigation of the charging process in a shell-and-tube LHTES module exposed to unipolar charge injection under three regimes, two applied voltages, and two orientations can be performed. The flow, heat transfer, and melting characteristics affected by the charged injection induced EHD flow can be determined. The 3D structures of the flow, thermal, liquid fraction, and charge density can be visualized. The total time taken for charging and net power storage capacity of the LHTES unit at different charge injection strengths, applied voltages, and orientations can be determined.

FIG. 1 is a set of diagrams of an electro-hydrodynamics-assisted, shell-and-tube latent heat thermal energy storage (LHTES) module, according to some embodiments. FIG. 1 illustrates (a) a graphical representation 102 of the electro-hydrodynamics-assisted, shell-and-tube LHTES module 100 and (b) a computational domain representation 104 of the electro-hydrodynamics-assisted, shell-and-tube LHTES module 100. In some embodiments, FIG. 1 illustrates a 3D schematic sketch of the shell-and-tube LHTES module and the computational domain thereof. As illustrated in FIG. 1, the electro-hydrodynamics-assisted, shell-and-tube LHTES module 100 can have height H=150 mm, and the shell and tube can have diameters DS=78 mm and DT=24.6 mm, respectively, though other suitable measurements for each of the foregoing are possible. The geometrical dimensions of the electro-hydrodynamics-assisted, shell-and-tube LHTES module 100 (herein after the LHTES unit 100) can be chosen based on historical data and to facilitate visualization of flow structures.

At the beginning of a simulation (time t=0 s), solid paraffin wax at approximately θ_∞=298.15 K can occupy the entire flow domain. Hot fluid at approximately 361.15 K can flow through the tube 106 of the LHTES unit 100. Given the small height of the LHTES unit 100 and the sufficiently high flow rate, the temperature difference between the inlet 108 and outlet 110 of the tube 106 may be negligible. Thus, the temperature boundary condition along the surface of the tube 105 can be a constant temperature of approximately θ_H=361.15 K. The walls are considered to have no-slip and no-penetration boundary conditions. Except for the wall of the tube 106, the walls can be thermally insulated. High electric potential (VH) can be applied over the wall of the tube 106, which acts as the emitter electrode. The wall of the shell 112 can be the grounded electrode (V=0 kV). The top surface 114 and the bottom surface 116 of the LHTES unit may be electrically insulated. For example, positive charges can be injected from the surface of the tube 106 (emitter electrode) radially outwards to the collector electrode (shell wall).

TABLE 1
List boundary conditions used in the simulations

Location	Velocity	Temperature	Voltage	Charge Density
Tube	{right arrow over (u)} = 0	θ = 361.15 K	V = VH	q = q0
Shell	{right arrow over (u)} = 0	∂ θ ∂ n = 0	V = 0 kV	∂ q ∂ n = 0
Top Wall	{right arrow over (u)} = 0	∂ θ ∂ n = 0	∂ V ∂ n = 0	∂ V ∂ n = 0
Bottom Wall	{right arrow over (u)} = 0	∂ θ ∂ n = 0	∂ V ∂ n = 0	∂ V ∂ n = 0

The list of boundary conditions used in the computational simulations is given in Table 1 produced above. A constant, uniform, and autonomous charge injection is assumed at the surface of the tube 106. Accurate experimental measurement of charge injection strength from the high voltage electrode surface may be difficult to measure. Thus, a dimensionless parameter C is used to quantify the charge injection strength from the high voltage electrode. The expression for C is as follows:

C = q 0 ⁢ ( D s - D T ) 2 ⁢ ε l ⁢ V H ( 1 )

Here, q₀ represents the volumetric charge density distribution emitted from the high-voltage electrode, and ϵ₁ represents the permittivity of the material. In experimental conditions, the charge injection strength can vary with time due to the deterioration of the electrode and dielectric properties of the working fluid. The charge injection parameter at weak (C=0.1) medium (C=1), and strong (C=10) regimes can be considered to evaluate the effect of charge injection strength. The charging performance of the LHTES unit 100 can be determined under all three charge injection regimes at applied voltages VH=0 kV, 5 kV, and 10 kV and in vertical orientation and horizontal orientation.

The combined set of governing equations to simulate the solid-liquid phase change assisted with an EHD flow generated by unipolar charge injection can be formulated based on various assumptions. For example, the flow may be a three-dimensional (3D), incompressible, Newtonian flow. Additionally or alternatively, the phase change material (paraffin wax) may be or include a dielectric material with poor electrical conductivity, so the electric current, magnetic field effects, and heating effects may be neglected. Additionally or alternatively, material properties may be independent of temperature, but the material properties can be considered different for the solid and liquid phases of the PCM. Additionally or alternatively, the density variation can be considered in the momentum equation according to the Boussinesq approximation. Additionally or alternatively, the PCM can be rigidly fixed to the shell wall. Additionally or alternatively, the solid thickness of the walls of the shell 112 and the tube 106 can be neglected for simplicity of the numerical model. The following equations of mass conservation (2), momentum conservation (3), energy conservation (4), electric potential (5), electric field (6), and charge conservation (7) can be or form the governing equations.

∇ · u → = 0 ( 2 ) ∂ u → ∂ t + ∇ ( u → ⁢ u → ) = - ∇ p ^ + v l ⁢ ∇ 2 u → - g → [ 1 - β ⁡ ( θ - θ ∞ ) ] + q ⁢ E → + 1 2 ⁢ E → 2 ⁢ ∇ ε + S u ( 3 ) ∂ θ ∂ t + ∇ · ( u → ⁢ θ ) = ∇ · ( α l / s ⁢ Δ ⁢ θ ) + S θ ( 4 ) ∇ · ( ε l / s ⁢ ∇ V ) = - q l / s ( 5 ) E → l / s = - ∇ V ( 6 ) ∂ q l / s ∂ t + ∇ · J → = 0 ; J → = ( K l / s ⁢ E → l / s + u → ) ⁢ q l / s - D l / s ⁢ ∇ q l / s ( 7 )

In Equations (2)-(4) and (7) {right arrow over (u)} and {right arrow over (E)} are the vectors of flow velocity and electric fields, respectively. {right arrow over (g)}[1−β(θ−θ_∞)] can represent the volumetric body force term acting on the fluid due to gravity, calculated as per the Boussinesq approximation. β denotes the thermal expansion coefficient of the PCM, and {right arrow over (g)} is the gravitational acceleration. The terms q{right arrow over (E)} and ½{right arrow over (E)}²∇ϵ in Equation (3) are the electric force terms due Coulomb and dielectric forces, respectively. In the energy equation, ν, θ and α are the kinematic viscosity, temperature, and thermal diffusivity of the PCM. The term S_θ in the energy equation is calculated as Equation 8.

S θ = - [ ∇ · ( u → · ( Lf ) + ∂ ( Lf ) ∂ t ) ] C p . ( 8 )

In Equation (8), L denotes the enthalpy of fusion, C_p represents the specific heat capacity, and f denotes the local liquid content (liquid fraction) in a computational cell. The value of f in each cell can be estimated based on temperature θ in Equation 9.

f ⁡ ( θ ) = ⁢ { 0 θ ≥ θ liq θ - θ sol θ liq - θ sol θ liq > θ ≥ θ sol 1 θ ≥ θ sol . ( 9 )

In Equation 9, θ_liq and θ_sol are the upper and lower bounds of the melting temperature range, respectively. In the momentum equation, S_u=A{right arrow over (u)} in Equation (3) can be added as a source term used to turn off the velocity in the computational cells identified to be in solid state. A is linear function f expressed as A=−ζ[1−f]. The value of ζ is set very high to nullify the contribution of all other terms in Equation (3), in the computational cells identified to be in solid or semi-solid state. In the cells which are in liquid state f=1, the value of S_u is zero, and the momentum equation without S_u is solved to obtain the velocity distribution. Equations (5) and (7) are the Poisson equation for electric potential and the Nernst-Planck equation for charge conservation in which V, ϵ, {right arrow over (J)}, q, K and D denote the electric potential, dielectric permittivity of the material, current density, charge density, ionic mobility and charge diffusivity of the PCM. The subscripts ‘l/s’ denote the properties of PCM in liquid or solid states. The properties of the phase change material (paraffin wax) are listed in Table 2 produced below. Walden's rule is used to calculate the liquid ionic mobility of PCM, and its value in the solid state is considered to be approximately 10 times its value in liquid state.

TABLE 2
List of properties of PCM (paraffin wax)

Property	Velocity

Liquid Density [ρl]	760.00	kgm−3
Solid Density [ρs]	870.00	kgm−3
Liquid Specific Heat Capacity [Cp, l]	2400.00	Jkg−1 K−1
Solid Specific Heat Capacity [Cp, s]	1800.00	Jkg−1 K−1
Liquid Thermal Conductivity [kl]	0.15	Wm−1 K−1
Solid Thermal Conductivity [ks]	0.24	Wm−1 K−1
Dynamic Viscosity [μl]	3.42 × 10−3	kgm−1s−1

Thermal Expansion Coefficient [β]

1 × 10−6−K−1

Latent Heat of Fusion [L]

192.00

kJkg−1

Solidus Temperature [θsol]	331.00K
Liquidus Temperature [θliq]	333.00K

Liquid Dielectric Permittivity [εl]	1.771 × 10−11	Fm−1
Solid Dielectric Permittivity [εs]	1.992 × 10−11	Fm−1
Liquid Ionic Mobility [Kl]	7.50 × 10−8	m2s−1V−1
Solid Ionic Mobility [Ks]	7.50 × 10−7	m2s−1V−1

Kinetic energy density can be a measure for the fluid motion intensity in the computational domain. The volume averaged kinetic energy density is calculated using Equation 10.

KE = 1 v f ⁢ ∫ v f u → 2 ⁢ dv ( 10 )

ν_f indicates the liquid volume in the computational domain at any given time moment. A higher value of KE may indicate stronger the fluid motion, and vice versa. The charging enhancement ratio can be the ratio of charging time taken by pure natural convection melting to the charging time in the presence of electric field. The time taken to reach 80% melting can be a reference. Thus, the expression for charging enhancement ratio is given as Equation 11.

Charging ⁢ enhancement ⁢ ratio = Time ⁢ taken ⁢ to ⁢ reach ⁢ 80 ⁢ % ⁢ melting ⁢ without ⁢ electric ⁢ field Time ⁢ taken ⁢ to ⁢ reach ⁢ 80 ⁢ % ⁢ melting ⁢ with ⁢ electric ⁢ field . ( 11 )

The total power stored in the 3D shell-and-tube LHTES device is calculated in Equation 12.

Net ⁢ power ⁢ stored ⁢ P = [ m l · ( C p , l · ( θ l - θ ∞ ) + L ) ] + [ m s · C p , s · ( θ s - θ ∞ ) ] Time ⁢ required ⁢ to ⁢ reach ⁢ 80 ⁢ % ⁢ melting . ( 12 )

Here, m_l and m_s are the liquid and solid mass at 80% melting. θ_l and θ_s are the volume averaged mean temperatures for the liquid and solid phases at 80% melting. The first term in the numerator denotes the net heat stored in the liquid phase of the PCM in which m_l·C_p,l[θ_l−θ_∞] is the sensible heat content and m_l. L is the latent heat content. The term [m_s·C_p,s·(θ_s−θ_∞)] denotes the sensible heat contribution of solid PCM. The extra electric power consumed to initiate charge injection and the resultant EHD flow can be determined using Equation 13.

P EL = V × I tot , ( 13 )

V is the applied voltage and I_tot is the total electric current in the system. The instantaneous total electric current at a given time can be calculated using Equation 14.

I tot = ∫ ∫ ( J → + ∂ E → ∂ t ) · n ^ ⁢ dS . ( 14 )

The governing equations for the solid-liquid phase change assisted with an EHD flow due to unipolar charge injection can be implemented in the open-source finite-volume method (FVM) framework of OpenFOAM. The governing equations can be discretized using the standard finite-volume procedures in OpenFOAM. The least squares method may be used to discretize the Laplacian terms in the governing equations. A second-order linear scheme can be employed to discretize the gradient terms. The advective terms in the momentum and energy equations can be discretized using central-differencing scheme. The term K{right arrow over (E)}q in Equation (7) can be the electromigration term. This term is transformed to the form of a traditional convective term to enable implicit discretization. Charge diffusion in dielectric materials may be negligible. The ionic diffusion coefficient D can have a very small value such as D=1×10⁻⁸. The negligible diffusion contribution in the charge transport equation can lead to strong convection dominance and may cause Equation (7) experience numerical instability.

To handle numerical instability, a total variation diminishing (TVD) scheme, along with a deferred correction approach, can be used for the advective derivatives in the charge conservation equation. The momentum source term S_u can be an explicit volumetric body force term and can be discretized using the standard FVM procedures. The enthalpy-porosity source term in Equation (4) may involve special treatment. The latent heat content in each control volume can be estimated as per the local temperature calculated in each control volume for every iteration. Then, the liquid fraction in each computational cell can be updated based on the local heat content, which differentiates the cells into solid and liquid states. The Euler backward method can be used to discretize the time derivatives. Residual convergence criteria of 1×10⁻⁸ is set for the variables in each time step.

FIG. 2 is a graph 200 of time evolution of liquid fraction using three grids in different

orientations, according to some embodiments. As illustrated, FIG. 2 includes (a) a plot 202 of a vertical orientation time evolution of liquid fraction using three grids and (b) a plot 204 of a horizontal orientation time evolution of liquid fraction using three grids. On the plot 202 and the plot 204, a horizontal axis 206 represents time, and a vertical axis 208 represents a fraction of liquid in the LHTES unit 100.

A 3D, uniform, structured mesh can be considered in the simulation. A grid sensitivity analysis can be performed to choose the appropriate grid. For example, three different grids, A (700×10³ cells), B (1.4×10⁶ cells), and C (2.1×10⁶ cells), can be considered in both orientations, horizontal and vertical. The growth of the total liquid fraction in the domain can be compared for the three grids in both orientations for cases with V=0 kV and 10 kV such as in FIG. 2. The variations between grids B and C may be less than 3% for the considered parameters. Thus, grid B with 1.4×10⁶ cells can be used, as in FIG. 3 that illustrates an example of a 3D view of the a mesh grid 300, according to some embodiments. Adaptive time stepping can be used in the simulations by setting maximum limits for Courant number and time step size as 0.5 and 0.001 s, respectively. Simulations for melting without any applied electric field may produce satisfactory results with much larger time steps (0.1 s). But, in the cases with electric field, the hyperbolic Nernst-Planck equation for charge transport may diverge with time step size 0.1 s. Furthermore, smaller time steps may be used to ensure the time independence of unsteadiness in the flow variables in the melting process with electric field. Thus, a time step size of 0.001 can be used for the simulations.

The physical problem under consideration deals with the charging process in a three-dimensional shell-and-tube LHTES unit assisted with unipolar charge injection. The present scenario is a combination of transport of electric charges an interface (solid-liquid), and melting process assisted with an EHD flow. The numerical solver can be validated to simulate pure natural convection melting (no electric field) in the LHTES unit 100.

FIG. 4 is a set of plots 400 of validation and results with respect to a numerical simulation involving the LHTES unit 100, according to some embodiments. As illustrated, FIG. 4 includes (a) a plot 402 of an example of validation of numerical simulations with respect to experimental data for time evolution of temperature at two spatial locations A and B in a vertical LHTES unit, and (b) a plot 404 of an example of experimental and numerical results of melt radius variation in horizontal direction of LHTES unit at two-time instances. In the plot 402, a horizontal axis 406 can represent time, and a vertical axis 408 can represent temperature. In the plot 404, a horizontal axis 410 can represent distance, and a vertical axis 412 can represent a melt radius.

A 2D numerical simulation can be performed considering dimensions and boundary conditions. The comparison of temperature profiles between the experimental data and present numerical simulations can be presented in the plots of FIG. 4 in (a). The root mean square (RMS) deviation of the numerical results as compared to experimental measurements at points A and B may be 2.22 K and 3.86 K. The comparison of present simulation data with experimental results of the interface radii at two-time moments (11 and 29 mins) can be presented in FIG. 4 (b). The maximum deviation between the experimental and present numerical simulation is 7.63% noted at 29 minutes and distance 0.2 m. Whereas, the mean deviation between the present simulations and the experimental data is only 2.82%.

EHD-assisted melting of n-octadecane in a cubical cavity with a cylindrical wire electrode placed at the center of cavity can be considered. 3D simulations can be set up. Two cases A and B are considered, with the case A having only the bottom wall as electrically grounded, while case B considers the four side walls (left, right, bottom, and top) to be grounded. The applied voltage in the cylindrical wire electrode is −25 kV. The occurrence of charge injection at electric field strengths in the order of 106 V/m can be determined. The charge density boundary condition at the cylindrical electrode surface is calculated based on the injection parameter C. The simulations can be run with three values for C such as 0.1, 1, and 10. The numerical results with C=1 agree well with the experimental data for time evolution of total liquid fraction such as in FIG. 5.

FIG. 5 is an example of a validation of numerical solver used in the simulation with respect to experimental results of unipolar charge injection induced electrohydrodynamic (EHD) flow assisted melting in a cubical cavity, according to some embodiments. As illustrated, FIG. 5 includes (a) a sectional view of the LHTES unit 100 and (b) a plot 502 of unipolar charge injection-induced EHD flow-assisted melting. As illustrated in FIG. 5, the LHTES unit 100 may be cubical, though other shapes are possible for the LHTES unit 100. For example, the LHTES unit 100 can include multiple adiabatic walls, such as adiabatic wall 504, and can include a high-voltage hot cylinder 506. The plot 502 can include a horizontal axis 508 that represents time, and the plot 502 can include a vertical axis 510 that can represent a liquid fraction. In some examples, the maximum deviation may be approximately 4.16%. The validation cases represented by the plot 502 may confirm the ability of the customized OpenFOAM solver to simulate natural convection melting aided with an EHD flow due to unipolar charge injection. The electric field strength in the 3D shell-and-tube LHTES setup may be in the order of 10^{6 V}/m.

The numerical simulations can be performed up to a total time t=5000 s. The roles of charge injection strength, applied electric potential and orientation of the LHTES unit on the charging process can be analyzed. The charge injection strength may be varied in weak, medium, and strong injection regimes with the injection parameter C=0.1, 1, and 10, respectively. Two applied voltages V=5 and 10 kV can be considered at each injection strength. At each combination of charge injection strength and applied voltage, the melting performance is evaluated at different orientations (Case 1—Horizontal and Case 2—Vertical). The charging performance of the LHTES unit 100 assisted with EHD flow is compared with the base case of natural convection melting (no electric field V=0 kV).

FIG. 6 is a set of plots 600 that illustrated a growth of liquid fraction with time in an LHTES unit 100 in different orientations. The EHD-assisted melting mechanism can be analyzed by monitoring the progress of liquid fraction, mean kinetic energy density, and mean temperature. The growth of liquid fraction with time in vertical (case 1) and horizontal (case 2) orientations at V=10 kV and different injection regimes (weak, medium and strong) are presented in (a) and (b) of FIG. 6, respectively. For example, FIG. 6 illustrates a plot 602 corresponding with the vertical case and a plot 604 corresponding with the horizontal case. The plot 602 can have a horizontal axis 606 representing time and a vertical axis 608 representing liquid fraction for the vertical case, and the plot 604 can have a horizontal axis 610 representing time and a vertical axis 612 representing liquid fraction for the horizontal case.

The melting rate curve for the natural convection melting process (solid black line) is provided for comparison in the set of plots 600. Initially, solid PCM can occupy the entire LHTES unit 100. The solid PCM begins to melt from the hot tube surface. The liquid fraction in the flow domain gradually increases with time, and for both cases 1 and 2 (at V=0 kV), the maximum liquid fraction achieved at t=5000 s is around 80%. In case 1, the slope of the melting curve steadily increases with the progress of time until it reaches 80.1% at t=5000 s. In case 2, the slope of the melting rate curve is higher in the early melting stages until t=2500 s. Nevertheless, at t>2500 s, the melting rate begins to decrease, and it reaches 79.9% at t=5000 s. In vertical orientation (case 1), the characteristic length L that determines the strength of natural convection is the entire height of the LHTES unit (H=150 mm). The Rayleigh number

R ⁢ a = g ⁢ β ⁢ θ ⁢ L 3 v ∝ ,

which is a non-dimensional number to characterize the strength of natural convection, is around 4×10⁹.

The characteristic length that determines the natural convection in horizontal orientation (case 2) may be the effective radius of the PCM cavity ((D_S−D_T)/2=26.7 mm). The resultant Rayleigh number may be around 2.3×10⁷. Additionally or alternatively, the natural convection in case 2 can take place in the region above the heated tube of the LHTES unit 100. The heat transfer in the region below the heated tube in horizontal orientation may primarily be governed by conduction. Natural convection in case I can take place along the entire height of the heated tube. In some examples, the role of natural convection in case 2 may be weaker, but still present, compared to case 1. For cases with electric field V=10 kV, the liquid fraction curve matches with the case without electric field for a short period in the preliminary stages of melting. During this period, the liquid volume in the computational domain may be less and may be closer to the tube wall. The low liquid volume and the proximity to the tube wall may restrict fluid motion. The melting process during this short period may be primarily due to conduction heat transfer.

At around t=100 s, the melting rate curves in the set of plots 600 for the cases with the electric field may deviate from those without the electric field. The point of deviation marks the onset of electroconvection due to the electric field. The onset of electroconvection leads to faster melting rates, as indicated by the increase in slopes of the melting curves. The melting rate in the presence of the electric field may be higher than in the case without an electric field. At weak injection regime (C=0.1) in case 1, 100% melting is reached at t=5000 s. In case 2, 93% melting at t=5000 s for C=0.1 is achieved. At medium and strong injection regimes (C=1 and 10), the melting rate curves are approximately identical irrespective of the orientation. This indicates that the melting process may be governed by electroconvection in medium injection and strong injection regimes, and the influence of buoyancy force is nullified. The melting rate in the strong injection regime may be marginally higher than in the medium injection regime. 100% melting is achieved at C=1 and 10 in both cases 1 and 2. The melting process increases in the presence of unipolar charge injection induced EHD flow.

FIG. 7 is a set of plots 700 of kinetic energy density for the LHTES unit 100 in a vertical orientation, and FIG. 8 is a set of plots 800 of kinetic energy density for the LHTES unit 100 in a horizontal orientation. As illustrated in FIG. 7, the set of plots 700 includes a plot 702 corresponding with no applied voltage and a plot 704 corresponding with applied voltage. As illustrated in FIG. 8, the set of plots 800 includes a plot 802 corresponding with no applied voltage and a plot 804 corresponding with applied voltage. The plot 702 includes a horizontal axis 706 that represents time and a vertical axis 708 that represents kinetic energy density, and the plot 704 includes a horizontal axis 710 that represents time and a vertical axis 712 that represents kinetic energy density. Additionally or alternatively, the plot 802 includes a horizontal axis 806 that represents time and a vertical axis 808 that represents kinetic energy density, and the plot 804 includes a horizontal axis 810 that represents time and a vertical axis 812 that represents kinetic energy density.

To further understand the effects of the electric field on the melting process, the time evolution of kinetic energy density KE in the computational domain, for cases 1 and 2, are plotted in FIGS. 7 and 8, respectively. The net fluid motion intensity in the computational domain can be quantified in terms of kinetic energy density. The kinetic energy density in the LHTES unit 100 at a given time moment can be determined in Equation 10. The time evolution plots of kinetic energy density in cases 1 and 2 and at C=0.1, 1, and 10 for applied voltage V=0 kV and 10 kV are illustrated in FIGS. 7 and 8. For each of the simulations, the kinetic energy density curve is close to zero initially, which may indicate that conduction dominates the melting process. The point where the kinetic energy density curves exhibit non-zero values may indicate the onset of fluid motion. For cases 1 and 2 (at V=0 kV, in the plot 702 and the plot 802), the kinetic energy density curves rise sharply after the onset of fluid motion and then fall down. During the initial melting process, the smaller liquid volume and the higher temperature gradient leads to stronger buoyancy induced fluid motion. Thus, a sharp rise in the kinetic energy density is present. However, with further progress of melting, the liquid volume in the domain increases, the temperature gradient lowers, and the net buoyancy force gets weaker.

This can lead to the fall of the kinetic energy density in the domain at later stages. The initial peak seen in case 2 is higher than case 1. Likewise, the fall in kinetic energy density in case 2 is also steeper than in case 1. This can be related to the variation of slope of the melting curve in case 2 for V=0 kV (refer to FIG. 6). This confirms that the natural convection may be stronger initially in the horizontal orientation, but the later stages of melting are marked with weak natural convection. For the cases with electric field (V=10 kV), the kinetic energy density curves are higher for higher charge injection strengths in both cases 1 and 2. Higher charge injection strength may lead to stronger Coulomb force. This can be understood by the expressions for charge injection parameter C and Coulomb force in Equations 1-3. From Equation 3, the charge density q₀ emitted from the electrode surface may be higher at higher values of C. Larger values of q can yield stronger Coulomb forces. Stronger injection regimes can result in higher mean kinetic energy density and can indicate stronger EHD flow at higher charge injection regimes. The kinetic energy density curves may be approximately identical and of the same magnitude for cases 1 and 2 at any given charge injection strength.

The kinetic energy density may not fall in the later stages of melting. The electro-convective fluid motion may be independent of the buoyancy force and may be strong enough to sustain the fluid motion throughout the entire time period of 5000 s. The kinetic energy density in the cases with an electric field may be approximately two to three orders of magnitude higher than in the case without the electric field. For instance, the maximum value of kinetic energy density in case 1 without electric field is around 450 μJ/m³, whereas the maximum value of kinetic energy density in case 1 with electric field is around 0.4 J/m³ (refer to FIG. 7). Additionally or alternatively, the maximum values of kinetic energy density in case 2 without and with electric field are around 630 μJ/m³ and 0.4 J/m³ (refer to FIG. 8), respectively. Despite the higher kinetic energy density in C=10 compared to C=1, the growth of the melting curve is approximately identical for medium injection and strong injection regimes.

Higher electric force in the strong injection regime can lead to stronger fluid motion than observed in the medium injection regime. The net melting rate may be restricted by the gap between the tube and shell. The stronger fluid motion at C=10 (compared to C=1) may result in a very marginal increase in melting rate. Both medium injection and strong injection regimes can exhibit much faster melting as compared to C=0.1 and pure natural convection melting. The oscillations present in the kinetic energy density curves of C=1 and C=10 can represent the chaotic periodic motion of the electro-convective flow.

FIG. 9 is a set of plots 900 of variation of mean temperature with respect to time for the LHTES unit 100. As illustrated in FIG. 9, the set of plots 900 includes a plot 902 corresponding with no applied voltage and a plot 904 corresponding with applied voltage. The plot 902 can include a horizontal axis 906 representing time and can include a vertical axis 908 representing a mean temperature with no applied voltage. The plot 904 can include a horizontal axis 910 representing time and can include a vertical axis 912 representing a mean temperature with applied voltage.

The variation of the mean temperature with respect to time in the computational domain for cases 1 and 2 at V=0 and V=10 kV are illustrated in (a) and (b) of in FIG. 9, respectively. Mean temperature in the domain may signify the heat transfer from the hot tube wall. A higher mean temperature in the domain may indicate a higher rate of heat transfer from the hot tube wall to the PCM. The initial stages of melting of the simulations have approximately the same mean temperature in the domain due to the dominance of conduction heat transfer. For the cases without electric field, the mean temperature smoothly increases with time. The mean temperature curve of case 2 without electric field shows a deflection point along the curve at around t=1500 s. This observation can be explained with respect to the variations in the flow intensity. In plot 702 of FIG. 7 (case 1 without electric field), at around time t=1500 s, the mean kinetic energy density in the domain begins to fall gradually with respect to time. The mean temperature curve of case 1 without electric field (sec plot 902 of FIG. 9) also increases smoothly without any noticeable deflection point. The kinetic energy density curve for case 2 without electric field drops close to t=1500 s (refer to plot 802 of FIG. 8). The corresponding mean temperature curve in the plot 904 of FIG. 9 may also exhibits a deflection point at around t=1500 s. After the onset of fluid motion, the mean temperature in the domain for the cases with the electric field is higher than the case of pure thermal convection melting.

Stronger charge injection may lead to stronger fluid motion and higher heat transfer. The difference between the medium injection and strong injection regimes may be negligible. The heat transfer may be limited by the gap between the shell and tube, irrespective of the stronger fluid motion in strong injection regime. A sudden rise in mean temperature curves of medium injection and strong injection regimes at time close to 2250 s may be present. At this point, the melt fraction is close to 97%, and the flow features may become similar to single-phase electro-convection in an annulus. Once the melting reaches close to 100%, the flow may become approximately identical to that seen in single-phase electro-convection in an annulus and the mean temperature reaches a steady state.

FIG. 10 is a diagram of liquid fraction, velocity, and temperature fields for natural convection melting in a vertical orientation for the LHTES unit 100. The liquid fraction, velocity, and temperature fields for pure natural convection melting (V=0 kV) are illustrated with respect to time in FIG. 10. The contours are presented in three circular planes at 0, 75, and 150 mm along the length and one central rectangular plane. Initially, the domain may be filled with PCM in a solid state, and there may be no fluid motion. Thermal conduction may dominate the heat transfer process. As time increases, the melting may begin, and liquid PCM may occupy the top part of the domain. The melt interface 1002 may not be vertical. The melt interface 1002, or the solid-liquid interface, may be closer to the hot tube wall at the bottom of the LHTES unit but may advance faster in the top part of the LHTES unit 100. This is illustrated in the liquid fraction contour, or shape of the melt interface 1002, at t=2500 s. The hotter fluid moves towards the top part of the LHTES unit 100, and the melting in the top region may be faster than in the bottom end of the LHTES unit 100. The corresponding velocity and temperature fields are presented in (b), (c), (c), (f), (h), and (i) of FIG. 10.

The temperature gradient may be lower at the top part of the LHTES unit 100, and the natural convective flow may be in the lower part of the domain where the temperature gradient is seen. Approximately the same flow pattern continues at t=5000 s. At t=5000 s, the volume of liquid in the domain increases, but the flow pattern and temperature field may approximately qualitatively follow the same trend as t=2500 s. The melting in vertical orientation (without an electric field) can be characterized by an inclined melt interface 1002, a higher temperature gradient, and a stronger natural convection in the lower region of the LHTES unit 100. At t=5000 s, approximately 80.1% volume fraction is reached in case 1 without an electric field.

FIG. 11 is a diagram of liquid fraction, velocity, and temperature fields for natural convection melting in a horizontal orientation for the LHTES unit 100. The instantaneous contours at time t=5, 2500, and 5000 s are illustrated in (a)-(c), (d)-(f), and (g)-(i), respectively, of FIG. 11. The natural convective fluid motion occurs above the surface of the hot tube. The heat transfer below the hot tube may be limited to thermal conduction. At t=2500 s, the melting occurs in the region above the hot tube, and the velocity distribution is also in the top region. In the horizontal orientation of the LHTES unit 100, the initial melting above the hot tube is faster due to natural convection. In later stages, the natural convection becomes weaker, and the melting process observed below the hot tube may be limited to thermal conduction. This can be related to the variation in the melting curve in case 2 at V=0 kV (refer to the plot 604 of FIG. 6).

FIG. 12 is a diagram of liquid fraction, charge density, velocity, and temperature fields for applied voltage and charge injection in a vertical orientation for the LHTES unit 100. The liquid fraction, charge density, velocity, and temperature fields in case 1 at V=10 kV and C=10 (strong injection regime) are illustrated in FIG. 12. At t=5 s, the melting may not have started, and the domain may be filled with solid PCM. The charge density distribution gradually decreases from the high voltage tube surface to the grounded shell surface. There may be no fluid motion, and the temperature distribution may be due to heat conduction. The contours at t=1250 s are illustrated in (e)-(g) of FIG. 12. The melt interface 1202 may be approximately vertical in case 1 in the presence of unipolar charge injection induced EHD flow. The diameter of the melted region in the bottom region may be marginally smaller, if at all, than in the top region. The charge density contours show charge void regions, a characteristic feature of electro-convection.

The charge density patterns in the three circular planes may not be identical because of the minor differences in the diameter of the melted region. However, electro-convection due to the charge injection from the hot tube surface can lead to fluid motion throughout the height of the LHTES unit 100 and can result in approximately uniform melting along the height of the LHTES unit 100. The temperature contour can indicate uniform heat transfer along the height as dictated by the electro-convection taking place along the height of the hot tube surface. At t=2500 s, 100% melting is achieved. Without an electric field, 100% melting may not be reached even at t=5000 s (refer to FIG. 10).

FIG. 13 is a diagram of contours of velocity distribution in the vertical orientation of the LHTES unit 100. 3D contours of velocity distribution in case 1 at 50% volume fraction for natural convection melting (left) and EHD assisted melting (right) are illustrated in FIG. 13. The velocity of the flow may be more substantial in the conical liquid region without an electric field.

FIG. 14 is a diagram of liquid fraction, charge density, velocity, and temperature fields for applied voltage and charge injection in a horizontal orientation for the LHTES unit 100. In some examples, unipolar charge injection from the hot tube surface can lead to a uniform cylindrical velocity field distributed throughout the height of the LHTES unit 100. The velocity magnitude in the presence of electroconvection can be approximately ten times higher than the pure natural convection melting. The distribution of liquid fraction, charge density, velocity, and temperature for the LHTES unit 100 in case 2 (horizontal orientation) at V=10 kV and C=10 (strong injection regime) is illustrated in FIG. 14. Electro-convection can lead to uniform fluid motion around the hot tube surface along the length of the LHTES unit 100. In the natural convection melting process in case 2, the melting occurs only in the region above the hot tube surface. Electro-convection can lead to uniform melting around the hot tube surface.

FIG. 15 is a diagram of contours of velocity distribution in the horizontal orientation of the LHTES unit 100. The 3D velocity distribution at 50% liquid fraction in case 2 without electric field (left) and with electric field (right) is illustrated in FIG. 15. The fluid motion may occur in the region above the hot tube in the case of pure thermal convection melting without any electric field. The case with charge injection shows uniform fluid motion around the hot tube surface. The electro-convection can result in a uniform flow field along the hot tube in cases 1 and 2. Furthermore, the velocity distribution is higher in the cases with the electric field. The melting process can become more uniform and faster due to charge injection induced EHD flow motion.

The performance of the LHTES unit 100 in both cases 1 and 2 and at all injection strengths (0≤C≤10) and applied voltages V=0, 5, and 10 kV can be analyzed. The charging performance can be evaluated in terms of charging time and enhancement ratio. The charging time is the total time each case takes to reach 80% melting. The enhancement ratio is the ratio of charging times taken by pure natural convection melting to the charging time in the presence of electric field (refer to Equation 11). The maximum liquid fraction achieved in 5000 s for the cases without electric field is 80%. The time to reach 80% melting is set as a benchmark to evaluate the charging performance. Case 1 and 2 without electric field may take 1.35 and 1.38 hours, respectively, to reach 80% melting. The time taken in case 2 is longer than in case 1 since the melting at later stages in the lower part of case 2 may be due to thermal conduction. For the cases with electric field, the charging time to reach 80% melting is lower than the pure natural convection melting.

Irrespective of the charge injection strength and the applied voltage, EHD flow may lead to shorter melting times. Increasing the applied voltage shortens the charging time within a given charge injection regime. For instance, when C=1, the total charging time to reach 80% melting is 0.71 hours and 0.53 hours at V=5 and 10 kV, respectively. In some examples, higher values of C indicate stronger charge injection regimes. Higher charge injection from the high voltage electrode can result in a stronger Coulomb force. The intensity of EHD flow is higher at higher values of C. However, the Coulomb force is a combined function of charge density and electric field intensity. In cases 1 and 2, the medium injection regime with V=10 kV may have a shorter charging time than V=5 kV in C=10. This indicates that the EHD flow induced in the medium injection regime at applied electric potential V=10 kV is more potent than that observed in C=10 and V=5 kV. Unipolar charge injection induced EHD flow can lead to an enhancement factor greater than 1. The EHD flow contributes positively to accelerating the charging process.

In case 1, a maximum enhancement ratio of approximately 2.62 may be achieved at strong injection regime with applied voltage V=10 kV. Likewise, in case 2, V=10 kV with C=10, a maximum enhancement ratio of approximately 2.70 may be achieved. At any given applied voltage and charge injection regime, case 2 may have slightly shorter charging times than case 1. The thermal convection in case 2 may be weaker than in case 1. It may be easier for the EHD flow to govern the melting process and to lead to shorter melting times. A summary of charging time and enhancement ratios is illustrated in Table 3 produced below.

TABLE 3
Summary of charging time and enhancement ratio at 0 <
C < 10 and 0 < V < 10 kV.

	Voltage	Charging	Enhancement
Injection Strength	[kV]	Time * [hours]	ratio

Case 1	No EHD	—	1.35	—
	C = 0.1	5	1.15	1.17
		10	0.77	1.75
	C = 1	5	0.71	1.90
		10	0.53	2.55
	C = 10	5	0.65	2.08
		10	0.52	2.60
Case 2	No EHD	—	1.38	—
	C = 0.1	5	1.13	1.22
		10	0.72	1.92
	C= 1	5	0.70	1.97
		10	0.52	2.65
	C = 10	5	0.60	2.30
		10	0.52	2.65
*Charging time to reach 80% liquid fraction.

FIG. 16 is a set of plots 1600 of power stored for the LHTES unit 100 at different charge injection regimes and at different applied voltages. As illustrated in FIG. 16, the set of plots 1600 can include a plot 1602 corresponding with the vertical orientation for the LHTES unit 100, and the set of plots 1600 can include a plot 1604 corresponding with the horizontal orientation for the LHTES unit 100. The plot 1602 can include a horizontal axis 1606 that represents the different charge injection regimes and the different applied voltages, and the plot 1602 can include a vertical axis 1608 that represents an amount of power stored by the LHTES unit 100. The plot 1604 can include a horizontal axis 1610 that represents the different charge injection regimes and the different applied voltages, and the plot 1604 can include a vertical axis 1612 that represents an amount of power stored by the LHTES unit 100.

The net power stored (see Equation 12) in each case when reaching 80% melting can be illustrated in FIG. 16. A higher applied voltage at a given charge injection regime can lead to higher power storage. Likewise, stronger charge injection can lead to higher power storage at a given applied voltage. The increase in power storage at V=10 kV at C=10 may be marginally higher than that in V=10 kV and C=1. The maximum heat transfer and melting rate can be limited by the gap between the shell and tube of the LHTES unit 100. The role of EHD flow in the total power stored is higher in case 2. A maximum of approximately 93.70 kW and approximately 115.85 KW net power storage can be achieved in case 1 and case 2, respectively.

The increase in charging performance of the LHTES unit 100 by applying an electric field may be associated with additional power requirements for generating EHD flow. The extra electric power consumed to initiate charge injection and the resultant EHD flow can be calculated using Equation 13. The average total electric current consumption (refer to Equation 14) for a period of 5000 s through the emitter electrode (hot tube connected to high electric potential) surface calculated as per Equation 14 for each case is listed in Table 4 produced below.

TABLE 4
Summary of time averaged electric power consumption
to generate EHD flow at 0 < C < 10 and
V = 5 over a time period of 5000 s..

	Applied
Injection Strength	Voltage [kV]	Electric power consumption [mW]

Case 1	C = 0.1	5	0.6069
		10	2.0138
	C = 1	5	0.8257
		10	3.5673
	C = 10	5	3.6181
		10	5.1432
Case 2	C = 0.1	5	0.8678
		10	2.0932
	C = 1	5	0.9500
		10	3.9450
	C = 10	5	5.6180
		10	9.9950

The additional electric power consumption ranges from a few μW in the weak injection regime to a few mW in the strong injection regime. However, the additional contribution by EHD flow to the net power storage is in the range of several kW. Thus, the increased power storage may be at least six orders of magnitude higher than the additional power expended to generate EHD flow. Unipolar charge injection induced electrohydrodynamic flow is a viable and efficient approach to increase the charging performance of the LHTES unit 100. FIG. 17 is a set of plots 1700 of liquid fraction over time for the LHTES unit 100 in

different orientations, according to some embodiments. The set of plots 1700 can include a plot 1702 corresponding with the vertical orientation of the LHTES unit 100, and the set of plots 1700 can include a plot 1704 corresponding with the horizontal orientation of the LHTES unit 100. The plot 1702 can include a horizontal axis 1706 representing time and can include a vertical axis 1708 representing liquid fraction, and the plot 1704 can include a horizontal axis 1710 representing time and can include a vertical axis 1712 representing liquid fraction.

In some embodiments, FIG. 17 represents the temporal evolution of liquid fraction over a total simulation period of 5 hours for vertical and horizontal orientations, respectively, under varying gravity levels (1 g, 0.5 g, and 0 g) and applied electric field conditions (0 kV and 10 kV). In the absence of electric field (0 kV), the melting process may be governed by buoyancy-driven natural convection, which strength is directly proportional to the gravitational acceleration. In the vertical configuration ((a) of FIG. 17), at 1 g, the liquid fraction reaches approximately 84% at 5 hours. When gravity is reduced to 0.5 g, the liquid fraction decreases to about 82%, while under microgravity (0 g), only around 80% melting is achieved within the same period. This reduction in melting rate with decreasing gravity may indicate the diminishing intensity of buoyancy forces, which arise from the small temperature-induced density differences within the liquid PCM and act vertically in the direction of gravity. At 1 g, these buoyancy forces drive an upward motion of warmer, lighter liquid PCM adjacent to the hot inner wall, coupled with a compensating downward flow of cooler liquid near the outer shell wall. This can establish a coherent circulation loop that can enhance convective transport of thermal energy away from the heated surface, accelerating the propagation of the melting front.

As gravity decreases, the magnitude of the buoyancy force may diminish, weakening the induced flow field and reducing convective mixing efficiency. In the 0 g case, the complete absence of gravity eliminates buoyancy altogether, resulting in a purely conduction-controlled melting process in which heat transfer is governed by thermal diffusion across the stationary PCM. A similar gravity-dependent trend can be present for the horizontal orientation ((b) of FIG. 17), but with systematically higher melting rates compared to the vertical orientation for each gravity level. For example, at 1 g, the horizontal case reaches approximately 86.5% liquid fraction at 5 hours, while at 0.5 g and 0 g, the liquid fractions reach about 84.5% and 80%, respectively. This superior melting performance in the horizontal orientation is attributable to the asymmetric development of natural convection. In this configuration, buoyancy acts perpendicular to the tube axis, causing the lighter liquid PCM adjacent to the hot tube to rise vertically and accumulate in the upper portion of the annular cavity. This can generate a large, stable, single convection cell in the upper domain that intensifies mixing and promotes more effective heat transfer in the region directly above the hot surface. While the lower portion of the domain remains largely stagnant and conduction-dominated, the dominant convection cell in the upper region can result in an enhancement in the overall melting rate compared to the vertical orientation, particularly during the early and intermediate stages of melting.

However, the influence of both gravity and orientation becomes completely negligible once the electric field is applied. Upon introduction of the electric field (10 kV), unipolar charge injection induces strong electrohydrodynamic (EHD) flow, and the melting dynamics become governed by electro-convective mechanisms. In both vertical and horizontal orientations, across all gravity levels, the EHD-assisted liquid fraction curves collapse onto nearly identical trajectories, achieving complete melting (100% liquid fraction) within approximately 0.55 to 0.6 hours. The applied electric potential establishes a radial electric field between the high-voltage emitter electrode (inner wall) and grounded shell, driving unipolar charge transport that generates Coulomb body forces. These forces produce radial and rotational electro-convective vortices throughout the PCM domain, vigorously mixing the liquid phase and rapidly distributing heat away from the heated surface into the bulk PCM. The magnitude of these electric forces exceeds the available buoyancy forces, even at full gravity, thereby dominating the flow field entirely and rendering gravity effects irrelevant. Moreover, because the applied electric field and the resulting flow structures are axisymmetric about the cylindrical electrode, the melting performance remains identical in both vertical and horizontal orientations regardless of the direction of gravity.

FIG. 18 is a set of plots 1800 of maximum liquid velocity over time for the LHTES unit 100 in different orientations, according to some embodiments. The set of plots 1800 can include a plot 1802 corresponding with the vertical orientation of the LHTES unit 100, and the set of plots 1800 can include a plot 1804 corresponding with the horizontal orientation of the LHTES unit 100. The plot 1802 can include a horizontal axis 1806 representing time and can include a vertical axis 1808 representing maximum liquid velocity, and the plot 1804 can include a horizontal axis 1810 representing time and can include a vertical axis 1812 representing maximum liquid velocity.

In some embodiments, FIG. 18 illustrates the temporal evolution of maximum velocity within the computational domain for vertical and horizontal orientations, respectively, under varying gravity levels and applied electric field conditions. The progression of maximum velocity can be a direct indicator of the strength and development of convective motion during the melting process, governed by natural convection, electrohydrodynamic (EHD) forces, or their combined interaction depending on the operating conditions. In the absence of electric field (0 kV), the flow field is controlled by buoyancy-driven natural convection. At 1 g and 0.5 g, the maximum velocity increases progressively during the initial and intermediate stages of melting, as the expanding liquid region allows larger convective cells to form and strengthen over time. The rate of this increase is directly correlated with the magnitude of the gravitational body force. For example, higher gravity results in stronger buoyancy forces, promoting more vigorous natural convection and hence higher flow velocities. Under microgravity (0 g), buoyancy forces are absent, and the velocity remains nearly flat throughout, indicating purely conduction-dominated heat transfer with negligible fluid motion. The inset plots in both figures present zoomed views. The insect plots illustrate the gradual onset and strengthening of buoyancy-driven flow as melting progresses at finite gravity levels, while the zero gravity case exhibits zero velocity. As the phase change approaches completion, the remaining solid phase becomes increasingly limited, reducing the thermal gradients within the liquid PCM. With diminished temperature differences available to sustain buoyancy forces, the driving potential for natural convection weakens, leading to a gradual reduction in fluid motion. This decay in velocity becomes more pronounced once the molten region saturates most of the domain and thermal homogenization takes place. This velocity decline may occur earlier and more sharply in the horizontal orientation compared to the vertical configuration. The horizontally oriented configuration can generate a dominant single convection cell in the upper portion of the domain, while the lower region remains largely stagnant. As melting progresses and the upper region fully liquefies, the driving thermal gradients across the active convection cell diminish more rapidly than in the vertical orientation, where a larger vertical fluid column sustains broader recirculation for a longer duration. As a result, in the horizontal orientation, the maximum velocity peak is reached sooner, followed by an earlier and sharper velocity decline toward the end of the melting process.

When the electric field is applied (10 kV), the evolution of maximum velocity changes. Upon activation of unipolar charge injection, the electrohydrodynamic (EHD) forces can generate radial electric body forces that induce development of electro-convective vortices. The vortices establish fluid motion throughout the liquid PCM, independent of gravity and orientation. The velocity rises during the early transient period and remains sustained until complete melting is achieved. Under EHD conditions, the maximum velocity curves may exhibit random fluctuations throughout the melting process such as after the initial transient phase. The irregular oscillations are indicative of the inherently unsteady and chaotic nature of electro-convective flow structures. The complex interaction between charge injection, Coulomb forces, and continuously evolving thermal gradients leads to dynamically changing flow patterns, which can be characterized by vortex merging, splitting, and fluctuating vortex intensity. The chaotic behavior can indicate electrohydrodynamic (EHD) convection under unipolar charge injection, such as at high electric field strengths, and contributes to the rapid and highly effective mixing that dramatically accelerates the melting process.

FIG. 19 is a set of plots 1900 of mean Nusselt number over time for the LHTES unit 100

in different orientations, according to some embodiments. The set of plots 1900 can include a plot 1902 corresponding with the vertical orientation of the LHTES unit 100, and the set of plots 1900 can include a plot 1904 corresponding with the horizontal orientation of the LHTES unit 100. The plot 1902 can include a horizontal axis 1906 representing time and can include a vertical axis 1908 representing mean Nusselt number, and the plot 1904 can include a horizontal axis 1910 representing time and can include a vertical axis 1912 representing mean Nusselt number.

In some embodiments, FIG. 19 illustrates the temporal evolution of the mean Nusselt number calculated along the inner hot wall for vertical and horizontal orientations, respectively, under varying gravity levels and applied electric field conditions. The Nusselt number reflects the convective heat transfer intensity at the inner wall, capturing both the influence of natural convection driven by buoyancy and electro-convective enhancement due to electric field application. In both orientations, during the early phase of melting (approximately first 0.5 hours), a sharp initial peak in Nusselt number can be present for all cases. The initially high Nusselt number corresponds to the purely conduction-dominated heat transfer stage when the PCM adjacent to the hot wall is in solid phase. The temperature gradient near the hot wall is steep due to the small thickness of the early melting layer, leading to local heat fluxes and elevated Nusselt numbers. As melting progresses, the temperature gradients begin to relax as the liquid layer grows, leading to a decline in the Nusselt number after the initial peak. For the cases without electric field (0 kV), after the initial stage, buoyancy-driven natural convection progressively develops and begins to contribute to heat transfer as liquid accumulates. In the vertical orientation ((a) of FIG. 19), as gravity increases from 0 g to 0.5 g to 1 g, the mean Nusselt number increases, reflecting stronger convective heat transfer via increasing buoyancy forces. In 1 g, the buoyancy-induced circulation becomes well-established, continuously enhancing heat transport across the liquid phase, and resulting in a relatively higher sustained Nusselt number compared to reduced gravity cases. At 0.5 g, the convective enhancement is present but weaker, while under microgravity (0g), natural convection is absent, and heat transfer remains conduction-dominated throughout the process, resulting in the lowest sustained Nusselt number after the initial decay. The inset plot in FIG. 19 provides a zoomed view of the post-conduction behavior for natural convection cases, clearly highlighting the sustained separation between the different gravity levels.

As melting continues and the liquid volume increases, the Nusselt number in natural convection cases declines. This decline can be attributed to a reduction in temperature gradients at the hot wall as the bulk PCM temperature rises and thermal stratification reduces. The decreasing temperature difference between the hot wall and the adjacent liquid weakens the driving potential for buoyancy and the resulting convective flow strength, leading to declining convective heat transfer rates during the later stages of melting. A similar behavior can be present in the horizontal orientation ((b) of FIG. 19), though with some quantitative differences. The horizontal configuration has a higher Nusselt number for corresponding gravity levels compared to the vertical orientation. This enhancement is a result of the asymmetric natural convection pattern that develops in the horizontal configuration. For example, buoyancy can cause the melted liquid PCM to rise and accumulate preferentially above the horizontally aligned hot tube, forming a single recirculation cell in the upper region. This localized enhancement of convective transport near the upper part of the hot wall maintains stronger temperature gradients for a longer duration during intermediate melting stages, resulting in higher mean Nusselt numbers relative to the vertical configuration.

The behavior may change with application of electric field (10 kV). For vertical and horizontal orientations, and across all gravity levels, the EHD-assisted cases exhibit a different evolution profile. After the initial conduction-dominated peak, the mean Nusselt number rises again due to the establishment of electrohydrodynamic convection. The Coulomb forces generated by unipolar charge injection initiate electro-convective vortices that actively mix the liquid PCM and enhance convective transport across the entire domain. As a result, the Nusselt number remains significantly elevated and sustained throughout the melting process, demonstrating the dominant role of electro-convective mixing in maintaining high convective heat transfer rates. The Nusselt number curves for each gravity level collapse onto a nearly identical profile under EHD conditions, indicating that the heat transfer enhancement provided by electrohydrodynamic (EHD) convection is independent of gravity and orientation. The gravity-neutral establishment of strong radial electro-convective vortices ensures uniform heat distribution across the liquid domain and maintains high temperature gradients at the hot wall throughout the phase change process. This sustained heat transfer performance directly correlates with the reduced melting times observed in the liquid fraction evolution results (refer to FIG. 17).

The LHTES unit 100 can be assisted with active EHD flow to enhance the charging performance of the LHTES unit 100. EHD flow due to unipolar charge injection can alter the flow structure, can increase the flow intensity, and can increase the fluid mixing and convective heat transfer within the melted PCM. The EHD flow intensity and heat transfer can be a direct function of charge injection strength. That is, stronger charge injection can lead to a more substantial EHD flow. In case 1 (vertical orientation), the presence of EHD flow can lead to uniform melting along the height of the LHTES unit 100. Unlike the inclined melt interface without an electric field, the melting is uniform from top to bottom of the LHTES unit 100 assisted by EHD flow. In case 2 (horizontal orientation), the natural convection may be limited to the region above the hot tube.

Thus, the melting may be limited to the upper region without an electric field. The electric field generates electro-convective flow around the hot tube, and the melting takes place uniformly around the hot tube.

Weaker buoyancy force in horizontal orientation may favor the action of EHD flow and may lead to slightly higher performance than case 1. The maximum additional power required to generate EHD flow may be approximately a few milliwatts, whereas the additional power storage achieved by generating EHD flow is approximately several kilowatts. Comparing the acceleration in charging time, higher power storage capacity, and the additional electric power consumption may show that horizontal orientation of the LHTES unit 100 with EHD flow generated in medium charge injection regime may have the highest charge enhancement (within the parameter space considered herein).