The perpendicular impingement of a gas stream on an electric arc, a configuration known as the arc in crossflow, is of primary relevance in the study of plasma–gas interactions as well as in industrial applications such as circuit breakers and wire-arc spraying. The flow dynamics in the arc in crossflow are the result of coupled fluid-thermal-electromagnetic phenomena accompanied by large property gradients, which can produce significant deviations from Local Thermodynamic Equilibrium (LTE) among electrons and gas species. These characteristics can lead to the establishment of distinct flow regimes depending on the relative values of the controlling parameters of the system, such as inflow velocity, arc current, and inter-electrode spacing. A two-temperature non-LTE model is used to investigate the arc dynamics and the establishment of flow regimes in the arc in crossflow. The plasma flow model is implemented within a nonlinear Variational Multiscale (VMS) numerical discretization approach that is less dissipative and, hence, better suited to capture unstable behavior than traditional VMS methods commonly used in computational fluid dynamics simulations. The Reynolds and the Enthalpy dimensionless numbers, characterizing the relative flow strength and arc strength, respectively, are chosen as the controlling parameters of the system. Simulation results reveal the onset of dynamic behavior and the establishment of steady, periodic, quasi-periodic, and chaotic or potentially turbulent regimes, as identified by distinct spatiotemporal fluctuations. The computational results reveal the role of increasing the relative arc strength on enhancing flow stability by delaying the growth of fluctuating and unstable flow behavior.

## I. INTRODUCTION

### A. The arc in crossflow

The arc in crossflow is a plasma configuration commonly encountered in industrial applications such as wire arc spraying and low-voltage circuit breakers, which involves the perpendicular impingement of a stream of gas onto an electric arc. The gas flow impingement causes the arc to bend and elongate downstream. The effects of the arc–gas flow interaction affect, for example, wire-arc spraying in terms of *process parameters*, such as particle dispersion,^{1,2} droplet formation,^{3} particle trajectories,^{2,4} atomization,^{1} and oxidation due to gas entrainment;^{5} as well as *electrode parameters*, such as electrode erosion,^{6,7} splitter erosion,^{8,9} and asymmetric melting.^{10} These phenomena may seem inter-dependent on each other (e.g., asymmetric melting can cause higher electrode erosion). Nevertheless, they are intrinsically the result of the strong plasma–gas flow interaction. This interaction is characterized by highly coupled and complex physical phenomena such as the imbalance between electromagnetic (Lorentz) and fluid dynamic drag forces, large temperature and density gradients, and marked deviations from Local Thermodynamic Equilibrium (LTE) between electrons and heavy-species (atoms, ions, molecules).^{11} In addition, the plasma–gas flow interaction can lead to the dynamic behavior of the arc, manifested by significant temperature and voltage fluctuations,^{12–14} which can markedly influence process characteristics.

The present work addresses the characterization of the intrinsic dynamics and the establishment of associated flow regimes in the arc in crossflow. These flow regimes depend on the relative values of the controlling parameters of the system, such as inflow velocity, arc current, and inter-electrode spacing. The present study adopts the Reynolds number, which characterizes the ratio of inertial to viscous forces, and the Enthalpy number, which characterizes the ratio of gas flow energy to electrical energy, as the main controlling parameters of the system. These controlling parameters are used to evaluate *qualitative* characteristics (i.e., the establishment of regimes such as steady, periodic, and quasi-periodic, as indicated by spatiotemporal property fluctuations) as well as *quantitative* characteristics (i.e., magnitude and frequency of fluctuations) of the arc in crossflow. The unraveling of these intrinsic flow characteristics elucidates the fundamental aspects of plasma–gas interactions and provides understanding that can aid in equipment and process design in arc in crossflow-based applications.

A schematic representation of the arc in crossflow model is presented in Fig. 1. The main parameters of the system are working gas type, total current *I _{tot}*, inlet gas velocity

*U*, and electrode gap

_{i}*H*. The model is applied to an arc discharge in pure argon. As indicated in Fig. 1(a), the interaction of plasma–gas flow interaction causes the bending of the arc, which consequently leads to large property gradients and coupled thermal-fluid-electromagnetic phenomena.

^{15,16}Figure 1(b) presents the computational domain for the solution of the model, including the domain's geometric parameters and boundaries. The present study builds on the prior computational characterization of the arc in crossflow presented in Refs. 11 and 17.

### B. Computational study of plasma dynamics

Plasma flow dynamics have been studied more thoroughly, compared to the arc in crossflow, in configurations such as plasma torches and pin-to-plate setups due to their direct relevance to industrial applications (e.g., plasma spraying and welding, respectively).^{2,18–24} The dynamics involved in these configurations are significantly involved, with several inter-related physical phenomena, such as continuous change in the arc length,^{25} azimuthal rotation of the arc,^{26} arc root movement,^{10} etc. Such dynamic interplay leads, for example, to different arc modes in non-transferred arc plasma torches.^{18,27} In contrast to the arc systems above, the arc in crossflow is significantly simpler, yet it acts as a canonical configuration for the study of plasma–gas flow interactions,^{28,29} offering fundamental understanding of arc plasma dynamics.

Computational investigations offer distinct advantages compared to experimental diagnostics, such as concurrent estimation of distributed solution (pressure, temperature, velocity, electric potential, etc.) and derived (magnetic field, Lorentz force, current density, etc.) fields. Computational studies of plasma dynamics necessarily require utilizing time-dependent models (to resolve temporal fluctuations) in three-dimensional (3D) spatial domains (to track the spatial evolution of fluctuations), which make them computationally expensive. The computational expense of the study of plasma–gas flow interactions is exacerbated considering that the associated large temperature, density, and property gradients along with complex multi-physical phenomena (such as fluid dynamics, heat transfer, chemical kinetics, and electromagnetics) ultimately leading to turbulence^{13}—the outmost computationally demanding fluid flow phenomena to simulate. Diverse research works have addressed the computational cost associated with the study of plasma flow dynamics at the expense of using simplified models, such as assuming the plasma to be in a LTE state.^{25,30} However, the LTE assumption is valid only in the plasma core, and significant deviations are expected in strong plasma–gas interactions, thermal boundary layers,^{2,18} etc. In this regard, non-LTE (NLTE) plasma flow models^{12,31} can be better suited than LTE models^{32,33} for the study of plasma–gas flow interactions. Moreover, the need to resolve dynamic behavior and instabilities imposes additional requirements in terms of the numerical resolution of the computational model, such as the use of fine temporal and spatial discretizations (computational grids or meshes, and time steps) and low-dissipation schemes.^{34}

### C. Outline

The article is organized as follows: Sec. II presents the mathematical model, the numerical approach, and the computational model setup. The controlling parameters of the system are presented in Sec. III. Section IV describes the approach used to identify and characterize the arc dynamics, the role of the numerical approach on the obtained solutions, and the effect of controlling parameters on the arc dynamics. Results of the identification of arc in crossflow regimes and their characteristics are presented in Sec. V. Summarizing and concluding remarks are presented in Sec. VI.

## II. COMPUTATIONAL MODEL

### A. Mathematical model

The plasma is considered as an optically thin, quasi-neutral, non-magnetized, compressible, reactive, and electromagnetic fluid in chemical equilibrium and thermal nonequilibrium (two-temperature NLTE). Ion diffusion, Hall currents, and electrode sheaths are neglected. The mathematical model is constituted by the equations of conservation of total mass, mass-average momentum, thermal energy of heavy-species, and thermal energy of electrons, and the equations for electric charge conservation and of magnetic induction. These equations are listed in Table I in the transient-advective-diffusive-reactive (TADR) form.

Equation . | Transient . | Advective . | Diffusive . | Reactive . |
---|---|---|---|---|

Mass conservation | $\u2202t\rho $ | $u\u22c5\u2207\rho +\rho \u2207\u22c5u$ | 0 | 0 |

Momentum conservation | $\rho \u2202tu$ | $\rho u\u22c5\u2207u+\u2207p$ | $\u2207\u22c5\mu (\u2207u+\u2207uT)\u2212\u2207\u22c5(23\mu (\u2207\u22c5u)\delta )$ | $Jq\xd7B$ |

Thermal energy heavy-species | $\rho \u2202thh$ | $\rho u\u22c5\u2207hh$ | $\u2207\u22c5(\kappa hr\u2207Th)$ | $Dtph+Keh(Te\u2212Th)\u2212\tau :\u2207u$ |

Thermal energy electrons | $\rho \u2202the$ | $\rho u\u22c5\u2207he$ | $\u2207\u22c5(\kappa e\u2207Te)$ | $Dtpe\u2212Keh(Te\u2212Th)\u22124\pi \epsilon r+Jq\u22c5(E+u\xd7B)+5kB2eJq\u22c5\u2207Te$ |

Charge conservation | 0 | 0 | $\u2207\u22c5(\sigma \u2207\varphi p)\u2212\u2207\u22c5(\sigma u\xd7(\u2207\xd7A))$ | 0 |

Magnetic induction | $\mu 0\sigma \u2202tA$ | $\mu 0\sigma \u2207\varphi p\u2212\mu 0\sigma u\xd7(\u2207\xd7A)$ | $\u22072A$ | 0 |

Equation . | Transient . | Advective . | Diffusive . | Reactive . |
---|---|---|---|---|

Mass conservation | $\u2202t\rho $ | $u\u22c5\u2207\rho +\rho \u2207\u22c5u$ | 0 | 0 |

Momentum conservation | $\rho \u2202tu$ | $\rho u\u22c5\u2207u+\u2207p$ | $\u2207\u22c5\mu (\u2207u+\u2207uT)\u2212\u2207\u22c5(23\mu (\u2207\u22c5u)\delta )$ | $Jq\xd7B$ |

Thermal energy heavy-species | $\rho \u2202thh$ | $\rho u\u22c5\u2207hh$ | $\u2207\u22c5(\kappa hr\u2207Th)$ | $Dtph+Keh(Te\u2212Th)\u2212\tau :\u2207u$ |

Thermal energy electrons | $\rho \u2202the$ | $\rho u\u22c5\u2207he$ | $\u2207\u22c5(\kappa e\u2207Te)$ | $Dtpe\u2212Keh(Te\u2212Th)\u22124\pi \epsilon r+Jq\u22c5(E+u\xd7B)+5kB2eJq\u22c5\u2207Te$ |

Charge conservation | 0 | 0 | $\u2207\u22c5(\sigma \u2207\varphi p)\u2212\u2207\u22c5(\sigma u\xd7(\u2207\xd7A))$ | 0 |

Magnetic induction | $\mu 0\sigma \u2202tA$ | $\mu 0\sigma \u2207\varphi p\u2212\mu 0\sigma u\xd7(\u2207\xd7A)$ | $\u22072A$ | 0 |

In Table I, $\u2202t=\u2202/\u2202t$ is the partial time derivative, $\u2207$ and $\u2207\u22c5$ are the gradient and divergence operators, respectively; *ρ* is mass density, *p* is the pressure, **u** is the mass-averaged velocity, *μ* is the dynamic viscosity, $T$ is the transpose operator, **δ** is the Kronecker delta tensor; **J**_{q} is the electric current density, and **B** is the self-induced magnetic field; *h _{h}* and

*h*are the enthalpies and

_{e}*T*and

_{h}*T*are the temperatures of the heavy-species and electrons, respectively; $\kappa hr$ is the heavy-species translational-reactive thermal conductivity, $\kappa e$ is the translational electron thermal conductivity, $Dtp$ is the pressure work, where $Dt\u2261\u2202t+u\u22c5\u2207$ represents the total derivative, and

_{e}*p*=

*p*represents the total pressure, with

_{h}+ p_{e}*p*and

_{h}*p*the heavy-species and electron pressure, respectively;

_{e}*K*is the energy exchange coefficient between heavy-species and electrons,

_{eh}**τ**is the stress tensor for a Newtonian fluid,

*ε*is effective net radiative emission coefficient;

_{r}**E**is the

*real*electric field,

*k*is the Boltzmann constant,

_{B}*e*is the elementary electric charge,

*ϕ*is the

_{p}*effective*electric potential, and

**A**is the magnetic vector potential.

The set of TADR equations is treated in a monolithic form as a single system given by

where $R$ is the residual vector, **Y** the vector of unknowns, **A _{0}**,

**A**

_{i},

**K**

_{ij},

**S**,

_{1}**S**are the coefficient matrices that describe the different transport processes,

_{0}*i*and

*j*are spatial indices, and Einstein's convention of repeated indexes is used. Equation (1) is solved for the set of primate variables

The model is complemented by the definition of thermodynamic (*h _{h}*,

*h*, etc.) and transport (

_{e}*μ*, $\kappa e$, etc.) material properties (data used from Ref. 35 and references within) as well as by constitutive relations [$Jq=\sigma (E+u\xd7B)$, $B=\u2207\xd7A$, etc.]. A detailed description of the model is presented in Refs. 11 and 35.

The evaluation of the above NLTE model against a LTE model^{32} and experimental results has been presented in Ref. 11. The investigation in Ref. 11 showed that results from the two-temperature model present better agreement with experimental observations than the LTE model, while providing quantification of the deviation between heavy-species and electron temperatures. Moreover, the use of the LTE assumption generally leads to overprediction of the arc voltage drop across the arc. For example, in twin torch systems^{36} and plasma torches^{31} the arc voltage drop obtained with an LTE model is significantly higher than that using an NLTE one. Therefore, the present work focuses on the exploration and identification of flow regimes using a NLTE model only.

### B. Numerical model

Computational thermal plasma flow simulations inherently involve large and nonlinear variations in solution fields and material properties, which make the associated models numerically stiff, and, hence, require the use of robust discretization and numerical techniques for their solution. Due to their inherent advantages, Finite Element Methods (FEMs) are widely used in diverse fields, including plasma flow modeling. Among FEMs, stabilized and particularly, Variational Multiscale (VMS) methods have demonstrated to be particularly effective dealing with highly multiscale phenomena.^{37} VMS methods have been used to solve scalar transport, radiation transport, incompressible, compressible, reactive, magnetohydrodynamics, and turbulent flow problems,^{35} as well as nonequilibrium plasma flows.^{11,12,27,31,38–41}

Recently, Modirkhazeni and Trelles^{42} have presented a VMS formulation that provides a better description of the intricate nonlinear coupling between the large-scales and the small-scales of the solution fields termed non-linear Variational Multiscale (VMS_{n}). (The large-scales are those actually captured by the discretization and the small-scales are unresolved sub-grid features.) VMS_{n} has demonstrated to be effective for the unified simulation of incompressible and compressible laminar and turbulent flow problems,^{43} as well as the flow from a non-transferred arc plasma torch,^{42} providing greater accuracy than a VMS approach. The greater accuracy of the VMS_{n} method has prompted its adoption in the present research of the arc in crossflow.

VMS methods start with a Galerkin formulation of the problem given by Eq. (1), i.e.,

where $\Omega $ is the spatial domain, **W** is the weight function, and $L$ is the transport operator defined from $R(Y)=LY\u2212S0$. Then, **Y** and **W** are decomposed into large- and small-scales, i.e., $Y=Y\xaf+Y\u2032$ and $W=W\xaf+W\u2032$. The use of traditional VMS methods assumes $L=L(Y)$$=L(Y\xaf+Y\u2032)\u2248L(Y\xaf)$, that is, the effect of the small-scales in the transport operator is neglected. The small-scales are subsequently calculated by an algebraic approximation as $Y\u2032\u2248\u2212\tau R(Y\xaf)$, where **τ** is the so-called *intrinsic time scales matrix* such that $\tau \u2248L\u22121$. Using this approximation, Eq. (3) can be treated as a function of large-scale *alone*. The VMS method's implementation and validation for the simulation of nonequilibrium thermal plasma flows, including a previous study of the arc in crossflow, are shown in Ref. 35.

The approximation $L(Y)\u2248L(Y\xaf)$ in traditional VMS methods is generally valid for many flow problems, but it is generally inadequate for problems in which the small-scales play a major role, such as in transitional and turbulent flows. The VMS_{n} method addresses the nonlinearity of the transport operator upfront, i.e., $L=L(Y\xaf+Y\u2032)$. This results in two equations to be solved concurrently, one for large- and one for small-scales each, i.e.,

where $*$ denotes the adjoint operator. The intrinsic time scale matrix is also a nonlinear function, i.e., $\tau =\tau (Y\xaf+Y\u2032)$. The nonlinearity of $L$ and **τ** requires the use of a nonlinear solution approach for the solution of Eq. (4). This is accomplished using a fixed-point (Picard) iteration given by

where *n* indicates the iteration counter and the iteration starts with $Y\u20320=0$. It is to be noted that if *n *=* *0, the VMS_{n} results in the VMS method, i.e., $L=L(Y\xaf)$, $\tau =\tau (Y\xaf),$ and $Y\u2032\u2248\u2009\u2212\tau (LY\xaf\u2212S\xaf0)$.

The final VMS_{n} formulation is given by concurrent iterative solution of the equations^{12,35}

and

where **N** is the basis function (i.e., FEM interpolation function) and the notation *f* (*a*;*b*) indicates that *a* is an argument of the equation and *b* is a parameter. In Eq. (6) a *discontinuity-capturing* operator is added to minimize the occurrence of spurious oscillations near unresolved high-gradient regions as well as to increase the robustness of the solution approach. Details about the numerical solution approach are found in Ref. 35.

VMS methods, by resolving the large-scales and modeling the small-scales of the solution, constitute coarse-grained simulation approaches, equivalent to Large Eddy Simulation (LES) methods for the simulation of turbulent flows.^{54–56} Yet, in contrast to traditional LES methods, which rely on Smagorinsky-like models to model small-scale interactions, VMS methods rely on a local approximation of small-scales problem (i.e., $Y\u2032\u2248\u2009\u2212\tau (LY\xaf\u2212S\xaf0)$). This characteristic makes VMS approaches particularly appealing for the computational modeling of highly nonlinear flow problems, such as plasma flows.

### C. Geometrical configuration and boundary conditions

The geometrical dimensions of the computational domain for the arc in crossflow in Fig. 1 are: *L _{x}* ×

*L*×

_{y}*L*= 2 mm × 10.1 mm × 37 mm. The electrodes are squares of size

_{z}*d*=

_{y}*d*= 3 mm. The inter-electrode distance (

_{z}*L*=

_{x}*H*) is the most important dimension for the arc in crossflow as it strongly influences the characteristics of the arc column. The dimensions

*L*and

_{y}*L*can be chosen freely as long as the plasma remains well within the computational domain (

_{z}*L*) and the extent of the domain does not alter the arc behavior (

_{y}*L*). The computational domain is discretized using a semi-structured tri-linear hexahedral elements grid of ∼250 000 nodes. This is in contrast to the simulations reported in Ref. 11, which used a significantly coarser grid of ∼85 000 nodes. Additionally, the current work used a model discretization based on the VMS

_{z}_{n}method described in Sec. II B; which is also in contrast to the simulations reported in Ref. 11 based on a classical VMS method. The higher numerical resolution and improved numerical discretization approach were required in the present work in order to resolve the onset of instabilities and to attain convergent solutions for a wider range of controlling parameters than those used in Ref. 11. A comparison between results obtained with the VMS and VMS

_{n}methods, for the same spatial discretization, is presented in Sec. IV C.

Table II lists the set of boundary conditions used to describe the problem. In Table II, *n* represents the outer normal direction to the corresponding boundary. The outflow_z boundary is set to be at atmospheric pressure, i.e., *p *=* p _{∞}* = 1.01325 × 10

^{5}Pa, while at all other boundaries a zero-pressure gradient normal to the surface is imposed. For the inflow boundary, the velocity profile is imposed such that it is maximum at the center and minimum at the edges, i.e.,

**u**= (

_{i}*u*)$T$ = (0 0

_{ix}u_{iy}u_{iz}*U*), where

_{i}*U*(

_{i}*x*) =

*U*[1 − (

_{imax}*x*/

*H*)

^{2}]. The no-slip condition (i.e.,

**u**=

**0**) is imposed over electrode surfaces and to the solid surfaces surrounding them (i.e., anode, cathode, and walls).

Boundary . | p
. | u
. | T
. _{h} | T
. _{e} | $\varphi p$ . | A
. |
---|---|---|---|---|---|---|

Inflow | $\u2202np=0$ | u = u _{i} | T = _{h}T_{0} | T = _{e}T_{0} | $\u2202n\varphi p=0$ | A = 0 |

Anode | $\u2202np=0$ | u = 0 | $\u2212kh\u2202nTh=hw(Th\u2212Tw)$ | $\u2202nTe=0$ | $\varphi p=0$ | $\u2202nA=0$ |

Cathode | $\u2202np=0$ | u = 0 | T = _{h}T _{c} | $\u2202nTe=0$ | $\u2212\sigma \u2202n\varphi p=Jqcath$ | $\u2202nA=0$ |

Wall | $\u2202np=0$ | u = 0 | $\u2212kh\u2202nTh=hw(Th\u2212Tw)$ | $\u2202nTe=0$ | $\u2202n\varphi p=0$ | $\u2202nA=0$ |

Outflow_z | $p=p\u221e$ | $\u2202nu=0$ | $\u2202nTh=0$ | $\u2202nTe=0$ | $\u2202n\varphi p=0$ | $\u2202nA=0$ |

Outflow_y | $\u2202np=0$ | $\u2202nu=0$ | $\u2202nTh=0$ | $\u2202nTe=0$ | $\u2202n\varphi p=0$ | $\u2202nA=0$ |

Boundary . | p
. | u
. | T
. _{h} | T
. _{e} | $\varphi p$ . | A
. |
---|---|---|---|---|---|---|

Inflow | $\u2202np=0$ | u = u _{i} | T = _{h}T_{0} | T = _{e}T_{0} | $\u2202n\varphi p=0$ | A = 0 |

Anode | $\u2202np=0$ | u = 0 | $\u2212kh\u2202nTh=hw(Th\u2212Tw)$ | $\u2202nTe=0$ | $\varphi p=0$ | $\u2202nA=0$ |

Cathode | $\u2202np=0$ | u = 0 | T = _{h}T _{c} | $\u2202nTe=0$ | $\u2212\sigma \u2202n\varphi p=Jqcath$ | $\u2202nA=0$ |

Wall | $\u2202np=0$ | u = 0 | $\u2212kh\u2202nTh=hw(Th\u2212Tw)$ | $\u2202nTe=0$ | $\u2202n\varphi p=0$ | $\u2202nA=0$ |

Outflow_z | $p=p\u221e$ | $\u2202nu=0$ | $\u2202nTh=0$ | $\u2202nTe=0$ | $\u2202n\varphi p=0$ | $\u2202nA=0$ |

Outflow_y | $\u2202np=0$ | $\u2202nu=0$ | $\u2202nTh=0$ | $\u2202nTe=0$ | $\u2202n\varphi p=0$ | $\u2202nA=0$ |

A heavy-species temperature *T _{h}* =

*T*

_{0}= 500 K is imposed at the inlet. A parabolic heavy-species temperature distribution is specified on the cathode surface $Tc=Tc0(1\u2212(y/ycath)2\u2212((z\u2212zoff)/zcath)2)$, where the maximum cathode temperature

*T*

_{c}_{0}= 3000 K (approximately equal to the melting point for Tungsten

^{28}) and the minimum value of

*T*occurs along the cathode boundaries with a value close to the wall-cooling reference temperature

_{c}*T*= 1000 K;

_{w}*y*and

*z*are the spatial coordinates,

*y*and

_{cath}*z*are the characteristic electrode lengths (

_{cath}*y*=

_{cath}*z*=

_{cath}*d*/2 =

_{x}*d*/2 = 1.5 mm), and

_{y}*z*is the offset distance from the origin to the center of the cathode. The anode and wall boundaries are described as experiencing forced water cooling over a metal surface with a convective heat transfer coefficient of Refs. 31 and 35,

_{off}*h*= 2 × 10

_{w}^{4}W m

^{−2}K

^{−1}, and reference wall temperature

*T*= 500 K. Such model is a rough approximation of the heat transfer from the plasma through the walls, yet an approximation extensively used in thermal plasma flow simulations.

_{w}^{27}A more adequate model involves including part of the solid wall within the computational domain, as done in Refs. 52 and 53. Finally, the electron temperature,

*T*is imposed with zero normal gradient condition over all the boundaries except at the inflow, where the inflow gas is assumed to be in LTE, i.e.,

_{e}*T*=

_{e}*T*=

_{h}*T*

_{0}.

^{44}

The current density *J _{qcath}* profile follows a quasi-Gaussian distribution over the cathode boundary, i.e., $Jqcath=Jqcath0\u2009exp\u2009(\u2212(r/Rc)nc),$ where

*r*= $(x2+y2)12$ is the radial coordinate;

*n*,

_{c}*R*and

_{c}*J*

_{qcath}_{0}are applied such that the total current $Itot=\u222bJqcathdS$, and

*S*represents the area of the cathode surface.

## III. CONTROLLING PARAMETERS

The main parameters of the arc in crossflow are the total current *I _{tot}*, the inflow velocity

*U*, the inter-electrode spacing

_{imax}*H*, and the type of gas [Fig. 1(a)]. The effects of the total current and the inflow velocity on the arc in crossflow have been studied experimentally,

^{45}as well as computationally.

^{32,33}A computational characterization of the arc in crossflow was presented in Ref. 11 using two non-dimensional numbers: the Reynolds number

*Re*and the Enthalpy number $\Pi h$. A similar approach is pursued in this study. The Reynolds number

*Re*characterizes the ratio of the magnitudes of advective transport to diffusive transport. In contrast to,

^{11}the definition of the Enthalpy number used in the present work $\Pi I$ characterizes the ratio of the magnitudes of electrical energy to thermal energy transport (i.e., $\Pi I=\Pi h\u22121$). The Reynolds number can be understood as characterizing the

*strength of the flow*, i.e., the stronger the flow, the larger will be the

*Re*. Similarly, the enthalpy number characterizes the

*strength of the arc*, i.e., the stronger the arc, the larger will be the $\Pi I$.

The Reynolds and the Enthalpy numbers for the present study are defined by

and

where the subscript “*r*” denotes the property evaluated at a reference discharge temperature and “*i*” denotes values evaluated at the inlet conditions. In the present work, the reference temperature is set to 16.7 kK, which corresponds to a recommended representative temperature for an argon plasma.^{15}

## IV. RESULTS AND DISCUSSION

### A. Solution fields

Representative three-dimensional (3D) steady-state solution fields are presented in Fig. 2 for two representative conditions: (left column) strong-flow–weak-arc, corresponding to *Re *=* *8630 and $\Pi I$ = 85.4, and (right column) weak-flow–strong-arc, corresponding to *Re *=* *3000 and $\Pi I$ = 800.

The resulting arc shape is the result of the imbalance between the electromagnetic force emerging due to the electromagnetic field distribution and the drag force caused by the gas flow. Experimental observations of emission intensity^{28} are often used to characterize the arc shape. In that regards, the value of an effective excitation temperature, responsible for most optical emission, can be expected to lie in between the values of *T _{e}* and

*T*. Therefore, the actual arc shape can be expected to be somewhat in between the distributions of the heavy-species temperature

_{h}*T*[Fig. 2(a)] and the electron temperature distribution

_{h}*T*[Fig. 2(b)]. For the strong-flow–weak-arc configuration (Fig. 2, left), the arc bends more profoundly at the arc–gas interface, leading it to become bow-shaped, whereas for the weak-flow–strong-arc configuration (Fig. 2, right), the arc becomes cusp-shaped. The arc–gas interface is characterized by a marked increase in the degree of thermodynamic nonequilibrium, characterized by the electron-to-heavy-species temperature ratio

_{e}*θ*=

*T*/

_{e}*T*. The distribution of the thermodynamic nonequilibrium parameter

_{h}*θ*and the effects of the controlling parameters on the arc shape in the arc in crossflow plasma were discussed in Ref. 11.

The results in Fig. 2 show that the electric potential distribution is more diffusive for the weak-flow–strong-arc configuration, while the strong-flow–weak-arc configuration shows an abrupt change. The abrupt change of the electric potential occurs at the region where the arc bends due to strong drag force and is characterized with greater magnitude of the electric field leading to a dramatic directional change in the current path (for very high gas flow rates, the bending angle can be close to 180°). It is interesting to note that the electric potential distribution is a maximum at the electrode center for the weak-flow–strong-arc configuration, whereas the maximum shifts upstream with increasing the gas flow rate. This behavior may be the result of the strong influx of electrons from the plasma core to the plasma–cathode–gas interface,^{29} which leads to a high degree of thermodynamic nonequilibrium.

To better capture the plasma cathode interaction, the grid size should be small enough to resolve the strong gradients of the solution fields (e.g., heavy-species temperature, velocity). Moreover, the plasma flow model should be complemented with the sheath model to better describe charge transport phenomena near the electrodes. A limitation of the arc in crossflow model used in the present study is that the model cannot appropriately describe near-electrode phenomena (i.e., electrode sheaths). The coupling of NLTE-plasma flow–plasma-sheath model is an active area of research,^{27,46,47} which is beyond the aim of the present work.

The velocity magnitude ǁ**u**ǁ depends strongly on the imposed electric current and the inflow gas velocity. The imposed current causes strong Joule heating and the gas inflow causes the cooling of the arc; their combined effect leads to a drastic increase in local velocity, as shown in Fig. 2(d). The velocity magnitude ǁ**u**ǁ distribution shows distinct characteristics in the electrode upstream and in the electrode downstream regions, i.e., the gradient of velocity magnitude is more intense downstream of the electrode than upstream. For high-power arcs (e.g., high current, *I _{tot}* > 200 A, and flow rates,

*U*> 100 m s

_{i}^{−1}), the gradient extends further downstream, causing the attachment to extend along the wall. The velocity magnitude in NLTE model simulations is significantly higher than those obtained in LTE simulations. For example, Kelkar and Heberlein reported maximum velocities of up to 1300 ms

^{−1}using a LTE model,

^{32}whereas the NLTE model used here leads to maximum velocities of ∼1560 ms

^{−1}for a coarse grid (∼85 000 nodes)

^{11}and velocities of ∼1800 ms

^{−1}with the finer grid used in the present study. The large discrepancies appear to be the result to the larger heavy-species temperatures obtained with NLTE models compared to those from LTE ones,

^{11,31,36}as well as the greater resolution of heavy-species temperature gradients, particularly near the cathode region. Higher heavy-species temperature

*T*leads to lower mass density

_{h}*ρ*, which induces larger velocities in order to maintain mass flow (

*ρ*

**u**) continuity. The Mach number $Ma=\u2009||u||/c=\u2009||u||/(\u2202p/\u2202\rho )1/2$, where

*c*is the speed of sound, indicates that the inflow conditions lead to a largely sub-sonic flow (i.e.,

*Ma*≪ 1). Despite the large acceleration of the flow in the arc core (maximum ǁ

_{i}**u**ǁ ∼1800 ms

^{−1}), the local

*Ma*remains less than one everywhere given that

*c*∼ 10

^{3}–10

^{4}ms

^{−1}within the arc.

^{35}

As seen in Fig. 2(e), the self-induced magnetic field increases at the arc–gas flow interface and reduces to a minimum in the plasma core. The magnitudes of the obtained magnetic fields are comparable to those reported in the computational simulation of circuit breakers, i.e., ∼60 mT in the simulations here compared to ∼45 mT reported in Ref. 9. The Lorentz force distribution is similar to the magnetic field distribution; nevertheless, the former presents a relatively larger gradient near the electrodes, particularly near the cathode, due to the high values and constriction of current density. The magnitude of Lorentz force is on the order of 10^{6} N m^{−3}, which is significantly larger than the values of ∼10^{4} N m^{−3} reported in the simulation of an inductively coupled plasma (ICP) torch in Ref. 48 (torch diameter of 8 mm). The magnitude of the Lorentz force decreases drastically in the region near the mid-plane (mid inter-electrode region) as expected due to the weaker magnetic field and more diffuse distribution of current density.

### B. Arc dynamics

There are significantly fewer investigations of the dynamics in the arc in crossflow compared to those found in other arc configurations, such as non-transferred plasma torches.^{10,21,22} Non-transferred arc plasma torches experience complex plasma–gas interactions that often lead to the azimuthal rotation and continuously elongation of the arc, its subsequent extinguishment, and re-establishment at a different, near-cathode, location. In contrast, the arc in crossflow is a canonical configuration that involves perpendicular plasma–gas interaction leading to a relatively constant arc, a characteristic that makes this configuration ideal for the study of arc–gas glow interactions. Nevertheless, the obtained simulation results show that the arc tends to develop temporal fluctuations depending on the values of the controlling parameter *Re* and Π_{I}. These fluctuations are the result of (i) the continuous competition between electromagnetic and flow drag forces, (ii) the occurrence of thermal nonequilibrium, and (iii) the large variation of plasma material properties.

The arc in crossflow simulations are run for an extent of time that is a function of the intrinsic characteristic time *t _{char}* =

*L*/

_{z}*U*. This intrinsic time can be understood as a characteristic transit time for a particle to navigate the whole axial extent of the domain. The dynamic characteristics of the flow are then represented as a function of the dimensionless time

_{imax}*t*/

*t*

_{char}. To discern representative dynamic behavior of the flow, four representative spatial observation points are defined, as indicated in Fig. 3. Figure 3(a) shows the placement of the observation points: one near the cathode (near-cathode) and three along the center of the channel (i.e., upstream, midstream, and downstream). The near-cathode location is 0.3 mm away from the center of the cathode [Fig. 3(b)]. The upstream observation point is placed at the midpoint in between the electrodes (center of the channel); and the midstream and the downstream points are placed 10 and 20 mm, respectively, downstream from it along the channel's center axis. To make the observation points data grid-independent (i.e., not function of the size of the discretization grid), as well as to make the results more consistent with experimental measurements (e.g., via probes), each observation point consists of a cubic element of side 0.2 mm, as shown in Fig. 3(b). The reported data from a given observation point correspond to the volume-average of the distribution of the corresponding data over the corresponding cube.

To identify the flow dynamics, fluctuation signals are extracted from the observation points data. Representative fluctuation signals are presented in Fig. 3(c) for the heavy-species temperature *T _{h}* and for the conditions

*Re*=

*8630 and $\Pi I$ = 85.4, i.e., the strong-flow–weak-arc studied in Fig. 2 (left). The fluctuation signal is given by*

*δΤ*= (

_{h}*Τ*− $T\xafh$Z)/$T\xafh$, where $T\xafh$ is the time-averaged value at steady-state. As observed in Fig. 3(c),

_{h}*T*in the upstream point varies by over 15% with respect to its mean, whereas at the near-cathode, midstream, and downstream locations, it varies by less than 5%. The drastic change in

_{h}*δΤ*along the electrode axis might indicate marked dynamics at the plasma-arc–gas-flow interaction region.

_{h}In addition to their relative magnitude, the fluctuation signals are characterized in terms of their frequency spectra. Figure 3(d) shows the frequency-domain counterpart to the temporal signals in Fig. 3(c). The power spectrum is obtained by normalizing the product of the Fourier domain data and its complex conjugate with the sample size. From Fig. 3(d), four main peaks are observed for the upstream condition, i.e., *f _{p}*

_{1}= 9.5,

*f*

_{p}_{2}= 14.3,

*f*

_{p}_{3}= 19.1, and

*f*

_{p}_{4}= 26.2 kHz; and their corresponding power (amplitude) is:

*P*

_{1}= 4.0,

*P*

_{2}= 8.8,

*P*

_{3}= 1.0, and

*P*

_{4}= 55.3 K

^{2}Hz

^{−1}, respectively. The low frequency peaks are usually due to large-scale variations, which are mostly closer to the fundamental frequency of the system. The higher frequency peaks correspond to small-scale property variations. The 4th peak frequency, i.e.,

*f*

_{p}_{4}= 26.2 kHz, shows that the frequency corresponding to the small-scale has the highest magnitude of the power spectrum (

*P*

_{4}= 500 K

^{2}Hz

^{−1}) at the near-cathode region. The small-scale variation in

*T*provides some notion into the sensitivity of the cathode attachment; particularly addressing the time scale and the intensity at which the electromagnetic and the drag forces continuously balance each other. Thus, to effectively analyze the plasma flow dynamics, adequate resolution of small-scale signals, as well as appropriate handling of the coupling among large- and small-scales, is necessary. The use of the non-linear VMS (VMS

_{h}_{n}) method aims to address these requirements, as described next.

### C. Role of the numerical solution approach

Numerical methods for thermal plasma flow simulations need to be *robust* in order to deal with the marked non-linearity of the associate models, as well as *accurate* in order to resolve the multi-scale characteristics of the solution fields, particularly if the aim is to describe dynamic and unstable behavior. The VMS method, although it has been satisfactorily applied to various plasma flow simulations,^{12,27,31,35,37,49} it is not inherently suitable for the description of multiscale phenomena in which the small- and large- scales interact in a highly coupled and nonlinear manner as it is the case in turbulent and instability phenomena. The non-linear VMS method presented in Sec. II B (i.e., VMS_{n}) addresses up-front the non-linearity of the coupling between small- and large-scales. The VMS_{n} method has been successfully applied to representative incompressible and compressible laminar and turbulent flow simulations and to the simulation of the flow from non-transferred arc plasma torches.^{42} Therefore, the VMS_{n} method is also expected to provide a more accurate description of the dynamic behavior in the arc in crossflow than the VMS method.

Figure 4 presents a representative comparison between the results obtained with the VMS and VMS_{n} methods for the same conditions used in Fig. 2(left), strong-flow–weak-arc (*Re *=* *8630 and $\Pi I$ = 85.4). The simulations were run until a (fluctuating) steady-state is achieved in all solution fields, i.e., moving time-averaged quantities ($T\xafh$, $T\xafe$, $\varphi \xafp$, etc.) remain constant in time. Temporal fluctuation signals and power spectra for electric potential are shown in Fig. 4(a), for heavy-species temperature in Fig. 4(b) and for electron temperature in Fig. 4(c) for each of the four observation points (i.e., near-cathode, upstream, midstream, and downstream).

The comparison shows the relative effectiveness of the VMS and VMS_{n} methods for capturing dynamic behavior. Specifically, the results in Fig. 4 show that the fluctuation signals obtained with the VMS_{n} method are more uniform (periodic) than those obtained with VMS. The electric potential fluctuations presented in Fig. 4(a) show the more uniform periodicity of the VMS_{n} results with respect to those from VMS, which are quantified by the different peaks in the frequency-domain signals. The VMS method results show five frequency peaks, i.e., *f _{p}*

_{1}= 5.6,

*f*

_{p}_{2}= 9.4,

*f*

_{p}_{3}= 13.2,

*f*

_{p}_{4}= 18.9, and

*f*

_{p}_{5}= 22.6 kHz, with power

*P*

_{1}= 3.5 × 10

^{−4},

*P*

_{2}= 8.7 × 10

^{−5},

*P*

_{3}= 2.1 × 10

^{−5},

*P*

_{4}= 1.1 × 10

^{−3}, and

*P*

_{5}= 2.2 × 10

^{−3}V

^{2}Hz

^{−1}, respectively for the near-cathode region. In contrast, the VMS

_{n}method shows only four peaks (named from 2nd to 5th for consistency with the VMS results), which are closely comparable to the peaks obtained from VMS method. The power of the 2nd, 3rd, and 4th peaks obtained using the VMS

_{n}method are remarkably low compared to those obtained with VMS. The 3rd frequency peak (

*P*

_{3}) of electric potential has lesser power (both in VMS and VMS

_{n}) compared to the other peaks (2nd, 4th, etc.). The results show that the 5th frequency peak (

*P*

_{5}) obtained using the VMS method is comparable to the 5th peak obtained using the VMS

_{n}method. The significant magnitudes of the 1st and 4th frequency peaks obtained with the VMS method are not seen in VMS

_{n}results.

The spectra of the electric potential fluctuations [Fig. 4(a)] are comparable to those for the heavy-species temperature [Fig. 4(b)] and the electron temperature [Fig. 4(c)]. The time-domain fluctuations of *T _{h}* obtained with the VMS method [Fig. 4(b), left] are quasi-periodic, whereas those obtained with VMS

_{n}[Fig. 4(b), right] are smooth and periodic. Although the peak frequencies of the fluctuations are approximately the same for

*ϕ*

_{p},

*T*, and

_{h}*T*, their corresponding power is significantly different. This is expected given that the magnitude of

_{e}*ϕ*

_{p}is ∼

*O*(10

^{1}), whereas for

*T*and

_{h}*T*it is ∼

_{e}*O*(10

^{3}–10

^{4}). Comparing the power spectra in Fig. 4(b), one can quantitatively assess the difference between the results from the VMS and the VMS

_{n}methods. The 2nd peak (

*P*

_{2}) and the 3rd peak (

*P*

_{3}) are comparable in magnitudes between the two methods (VMS and VMS

_{n}), i.e., for the near-cathode region,

*P*

_{2}= 1.5 and 4.8 K

^{2}Hz

^{−1}and

*P*

_{3}= 12.6 and 9.0 K

^{2}Hz

^{−1}for VMS and VMS

_{n}, whereas the 5th peak (

*P*

_{5}) is remarkably different, i.e.,

*P*

_{5}= 87.0 and 491.4 K

^{2}Hz

^{−1}obtained using VMS and VMS

_{n}, respectively. This marked deviation in the peak frequency can be analyzed considering that the low (fundamental) frequency peak corresponds to the large-scale variations and the high frequency peaks correspond to the small-scale variations; hence, the sequence

*f*

_{p}_{1}>

*f*

_{p}_{2}>

*f*

_{p}_{3}>

*f*

_{p}_{4}>

*f*

_{p}_{5}can be correlated with the resolution from large-scales to progressively smaller-scales. The relative agreement between

*P*

_{2}

_{VMS}≈

*P*

_{2}

_{VMSn}and

*P*

_{3}

_{VMS}≈

*P*

_{3}

_{VMSn}indicates that VMS

_{n}can capture all the large-scale variations similarly as VMS, whereas the observed

*P*

_{5}

_{VMS}≪

*P*

_{5}

_{VMSn}indicates that VMS

_{n}captures the small-scales significantly more accurately than VMS. Given the large power of the frequency signals at the arc-gas interface (i.e., both the near-cathode and the upstream locations), it can be concluded that the VMS

_{n}provides a better description of arc–gas flow interactions than VMS. This finding motivates the use of the VMS

_{n}method for the identification of flow regimes.

## V. ARC CHARACTERISTICS AND REGIMES OF THE ARC IN CROSSFLOW

### A. Arc shape and characteristics

The arc shape and the arc behavior are the result of the energy balance between the supplied electrical energy (Joule heating) and the energy lost as a result of convective cooling to the stream of working gas. The phenomenological characterization of the arc in crossflow with respect to the arc shape was presented in Ref. 11 through the computational simulations of 19 cases for the range of Reynolds and the Enthalpy dimensionless numbers, 0 < *Re *<* *10 000 and 2 < $\Pi I$ < 350. As part of the present research, the arc shape (Fig. 5) and the arc characteristics (Fig. 6) were obtained for an extended *Re*-$\Pi I$ range (i.e., 0 < *Re *<* *10 000 and 50 < $\Pi I$ < 800).

Consistent with experimental and numerical studies of the arc in crossflow,^{11,29,45} the arc shape changes from cusp- to bow-shaped for increasing *Re* and decreasing $\Pi I$, resulting in increased arc length (e.g., Fig. 5: case B_{1} to case D_{1} or B_{4} to case D_{1}) and increased electric potential [Fig. 6(a)]. For increasing *Re* and $\Pi I$ (Fig. 5: case B_{1} to case D_{4}), the arc temperature [Fig. 6(b)] and arc power [Fig. 6(c)] increase due to strong electromagnetic pumping. The overall increase in arc temperature leads to large temperature gradients as well as increased local acceleration across the electrodes. The variation of *Re* and $\Pi I$ among the cases is accommodated by modifying the inlet velocity *U _{i}*

_{max}(i.e., from 1.2 to 117.0 m s

^{−1}) and the imposed current

*I*

_{tot}(i.e., from 16 to 657 A).

The time step size for all the cases (B1–D4) is on the order *O*(10^{−6}) to *O*(10^{−7}), similar to our previous works,^{11,17} except for the high Reynolds and low Enthalpy number cases (e.g., D1) for which the time step is on the order of *O*(10^{−8}). The smaller time step size was necessary due to the large elongation and the highly dynamic nature of the arc. Additionally, the case D1 involved a greater number of time steps to obtain a converged solution.

The controlling parameters *Re* and $\Pi I$ determine not only the overall arc shape but also the arc behavior (i.e., the characteristics of property fluctuations). The effects of the Reynolds and the Enthalpy numbers on the arc behavior are analyzed next.

### B. Effect of Reynolds and Enthalpy numbers on the arc dynamics

The characteristics of the arc as a function of the Reynolds number are depicted in Fig. 7 and those as obtained as a function of the Enthalpy number in Fig. 8. Figure 7 (left) shows that the arc evolves toward having an anode attachment further downstream with increasing *Re*. The magnitude and frequency of the temporal fluctuations of *T _{h}* also increase with increasing

*Re*, as depicted in Fig. 7 (center). The fluctuations change from steady [Fig. 7(a)], to quasi-periodic [Fig. 7(b)], and then to near-periodic with higher amplitudes [Fig. 7(c)]. The increased characteristic fluctuations with increasing

*Re*results in an increase in the number of peak frequencies

^{21}and an increase in the magnitude of their power, as seen in right of Fig. 7 (right). Low

*Re*simulations are characterized by the presence of a single frequency peak due to negligible variation temporal fluctuations [i.e., in Fig. 7(a), for

*Re*=

*3000,*

*δΤ*< 1% and

_{h}*f*= 1.3 kHz], and the amplitude fluctuations and the number of peak frequencies increase with

_{p1a}*Re*(i.e., for

*Re*=

*7500 [Fig. 7(b)],*

*δΤ*≈ ±13% and presents two dominant frequencies:

_{h}*f*= 11.1 and

_{p2a}*f*= 20.6 kHz; and for

_{p2a}*Re*=

*10 000 [Fig. 7(c)],*

*δΤ*≈ ±25% and presents four dominant frequencies:

_{h}*f*= 13.2,

_{p3a}*f*= 26.4,

_{p3b}*f*= 37.7, and

_{p3c}*f*= 49.0 kHz. The increase in the number of dominant frequencies indicates the increasing complexity of the temporal fluctuations and, hence, of characteristic temporal scales (each scale given by the inverse of a characteristic frequency), as well as spatial (small-) scales (associated with the amplitude of the fluctuations).

_{p3d}Additionally, at low Reynolds number [i.e., *Re *=* *3000, Fig. 7(a), right], the power of the fluctuation signals in the near-cathode region is higher than at any other locations (e.g., *P*_{1} = 29.3 and 16.2 K^{2} Hz^{−1} at the near-cathode and upstream regions, respectively). However, for increasing Reynolds, the power of the fluctuation signals increases at the Upstream observation point and surpasses the fluctuations at the near-cathode at higher Reynolds (i.e., *P*_{2} = 6.8 × 10^{4} and 6.6 × 10^{4} K^{2} Hz^{−1}, *P*_{3} = 7.8 × 10^{4} and 10.2 × 10^{4} K^{2} Hz^{−1}; at the near-cathode and upstream regions, respectively). The strong fluctuation at the Upstream observation-location is due to the downstream arc anode attachment and the strong nonequilibrium at the plasma–gas interface.

The results in Fig. 8 show that with increasing Enthalpy number, the anode attachment moves upstream, the arc shape changes from bow- to cusp-shaped (i.e., the arc gets bulgy), and the amplitude and frequency of fluctuations decrease [i.e., Figs. 8(a)–8(c), *δΤ _{h}* reduces from ±0.5% to an insignificant value], and hence the arc becomes relatively more stable. The relative arc stability can be assessed by observing the frequency spectra, where the power of the signals gets reduced drastically, by over four orders of magnitudes, with increasing $\Pi I$ (i.e.,

*P*

_{1}= 29.3, 16.2, 18.9 K

^{2}Hz

^{−1};

*P*

_{2}= 0.4, 1.6, 1.4 K

^{2}Hz

^{−1}; and

*P*

_{3}= 1.6 × 10

^{−4}, 8.6 × 10

^{−3}, 4.9 × 10

^{−3}K

^{2}Hz

^{−1}at the near-cathode, upstream, and downstream observation points, respectively; the subscripts 1, 2, and 3 correspond to $\Pi I$ = 200, 350, and 800, respectively). The drastic decrease in the power of the signals indicates the establishment of an increasingly stable configuration of the plasma, and the signals can be considered as

*white noise*.

At higher Reynolds, increasing the Enthalpy number not only increases arc stability but also increases the maximum arc temperature [Fig. 6(b)] and the overall arc power [Fig. 6(c)]. The increase in the arc temperature causes large temperature gradients, specifically in the region near the cathode (e.g., see Fig. 5, case D4). To effectively study arc dynamics at high *Re* and $\Pi I$, the grid resolution should be fine enough to appropriately resolve the associated large gradients present in the solution fields. The use of a finer discretization in the present study compared to the one used in our previous study^{11} has significantly improved the convergence of the numerical simulations. Nevertheless, the use of the finer grid resulted in a significant increase in the computational cost (∼1.5–2 times greater than in Ref. 11). The present study indicates that the magnitude of the peak frequencies (*f _{p}*) and the increase in the number of peak frequencies with Reynolds number (e.g., from 1 for Re = 3000 to 4 for Re = 10 000, as seen in Fig. 7) for the studied

*Re*-$\Pi I$ range is grid independent and that even finer discretizations are required to analyze the arc in crossflow for a wider

*Re*-$\Pi I$ range.

### C. Spatiotemporal fluctuations and regimes identification

The dynamics of the arc in cross flow as function of the Reynolds and Enthalpy numbers can be summarized using three quantities: the time-averaged values, the peak frequency of fluctuations, and the variance of those fluctuations. These quantities are presented in Figs. 9 and 10 for *ϕ _{p}*,

*T*, and

_{h}*T*as representative flow properties, and for each observation point.

_{e}Figure 9 presents the variation of the time-averaged electric potential $\varphi \xafp$, heavy-species temperature $T\xafh$, and electron temperature $T\xafe$ as a function of *Re* and $\Pi I$ at each observation point (i.e., near-cathode, the upstream, the midstream, and the downstream). These quantities are defined using a moving-average procedure, i.e., $\varphi \xafp$ for a given observation point is defined by

where **x**_{op} represents the spatial location of the given observation-point, *t _{f}* is the final time of the simulation (

*t*≫

_{f}*t*), and Δ

_{char}*t*is an arbitrary time interval used for the analysis such that Δ

_{f}*t*>

_{f}*t*. The expression in Eq. (10) is evaluated for a temporal extent that ensures that steady-state has been achieved such that $d\varphi \xafp/dt=0$ independently of the specific value of Δ

_{char}*t*. Similar expressions are used to define $Th\xaf$ and $Te\xaf$.

_{f}In Fig. 9, it can be observed that for high *Re* and low $\Pi I$ (i.e., near case D_{1}, see Fig. 5) at the near-cathode location, $\varphi p\xaf$ is maximum (i.e., 33.8 V), as expected due to the large arc length [consistent with the maximum electric potential distribution in Fig. 6(a)]. This implies that the high energy loss by advective cooling has to be compensated by greater electrical energy consumption. Also, at this point $Th\xaf$ is minimum (i.e., 12.2 kK) due to the strong cooling, whereas $Te\xaf$ remains largely unaffected with varying *Re* for the same value of $\Pi I$ (i.e., 16.4 kK).

For high *Re* and $\Pi I$, also at the near-cathode location, $T\xafe$ increases due to strong Joule heating (i.e., ∼37 kK for case D_{4} in Fig. 5). At the upstream location, the values of $\varphi p\xaf$ and $Th\xaf$ are similar to those at the near-cathode location, whereas the $Te\xaf$ distribution remains unaffected for increasing *Re*. At the midstream and the downstream observation points, the distributions of $\varphi p\xaf$, $Th\xaf$, and $Te\xaf$ are somewhat analogous to their corresponding instantaneous maxima (depicted in Fig. 6). The distributions of the maximum heavy-species temperature *T _{h}*

_{max}[Fig. 6(b)] and the average heavy-species temperature $Th\xaf$ at the midstream and the downstream (Fig. 9) are analogous to the power distribution [Fig. 6(c)], which indicates that the increase in arc power increases the overall temperature. In contrast, the heavy-species in the near-cathode and the upstream locations undergo drastic reductions in temperature due to the strong interaction with the gas flow stream. The reduction in

*T*is a result of the strong cooling of the arc by the gas flow, and the increased degree of nonequilibrium, which can potentially lead to turbulence.

_{h}^{13}

The phenomenon of abrupt drop in the heavy-species temperature with increasing *Re* and for low values of $\Pi I$ (Fig. 9, upstream and near-cathode observation points) is somewhat analogous to the drastic reduction in the plasma jet length in a non-transferred plasma torch due to the increase in gas flow rate, which results in the transition from laminar to turbulent flow.^{22,50,51} This behavior can be better understood by inspecting the frequency and variance of fluctuations.

Figure 10 presents the variation of the peak frequency *f*_{max} of the fluctuations electric potential, heavy-species, and electron temperatures as a function of *Re* and $\Pi I$ at each observation point. The peak frequencies correspond to the frequencies with the highest power, e.g., as shown in Figs. 7 and 8.

In the near-cathode region, for increasing *Re*, the peak frequencies increase by ∼1–2 orders of magnitude, whereas increasing $\Pi I$ at constant *Re*, the relative change in peak frequency is subtle. This is somewhat consistent with observations in plasma torches,^{19} in which the peak frequency of voltage fluctuations increases mostly with increasing gas flow and varies little with increasing total current. The peak frequencies of *T _{h}* and

*T*in the upstream location are similar to those at the near-cathode. These peak frequencies increase when moving downstream (i.e., from the upstream to the midstream, and then to the downstream observation points). In contrast, no significant change is observed for

_{e}*ϕ*among the different observation points.

_{p}The observed trends for the peak frequency *f _{max}* of the fluctuations as a function of

*Re*and $\Pi I$ have to be contrasted against their amplitude. In this regard, Fig. 11 presents the distribution of the variance

*s*

^{2}of the fluctuations of electric potential, heavy-species temperature, and electron temperature as a function of

*Re*and $\Pi I$ at each observation point.

The variance quantifies the squares of the amplitudes of the fluctuation signals. Specifically, the variance of the electric potential $\varphi p$ is given as

Similar expressions are used for the variances of heavy-species and electron temperatures.

It is to be noted that the variance plots in Fig. 11 are significantly non-smooth, showing distinct faceted regions. Non-smoothness of the distributions of peak frequency with *Re* and $\Pi I$ can also be observed in the results in Fig. 10, but to a lower degree. The lack of smoothness in *σ*^{2} is the results of the very limited amount of points used to probe the *Re*-$\Pi I$ map (i.e., 12 points, as indicated in Fig. 5). The small number of *Re*-$\Pi I$ sets studied is due to the large computational cost of the simulations, each involving the fully-coupled solution of the 10-variable equations-set given by Eq. (1), for a total of ∼2.5 × 10^{9} unknowns and for over ∼6000 time steps (∼0.1–3.0 *μ*s) after achieving steady-state, and requiring approximately 5 weeks of computing time per case using 16 Intel E5–2600 cores.

The magnitudes of $s\varphi 2$ compared to those of $sTh2$ and $sTe2$, as shown in Fig. 11, indicate that the fluctuations in electric potential are significantly smaller than those for the heavy-species and electron temperatures [i.e., the former vary from *O*(10^{−4}) to *O*(10^{2}), whereas the latter from *O*(10^{1}) to *O*(10^{6})]. This can be expected given that temperatures distributions are more dependent on the plasma flow (e.g., heat generated by Joule heating vs heat loss by convection) than the electric potential distribution (which depends mostly on the electrodes configuration and total imposed current). Additionally, of $sTh2$ and $sTe2$ rapidly increase with increasing *Re* and $\Pi I$ at the near-cathode and upstream observation points, whereas the slowly increase with increasing *Re* and $\Pi I$ at the midstream and the downstream locations.

The distributions of peak frequency (Fig. 10) and variance (Fig. 11) of the heavy-species temperature fluctuations at the midstream and downstream observation points suggest they strongly correlate with the arc dynamics and the degree of plasma–gas flow interaction. Between these two locations, the midstream observation point remains more consistently within the bulk plasma, and it is therefore preferred for the assessment of flow characteristics. Therefore, *f _{Th}*

_{max}and of $sTh2$ can be used to assist in determining characteristic regimes of the arc in crossflow. These quantities are employed to determine the regimes of the arc in crossflow, as shown in Fig. 12.

Figure 12 shows the distribution of the product $fThmaxsTh$, where $sTh=(sTh2)12$ is the standard deviation, at the midstream observation point, which is used as a metric for regime identification. The results in Fig. 12 show the regions delineating four distinct flow regimes in the *Re*-$\Pi I$ map, i.e.,: I—static, II—periodic, III—quasi-periodic, and IV—chaotic. These regions approximately present different characteristics in terms of the frequency of fluctuations and their magnitude: region I—static presents the low-frequency fluctuations and of very small magnitude; region II—periodic presents the well-defined fluctuations of higher frequency and amplitude; region III—quasi-periodic presents more fluctuations, with several peaks in the frequency-domain, and increased amplitude of fluctuations; and finally, region IV—chaotic shows more complex frequency-domain signals and of greater amplitude. The transition from a laminar to a (potentially) turbulent regime of the arc in crossflow occurs from region II to region III (can be termed as the laminar–turbulent transition region) or region IV; and is governed by both the Reynolds number and the Enthalpy number. It is to be noted that region IV comprises only one set of results (i.e., case D_{1} in Fig. 5). This simulation was significantly more computationally expensive compared to the others in terms of the large number of non-linear iterations per time step required for convergence and the large number of time-steps needed to achieve steady-state. Such behavior in a numerical simulation is often indicative of the onset of instability and potentially turbulent behavior,^{34} which requires significantly higher spatial and temporal resolution for their description (i.e., finer grid and smaller time steps). Nevertheless, for consistency, the simulation of case D_{1} used the same discretization grid as for all the other cases (i.e., 250k nodes).

The approach used to delineate the regions for each regime is similar to that used in Ref. 51, and consists on the determination of the separation lines **v**_{1}, **v**_{2}, and **v**_{3}. These lines express the Enthalpy number as a function of the Reynolds number, and each line is of the form

where *k* indicates the index of the line (*k *=* *1, 2, or 3), and *α* and *β* are coefficients defining the line,. The extent of each regime can be determined for a given set of Re and $\Pi I$ as

Even though the followed approach is at most approximate for the identification of distinct flow regimes (i.e., there are not unambiguously distinct regions in the distribution of $fThmaxsTh$ through the *Re*-$\Pi I$ map), a distinct advantage of it is that the approach could be adopted experimentally, i.e., by recording *in situ* spatial-temporal data of heavy-species temperature at a location downstream from the electrodes yet within the bulk arc plasma.

## VI. SUMMARY AND CONCLUSIONS

The flow dynamics and the establishment of different flow regimes in the arc in crossflow are studied using a 3D time-dependent NLTE plasma flow model. A nonlinear Variational Multiscale (VMS_{n}) numerical discretization approach is used to describe the coupling between the large- and small-scale features of the flow, providing greater numerical accuracy than traditional VMS methods. The Reynolds number and the Enthalpy number are used as control parameters. The dynamic behavior of the flow is characterized using the property variations gathered at four representative observation points: one near the cathode and three along the main axis for the flow. The time-averaged values of the temporal signals, as well as the frequency and variance of the signal fluctuations, are used to quantify the flow dynamics. The product of the standard deviation and the peak frequency of the heavy-species temperature at the midstream observation point is used to identify the distinct regimes, i.e., steady, periodic, quasi-periodic, and chaotic. The computational results reveal the role of increasing the relative arc strength in enhancing flow stability by delaying the growth of fluctuating and unstable flow behavior.

## ACKNOWLEDGMENTS

The authors gratefully acknowledge the support from the U.S. National Science Foundation through award CBET-1552037 and the U.S. Department of Energy through award DE-SC0018230.