Progress of computational plasma fluid mechanics

This article reviews and discusses the recent progresses of studies with the concept of “Computational plasma fluid mechanics.” Computational demonstrations show that the inhouse simulation codes such as PLasma All-Speed Turbulence with Implicit Pressure Code have captured hydrodynamic instabilities and reproduced flow dynamics in thermal plasma—nonionized gas coexisting systems. A unique method has made it feasible to study collective growth of binary alloy nanoparticles by numerical analysis. Smoothed Particle Hydrodynamics method with incompressibility modification has achieved complex behaviors of molten metal involving phase change, flow, heat transport, material mixing, and large deformation during arc welding. It is essential to study thermal plasma processes as comprehensive fluid systems in which hot plasma, cold nonionized gas, and materials coexist. The viewpoint and approaches of fluid mechanics as well as plasma physics are indispensable. Computational study will play a more important role in giving us new and deeper insights.


Introduction
Plasma, an ionized gas state, exhibits a variety of characteristics depending on the generating conditions. Low-pressure plasmas are usually understood on a molecular scale because the motions of each molecule, atom, ion, and electron play important roles individually. Meanwhile, high-pressure plasmas such as thermal plasmas and fusion plasmas are dominated by "fluid" dynamics.
One more keyword is "mechanics." Herein, understanding of only plasma is not enough. With intention of plasma engineering, it is essential to consider the physics and mechanisms in a comprehensive system where electrically conductive high-temperature plasma and nonconductive lowtemperature nonionized gas coexist. Their interface regions have remarkably large differences of temperature, density, and transport properties caused by dissociation and ionization. Furthermore, the system occasionally involves interaction with materials, which drives significant transports of mass, momentum, and heat. This is what makes a significant difference between "plasma fluid mechanics" and fusion plasma research progressed by the great efforts of many scientists in the long history.
From this perspective, "plasma fluid mechanics" mainly targets thermal plasma flow systems in which nonionized gas, liquid metal, or metal aerosol coexist. Thermal plasmas, such as arc plasma and RFICP (RF inductively coupled plasma), are generated at atmospheric pressure, where the energy input to electrons by the discharge is transferred to heavy particles (ions, neutral atoms and molecules) by HF collisions, and consequently their temperatures are almost the same. The temperature of the heavy particles is also higher than 10 000 K, resulting in a high enthalpy state. In other words, thermal plasma is an electromagnetic thermofluid that is generated and sustained through not only dissociation, ionization, and temperature increase, but also volume expansion with decrease in mass density.
The volume expansion is achieved by the increase in flow velocity. Strong velocity shear is generated at the interface region between the fast, low-density thermal plasma and the slow, high-density nonionized gas, which produces vortices due to the Kelvin-Helmholtz instability. Those multiple vortices interact with each other, deform, merge, and break up, eventually forming a turbulent field composed of multiscale vortices. Such vortices of thermal plasma-induced turbulence were actually captured in the experimental studies for a DC non-transferred arc plasma jet of argon gas by Pfender et al. 1) and a DC transferred arc plasma sustained in SF 6 (sulfur hexafluoride) gas by Inada et al. 2,3) Plasma spraying and nanoparticle mass-production are typical applications of thermal plasma. Because vortices convey thermal energy and substances, a turbulent field composed of vortices dominates the transport phenomena in the whole flow system. The mixing process of the turbulent vortices considerably affects the heat transfer from plasma to raw materials and the environment for nanoparticle formation.
In arc welding, even if the flow is less turbulent, vortices can also deteriorate welding quality by entrain oxygen into the shielding inert gas. 4,5) Conversely, vortices are generated in the weld pool and strongly affect the heat transfer, which determines the depth of weld penetration. 6,7) In addition to such vortices, the electromagnetic-thermal-hydrodynamic interaction between the arc plasma and the deforming weld pool or the melting electrode determine the welding process as well.
As described above, the importance of vortices and electromagnetic-thermal-hydrodynamic interactions in thermal plasma processes has been suggested. However, thermal plasma flow systems have high temperatures exceeding 10 000 K or even 20 000 K with strong light emission. Furthermore, the processes such as vortex generation, arc plasma quenching, nanoparticle formation, and electrode deformation take place in time scales from microseconds to milliseconds. Those conditions make direct observation and measurement arduous. The experimental approaches are faced with technological limitations, and the research and development still rely on the experience of engineers and researchers.
In order to overcome those difficulties, the computational analyzes based on the physical and mathematical theories are expected to advance the fundamental and comprehensive understandings of the complex nonlinear phenomena in thermal plasma flow systems.
This paper presents and reviews the recent progresses of the computational studies based on the concept of "plasma fluid mechanics" from the viewpoints of fundamental fluid dynamics of thermal plasma-nonionized gas coexisting systems and mechanical engineering for their applications.
2. Computational plasma fluid dynamics 2.1. Outline of thermal plasma-nonionized gas coexisting flow simulation 2.1.1. Mathematical descriptions. The flow field where an electrically conductive high-temperature thermal plasma and a nonconductive low-temperature gas interact with each other has great variations of the thermodynamic and transport properties such as density, viscosity, thermal conductivity, and specific heat at constant pressure. The entire flow field is expressed by electromagnetic thermofluid approximation under the fundamental assumptions: (i) the pressure is almost atmospheric, (ii) local temperatures of electrons, ions, and neutral atoms and molecules are identical, and (iii) the plasma is optically thin. Therein, the governing equations are given by the conservations of mass, momentum, chemical species, and energy as r u, p, , h j, B, and g respectively stand for time, position, mass density, velocity, pressure, viscosity, electric current density, magnetic flux density, and gravitational acceleration for Eqs. (1) and (2). In Eq. (3), Y , D, and m  denote mass fraction, diffusion coefficient, and a net rate of mass production, respectively. Subscript a means chemical species. Subscripts i, j, and k denote directions, abiding by Einstein summation convention. ij d and ijk e stand for Kronecker's delta and Eddington's epsilon, respectively. In Eq. (4), h, , l C , p k , B e, T , E, and R  denote enthalpy, thermal conductivity, specific heat at constant pressure, Boltzmann constant, elementary charge, temperature, electric field, and net radiation rate, respectively. S ij is the velocity strain tensor defined as (3), and (4) holds for compressible flows as well, although they are often described as The derivation is described in Because thermal plasma is an electrically conductive fluid, it interacts with the electromagnetic field which contributes to generating and sustaining the plasma state as expressed by Eqs. (2) and (4). Therein, the electric current density j in the plasma is given by the generalized Ohm's law as where s signifies electrical conductivity. The electromagnetic field is totally described by the Maxwell's equations as where 0 m is permeability of vacuum. Many simulations for thermal plasma flows have been carried out with a steady-state restriction by neglecting the transient terms in the governing equations. [8][9][10][11][12][13][14][15][16][17][18] However, fluid phenomena, even for air flow at RT, are naturally unsteady and often accompanied with fluctuations and vortex generation due to hydrodynamic instabilities. [19][20][21][22][23][24] To capture those features, time-dependent simulation is to be performed with a computational method capable of solving the governing equations accurately, as discussed and demonstrated in the following subsections. 2.1.2. Problem of turbulent-laminar coexistence. In general, a flow with a high Reynolds number tends to be turbulent, in which multiscale vortices are mixing the flow itself and transporting substances and thermal energy at a much higher rate than gradient-driven diffusions. On the other hand, a flow with a low Reynolds number hardly becomes turbulent, in which vortices are attenuated and eventually disappear. As a result, the flow tends to be laminar and convective mixing hardly occurs. 25) During plasma discharge, a part of cold nonionized gas touches and/or enters hot plasma, and then ionizes and becomes plasma. A part of hot plasma cools down and becomes back to cold gas again. Herein, the fluid experiences both a hightemperature state with high kinematic viscosity and a lowtemperature state with low kinematic viscosity (see Fig. B·1 in Appendix B). A low-temperature flow tends to have a high local Reynolds number and be turbulent having multiscale vortices larger than a few tens of microns. Meanwhile, a high-temperature plasma tends to have a low local Reynolds number and can have only vortices larger than sub-millimeters 26,27) which are usually comparable to or larger than the computational grid scale. In other words, when the cold gas flow with small vortices is ionized to plasma or enters the pre-existing plasma, the small vortices are attenuated; consequently, the turbulence in the fluid is suppressed. The turbulence model used in simulations should be able to express the turbulent-laminar coexistence including the turbulence transition and re-laminization processes.
To meet the unique condition, the computational demonstrations in this section adopted an Large Eddy Simulation (LES) approach, which implemented the coherent structure model 28,29) for the effects of vortices unresolvable by the computation grid. It is noteworthy that the use of the conventional Smagorinsky model, 30) which is often implemented in LESs, is not recommended because it gives a turbulent viscosity even in a local laminar region having a velocity gradient. The detailed discussion about turbulence models applicable to thermal plasma-nonionized gas coexisting flows is found in the literature. 27) 2.1.3. Problem of density variation. Thermal plasmanonionized gas coexisting flows have remarkable variations of mass density, e.g. approximately 45 times between 300 K and 12 000 K for argon gas, and over 1000 times between 300 K and 25 000 K for SF 6 gas including their plasma states. Nonetheless, the flows are almost incompressible or slightly compressible because of the low Mach numbers, e.g. 0.003-0.045 in an argon RFICP discharge system, and 0.04-0.5 in a DC transferred arc plasma in an SF 6 discharge system. Although those plasmas offer high flow speeds due to thermal expansion, their acoustic speeds are also high due to the high temperatures.
For those flows, the use of compressible flow simulation methods is not practical because they require a sufficiently small time-step width to resolve acoustic wave propagation and a large amount of computation time to reproduce the dynamics of turbulent vortices on an engineering-meaningful time scale even using high-performance computers. Therefore, an approach that is based on the incompressible flow simulation method and incorporates corrections to express compressibility is effective and practical regarding The Japan Society of Applied Physics by IOP Publishing Ltd a thermal plasma-nonionized gas coexisting system as an "almost incompressible flow with large density variation." However, under such conditions, numerical calculations with flow simulation methods used commonly are generally difficult to obtain the convergent solutions because of the numerical instabilities causing easy computational breakdown. For this reason, almost all previous thermal plasma flow simulations have adopted numerically stable methods that can obtain at least a solution with low accuracy. Those conventional methods discretize the convection terms in the governing equations by a low-order scheme which produces numerical viscosity damping numerical instability and stabilizing the calculation. Simultaneously, the effect also damps vortices and blurs the gradients of temperature and concentration. To the author's knowledge, there had been no reports of simulating vortexinduced entrainment of cold gas into a thermal plasma jet, complex fluid mixing associated with vortex breakup, and multiscale vortices spreading far from the plasma, after the 25 years since they were observed in the experiment 1) in 1991.
The author has been working on the development of computational methods that can perform stable and long-time flow simulations even with relatively large time-step widths while capturing turbulent vortices induced in thermal plasmanonionized gas coexisting flow systems. Using them, the computational methods demonstrated successful simulations of RFICPs in three-dimensional (3D) argon flow systems 26,31) in 2012 and 2013, a DC non-transferred arc plasma jet in an argon flow system in a two-dimensional (2D) frame 27) in 2016 coupled with nanoparticle formation-transport processes 32) in 2019, and in a 3D frame 33,34) in 2018 and 2020. Section 2.2 presents other sets of computational results related to these works.
2.2. Computational demonstrations 2.2.1. SF 6 flow system with DC transferred arc plasma. SF 6 (sulfur hexafluoride) gas has been used as a quenching gas for an arc plasma generated in high-voltage gas-blast circuit breakers when a short-circuit fault occurs in electric power systems. Compared to other gases (e.g. air), SF 6 gas offers higher quenching performance, i.e. extremely shorter time for arc plasma extinction. To seek highlydemanded alternative gas, the mechanism of fast plasma quenching in SF 6 gas should be elucidated.
For an arc plasma in SF 6 gas flow, the experiment 3) visualized a turbulent-like flow structure around the arc plasma sustained in a converging-diverging cylinder nozzle as a model-type gas circuit breaker. The experiment result suggested that turbulent vortices might play a significant role for arc plasma quenching by mixing cold nonionized gas into arc plasma effectively.
As mentioned in Sect. 2.1.3, the arc plasma has Mach numbers ranging from 0.04 to 0.5 which suggests that the flow is almost incompressible or weakly compressible in the arc plasma. On the other hand, the low-temperature SF 6 gas is an incompressible flow when it is injected, but it might locally become a compressible and/or supersonic flow around the arc plasma because it can receive the large momentum from the arc plasma and be accelerated strongly while it has a low acoustic speed.
To reproduce such an electrically conductive/nonconductive, incompressible/compressible, and subsonic/supersonic turbulent flow in an SF 6 arc-plasma-discharge system, the author has been developing an inhouse simulation code PLasma All-Speed Turbulence with Implicit Pressure Code (PLASTIPC) by combining the previous method Method- III 27) and a compressibility correction 35) with an additional unique modification. The details of PLASTIPC are described in Appendix C.
The numerical simulations were performed using two inhouse codes: the basic code and the PLASTIPC as described on Table I, under the same conditions as the experiment. 2) Figure 1(a) shows the geometry of the arc plasma discharge system composed of a copper cathode and a copper anode with their distance of 50 mm and a converging-diverging cylinder made of polytetrafluoroethylene (PTFE). The DC transferred arc plasma was generated around atmospheric pressure with the current of 54 A between the cathode and the anode. Therein, the following equation of current density conservation was solved with the aforementioned governing equations Here, the physical model assumed that the entire flow including the arc plasma was in an local thermodynamic equilibrium (LTE) state to exclude nonequilibrium effects which might contribute to the fast plasma quenching in SF 6 atmosphere. Furthermore, the computational domain was set intentionally in a 2D axisymmetric coordinate frame. This restriction excludes effects from 3D deformation of the arc plasma due to bending or twisting. It is also noteworthy that turbulent flows consist of vortices with three-dimensionally various orientations. Nonetheless, the 2D restriction would help to reveal the minimal factor of vortex formation at the arc plasma fringe which would trigger turbulence growth in the arc-plasma-discharge flow system. We believe that these sophisticated simulation conditions could highlight the fundamental features of the arc-plasma-induced turbulence of SF 6 flow and elucidate the underlying physics and mechanism of the vortex-related fast plasma quenching process.
The temperature-dependent data of thermodynamic and transport properties of SF 6 from 300 to 30 000 K considering dissociation and ionization 36) were provided by Dr. Yasunori Tanaka of Kanazawa University and implemented in the computation codes. The origin of the coordinates was set on the midpoint between the electrodes on the central axis. The computational domain was discretized uniformly by grid points with their distance of 0.1 mm in both directions. The time-step width was set to be 0.01 ms. Boundary conditions between arc plasma and electrodes were similar to the literature. 37) The flow near the angled wall was treated by the immersed boundary method with direct forcing. 38) Figures 1(b)-1(f) show the comparisons of the computational results obtained by the different codes for instantaneous distributions of the temperature, the vorticity, the flow vectors with colors indicating the Mach number, the acoustic speed, and the extracted edges of approximated refraction index (RI) gradient, respectively. The corresponding animations are also presented in Amination 1 as the supplementary data. It is noteworthy that the colors were assigned to the values on the logarithm scales for the acoustic speed and vorticity because of their large variations of the magnitudes. Here, fluid with its temperature over 8000 K was regarded as plasma.
The basic code and PLASTIPC exhibited the distinctively different results. PLASTIPC captured the small structures in and around the arc plasma in the cylinder for all the distributions. On the other hand, the basic code did not because of its damping effects of the 1st-order schemes as explained in Sect. 2.1.3. As shown in Fig. 1(b), PLASTIPC obtained the thinner arc plasma over 20 000 K, whereas the basic code obtained the wider arc plasma around 15 000 K. Under the LTE condition, the temperatures are convertible to the electron number densities which can be measured. Figure 2 presents the map of the arc diameters and the maximum electron number densities at the axial positions in −5 to 5 mm for the representative instantaneous data. The experiment data was provided by Dr. Yuki Inada of Saitama University and identical to those in the literature. 2) The plot points of the simulation by PLASTIPC distributed in the range of the experiment results for the arc diameters, whereas those of the simulation by the basic code did not. Although the maximum electron number densities of PLASTIPC showed values little higher than those of the experiment, the discrepancy was acceptable because it was caused by the 2D restriction to the arc plasma which resulted in the constriction of the arc plasma. Figure 3 also compares the time evolutions of the electron number density after the power supply was turned off. The multiple curves of the simulations were drawn from several different positions on the central axis. PLASTIPC obtained the similar decay curves to the experiment 2) which showed that the electrons of the arc plasma disappeared within 30 μs. Meanwhile, the basic code obtained much slower decay curves. These results indicate that PLASTIPC could simulate the mixing of high-and low-temperature fluids because it reproduced even small-scale vortices which contribute to turbulent-like mixing, as shown in the right part of Fig. 1(c) of a complex distribution with intermittent negative and positive vorticities in and near the arc plasma. Figures 2 and 3 rationally support the validity of the simulation by PLASTIPC. It is remarkable that those agreements between the experiment 2) and the simulation by PLASTIPC were obtained under the LTE assumption and 2D axisymmetric restriction. This fact suggests that the turbulent mixing effect following vortex generation in and near the arc plasma can dominate the fast quenching of the arc plasma itself rather than nonequilibrium or 3D effects. Figure 1(f) shows the extracted edges from the approximated RI gradients, which were not the same but similar to Schlieren patterns. They were depicted to emphasize the density nonuniformity attributed to the spatial variations of temperature and pressure. Those patterns highlighted that the complex flow structure was formed in the narrowest and diverging zones in the results by PLASTIPC. As the time passed, those complex patterns were conveyed downstream as shown in Animation 1 (f). The advection speed obtained by PLASTIPC were estimated to be approximately 20-30 m s −1 . The speed of the band-pass Schlieren patterns in the experiment 3) ranged from 20 to 26 m s −1 . They also agreed with each other.
The computational demonstration and discussion in this section support that PLASTIPC can be a useful computationally-analysis tool for turbulent-like thermal plasmanonionized gas coexisting flows in all-speed range from incompressible to compressible and from subsonic to supersonic conditions.

2.2.2.
Water DC non-transferred arc plasma jet for metal separation as nanoparticles. Thermal plasma of water-originated species, called "water plasma," is also an attractive medium for innovative green technology. By virtue of its high enthalpy and reactive components, it can be used for waste treatments such as decomposition of harmful substances and recovery of useful substances from even persistent wastes. Watanabe and his colleagues have reported their experimental works about successful decomposition of water-soluble organic compounds [39][40][41] and biomass utilization. 42) Water-solution waste includes not only organic compounds but also metals occasionally. When a water plasma jet is generated directly using such solution waste as the source, it is presumed that metal substances are transported with the high-  The Japan Society of Applied Physics by IOP Publishing Ltd temperature plasma flow and form aerosol at the plasma fringe. Aerosol particles formed in such conditions tend to be nanoparticles. 43) Therefore, this metal separation process can also be regarded as a nanoparticle production process.
Not only nanoparticles but also vortices start to be generated at the plasma fringe, where strong velocity shear and a steep temperature gradient occur. The high temperature for rapid vaporization of the raw material is one of the advantages of thermal plasma; in addition, the steep temperature gradient is also an important advantage for nanoparticle production. The temperature gradient corresponds to the cooling rate at which the material vapor is converted into nanoparticles. The various studies indicated that smaller nanoparticles were obtained when the cooling rate was higher. [15][16][17]44,45) Thermal plasma-induced vortices play one of the dominant roles in the nanoparticle formation process, because the vortices entrain cold ambient gas into the hot plasma and affect the local temperature and number densities of the material vapor atoms and nanoparticles.
The following is an overview of the nanoparticle formation at the plasma fringe: (i) the material vapor in the plasma is transported to the plasma fringe where the temperature decreases, (ii) the material vapor becomes supersaturated, (iii) critical nuclei are formed (homogeneous nucleation), (iv) the material vapor condenses on the nuclei (heterogeneous condensation) and nanoparticles grow, and (v) simultaneously the nanoparticles collide and unite forming larger nanoparticles (interparticle coagulation).
In engineering time-and spatial-scales, the aerosol dynamics approach effectively describes the growth and transport processes of nanoparticles. We have several types of the mathematical models. Here, one of the sophisticated models [32][33][34]46) is used. It assumes: (a) nanoparticles are liquid spheres, (b) the temperature of nanoparticle is identical to that of bulk fluid surrounding them, and (c) electric charge effects are neglected although particle charging are occasionally considerable. 47) The governing equations of the model are written as where n is number density, D is diffusion coefficient.
Subscripts v, c, p, and s denote vapor, critical state, particle, and saturated state, respectively. The variable f is defined as f n k, where k is the average monomer number in a nanoparticle. The diffusion coefficients of vapor D v and nanoparticles D p are derived from the molecular theory 48) and the literature, 49) respectively. The homogeneous nucleation rate J u and the number of monomers composing a nanoparticle in a critical state k c are estimated by the modified selfconsistent nucleation theory. 50) 0 b is the parameter related to collision frequency derived in the literature. 51) K th is the thermophoresis coefficient. 52) This section presents a computational demonstration of a DC non-transferred arc plasma jet generated from water vapor mixed with gold vapor diluted, as shown in Fig. 4(a). The plasma jet with a maximum temperature of 5754 K and a maximum velocity of 1119 m s −1 with a parabolic distribution are steadily ejected from a nozzle exit with a 1.8 mm inner diameter and 3.6 mm outer diameter into the computational domain with water vapor atmosphere at 400 K. Gold is supposed to be already vaporized and included as vapor atoms in the plasma jet at a rate of 0.1 g min −1 . Under this sophisticated condition, the chemical reactions such as oxidation can be neglected. This process will obtain gold nanoparticles which have high biocompatibility and are expected to be applied to photodynamic therapy, photothermal therapy, X-ray computed tomography imaging, drug delivery, and biosensors. 53) The origin of the coordinate system was set at the center of the nozzle exit. A staggered grid system with the width of 0.02 mm in each axis direction was used. The time-step width was set to be 0.02 ms. Because the Mach numbers in the flow field were on the order of 10 −3 to 10 −2 , the flow was incompressible in the whole computation domain. Therefore, the simulation was performed using PLASTIPC without compressibility correction. For the transient terms, the 3rdorder Adams-Bashforth-Moulton scheme 54) was used. The temperature-dependent data of the net radiation rate and transport properties implemented in the simulation were obtained from the literatures. 55,56) Figure 4 and Animation 2 show the computation results. The dynamic behavior of the temperature field showed the obvious entrainment of ambient cold gas into the plasma jet. The vortex structures were represented by Q* which was the isosurfaces of the second invariant of the velocity gradient tensor normalized by the mean flow velocity at the nozzle exit and its diameter. The water plasma jet was straighter, and the flow state was less turbulent than the argon thermal plasma jet which was obtained by the previous simulations, 33,34) because of the difference of the Reynolds numbers. The vortex structures indicated that vortex rings were generated regularly due to the Kelvin-Helmholtz instability in the upstream region, an upstream vortex caught up with a downstream vortex, and the vortices deformed and merged with each other in the downstream region.
Supposed that the temperature distribution was considered as streak patterns, the interface region appeared to form a single layer of the vortex rows. However, the patterns on the cross-sectional plane of x 2 = 0 showed clearly that this water plasma jet generated a two-layered structure with thick vortex rings from the inside of the high-temperature plasma jet to the interface region and thin vortex rings in the surrounding lowtemperature region.
In the high temperature region of the water plasma jet, gold nanoparticles did not exist but were in the vapor state. Because the gold vapor was converted to gold nanoparticles at the plasma fringe, many gold nanoparticles existed outside the plasma. Meanwhile, the number of the nanoparticles SL0801-6 © 2023 The Author(s). Published on behalf of The Japan Society of Applied Physics by IOP Publishing Ltd decreased in the downstream region. This was attributed to not only a dilution effect due to diffusion but also collisional growth through interparticle coagulation.
The distribution of gold vapor exhibited island-like high number density zones. These zones were formed because the mass density of fluid including gold vapor increased by the decrease in temperature as a result of the entrainment of the ambient low-temperature gas. Furthermore, the two-layered structure of the vortex rings in and out of the plasma jet indicated that this flow system had two convective transport regimes with different scales: one was the gold vapor transport by large vortices in and near the plasma, and the other was the gold nanoparticle transport by small vortices outside the plasma.

2.2.3.
Argon RFICP with an off-axis DC non-transferred arc plasma jet. An RFICP which is produced under atmospheric pressure offers a large volume thermal plasma. Moreover, a higher performance plasma is obtained by an injection of a DC non-transferred arc plasma jet into an RFICP. This so-called "DC-RF hybrid plasma" has been used as a powerful heat source to achieve high throughput of nanoparticle production. 57,58) However, the plasma is still too difficult to control; actually, it forms a very complex flow including multiple vortices as shown by a particle tracking velocimetry. 59) Timedependent 3D simulations also revealed the dynamic motions of multiscale vortices in the torch. 26,31) The results showed additionally that the thermofluid field had a non-axisymmetric profile with separate high-and low-temperature regions even though a thermal plasma jet was injected along the central axis of the torch.
This section presents another computation result to give new insights into a 3D complex vortex structure in a torch where an off-axis thermal plasma jet was injected intentionally to a cold gas region beside an RFICP, because a vortex structure was an indispensable information for effective flow control from a viewpoint of fluid mechanics. The Japan Society of Applied Physics by IOP Publishing Ltd Figure 5(a) shows a schematic illustration of the DC-RF hybrid plasma torch. The RFICP was generated with a frequency of 4 MHz by a spiral coil under atmospheric pressure. From the torch inlet, argon sheath gas with 300 K was injected at 23 l min −1 along the cylindrical quartz wall whose radius was 25 mm. An argon DC non-transferred arc plasma jet with 12 000 K was injected at 3 l min −1 along a 16 mm off-axis location of the torch inlet. The jet injector had a radius of 4 mm.
The flow field, where the plasma at a high temperature and a cold gas at RT coexisted, was treated simultaneously considering large variations of the transport properties and the density with several orders of magnitude. Meanwhile, the Mach number was estimated to range from 0.003 to 0.045. Therefore, the partially ionized flow in the torch was simulated as an incompressible flow to obtain the solution in a practical time scale.
The computational domain was represented with approximately 2 million control volumes by a Finite Volume Method. The so-called "blending scheme" was used to discretize the convection terms. 60) This scheme blended the central differencing scheme with the second-order accuracy and the upwind differencing scheme with the first-order accuracy as g Herein, the value of 0.9 b g = was set. It is noteworthy that this "blending scheme" is different from "hybrid differencing scheme" proposed in the literature. 61) Figure 5(b) presents the instantaneous temperature profile in the torch. Animation 3(b) shows the transient behavior of the plasma mainly after the DC non-transferred arc plasma jet at the time of 0.05 ms. The high-temperature plasma jet mainly hit and merged with the pre-existing RFICP. Meanwhile, a part of the plasma near the torch wall absorbed the electromagnetic energy, increased its temperature due to the Joule heating, and touched the torch wall instantaneously. Figure 5(c) and Animation 3(c) show the vortices generated in the torch. Even before the jet injection, larger vortices with high temperatures were formed by the Lorentz force in the plasma. On the other hand, cold vortices exhibited small or thin forms. Those features agree with the prediction from the Kolmogorov theory as explained in the previous reports. 26,27) After the plasma jet was injected, the Lorentz force bent the plasma jet inward; consequently, large hightemperature vortices which stayed in the coil zone were stretched downstream. Around the time was 0.20-0.25 ms, a row of cold thin vortices was generated surrounding the plasma jet because of Kelvin-Helmholtz instability. This "rib-like structure" is a key feature observed at turbulence transition in nature. After the time was 0.30 ms, many cold vortices were generated and formed a turbulent-like field.

Binary alloy nanoparticles
Nanoparticles have attracted interests of considerable attention in the scientific community and industry because they offer unique characteristics that differ greatly from those of bulk materials or powders composed of larger particles. 62) In particular, nanoparticles of binary alloys are anticipated as potentially useful materials. For instance, nanoparticles of molybdenum-disilicide MoSi 2 and titanium-disilicide TiSi 2 have been synthesized 63) for extremely small electronic and mechanical applications such as solar-control windows, electromagnetic shielding, and contact materials in microelectronics. Titanium-diboride TiB 2 displays high strength, durability, mp, hardness, and wear resistance. Therefore, it has been applied as crucibles, cutting tools, impact-resistant armors, and wear-resistant coatings. [64][65][66] Since nanoparticles of TiB 2 have electrical conductivity 67) and low work function, 68) they are also expected to be used for the electromagnetic shielding and solar control windows interacting with IR and UV lights. Furthermore, magnetic nanoparticles also play an important role in numerous applications such as magnetic data storage, ferrofluids, catalysis, and biotechnology. [69][70][71][72][73][74] Because raw materials have high mp or bp, high-rate synthesis of binary alloy nanoparticles is almost impossible using conventional methods such as liquid-phase deposition or grinding techniques. Even combustion processes are ineffectual because the combustion flame cannot reach sufficiently high temperatures to vaporize the raw materials; moreover, the oxidation atmosphere for combustion causes unfavorable production of contaminants.

Model description
Binary alloy nanoparticles are formed from two kinds of metal vapors through binary-component homogeneous nucleation and heterogeneous co-condensation under a cooling environment of the plasma tail. The nucleation and condensation take place simultaneously, but their characteristic time scales are different by several orders of magnitude. Furthermore, the nanoparticles also grow collectively, coagulating among themselves, with different diameters ranging by two to three orders of magnitude. For such a spatiotemporal multiscale process, numerical calculation based on molecular dynamics is practically impossible to treat the growths of all nanoparticles due to the limitations of current computer performance, because they can only track the formation process of a few dozen nuclei for a few nanoseconds. 90) Therefore, we have constructed a binary aerosol formationgrowth model mathematically and developed a unique computational method to calculate the model. 84 The Japan Society of Applied Physics by IOP Publishing Ltd model was based on the following assumptions: (i) nanoparticles are spherical; (ii) temperatures of material vapors and nanoparticles are identical; (iii) heat generated by condensation and electric charge of nanoparticles are negligible; (iv) material vapors are regarded as ideal gases; and (v) the mp of a nanoparticle depends on its size and composition which is supposed to be uniform in the particle. The profile of numerous nanoparticles was expressed by a particle sizecomposition distribution (PSCD). During computation of the time evolution of the PSCD, the distribution was treated with the manner of two-directional nodal discretization. 91) Therein, the net production rate of nanoparticles with the volume v k and the second-metal content The first and second terms in the right-hand side mean the contributions of homogeneous nucleation and heterogeneous condensation, respectively. The third and fourth terms represent the gain and loss by interparticle coagulation. x and y are splitting operators in the size and composition directions, considering conservation of particle volume. 88 where d is diameter, D is diffusion coefficient, and Kn is Knudsen number defined as the ratio of the gas mean free path to the particle radius. Subscript vap signifies vapor. This equation is applicable to condensation on multiscale particles from the sub-nanometer to a few micrometers. a denotes accommodation coefficient that can depend on the materials' mixing ratio and curvature at the particle's surface. Because the exact values of a are still unknown, the values were set to be 0.9 on a liquid particle and 1.0 on a solid particle. N s M ( ) signifies saturated vapor concentration considering the effects of material mixture and surface curvature. 93) The population balance equations of material vapors are also calculated simultaneously because the number density of material vapor atoms directly affects the rates of nucleation and condensation: Additionally, a model of mp depression due to the nanoscale size and mixture effect was combined with the formation-growth model. For inorganic materials, the mp T MP of a particle decreases linearly with the particle diameter d, and this effect is considerable especially at nanometer scales: 94,95) where T MP,bulk represents a mp in a bulk state, and parameter e is a characteristic property determined by the solid and liquid surface energies and bulk melting enthalpy.

Computational demonstration
This section presents the recent work on the computational analysis for a collective formation process of Co-Sm (cobalt-   94,95) X denotes molar fraction. This model also considered the dependency of the material composition on T .

MP,bulk
The detailed discussion for that treatment is found in the leterature. 78) The material properties of Co and Sm were basically obtained from the data book. 96) The temperaturedependent density and saturation pressure of Sm were obtained respectively from the literatures. 97,98) The computation was performed under a typical cooling condition at the plasma tail, where the temperature monotonically decreases at a rate of 5.0 × 10 4 K s −1 . Corresponding to that cooling rate, the time-step width t D for integration was set to 2.0 × 10 −6 s, which gave sufficient resolution for the process. The total feed rate of the precursory mixed powder with Co: Sm = 5: 1 was set at 0.3 g min −1 in the carrier argon gas at 3.0 l min −1 . It was assumed that after the raw material powder was vaporized completely by plasma, the vapor molecules were transported and diluted by half at the position of 2100 K, which was the starting point of the computation. At 1585 K, most of the particles with diameters of around 50 nm have Sm contents of 1 to 2 at%. In this stage, the particles smaller than 150 nm were in a liquid state and those larger than 150 nm were in a solid state. Because vapor molecules of Co and Sm condensed simultaneously on both solid large particles and liquid small particles as indicated in Figs. 6(e) and 6(f), most of the particles were created as liquid Co-Sm droplets with condensed Co and Sm as shown in the inset schematic illustration in Fig. 6(a). At 1500 K in Fig. 6(b), the particles with diameters of around 40 nm and Sm contents ranging from 3 to 4 at% were in a solid state; at the same time, liquid particles with size of approximately 35 nm existed as well. Condensation of Co was practically completed at 1500 K, and only Sm continued to condense on the liquid and solid particles as indicated in Figs. 6(e) and 6(f). Consequently, liquid or solid Co-Sm particles with condensed Sm were produced. At 1400 K in Fig. 6(c), almost all particles were solidified. Most of the particles had diameters from 50 to 80 nm and Sm contents from 7 to 8 at%. In Fig. 6(f), the conversion ratio of Sm was 30%, indicating that Sm molecules were still in the vapor phase and would condense on the solid particles which had already formed. Hence, the solid Co-Sm particles were covered with thick condensed Sm at this temperature. At 1100 K in Fig. 6(d), condensation of Sm was completed as indicated in Fig. 6(f), and all the particles were solid. Most of the particles had diameters from 50 to 80 nm and Sm contents from 16 to 17 at%. At temperatures below 1100 K at which condensations of both Co and Sm were completed, it was presumable that 90% of the particles were Co-Sm-core/Sm-shell particles and 10% were Co-Sm alloy particles as the final products. Figures 7(a) and 7(b) show the particle size distributions obtained from the computation and the experiment, respectively. It is noteworthy that both counted the particles larger than 13 nm. They agreed with each other for their average diameters, standard deviations, and shapes of the profiles, which supported the validity of this computational work. The original literature 78) discussed the discrepancy between them as well, suggesting that this model still had rooms to be improved especially for treating core-shell-structured nanoparticles. Nonetheless, it is noteworthy that this model could obtain reliable calculation results for conditions under which two vapors completed their condensations before nanoparticles were solidified. 77,84,88) 4. Computational fluid mechanics for arc welding 4.1. Transport phenomena in arc welding Arc welding is an indispensable manufacturing technology used in a wide variety of industrial fields, including automobiles and other transportation equipment, civil engineering, construction, and plant facilities. There are several types of arc welding: Gas Tungsten Arc (GTA) welding, in which tungsten with a high mp is used as an electrode, Gas Metal Arc (GMA) welding, in which the metal wire is used as an electrode and it is intentionally melted to weld and join the base metal, and Submerged Arc (SA) welding, in which an electrode is placed in a flux powder.
In any of these methods, the key to welding quality is the heat transport process by molten metal flows of the electrode and the weld pool interacting with arc plasma at temperatures exceeding 10 000 K. The molten metal flows involve complex phenomena such as phase change, deformation due to interaction with arc plasma, electromagnetic force (Lorentz force), resistance heating (Joule heating) generated by electric current, and surface shear force (Marangoni effect) caused by the gradient in surface tension. 7) Precise understanding of the heat and mass transfer processes in the molten metal flows as well as arc plasma are required to achieve higher-performance welding with accurate control.
The characteristics of arc plasmas during welding have been elucidated by experimental and computational studies. 99) In particular, the imaging spectroscopy techniques revealed the temperature and concentration of species in arc plasmas even contaminated with metal vapors. [100][101][102][103][104][105][106][107][108][109] Moreover, oxygen contamination was measured on a water-cooled cupper 110) and investigated by optical visualization techniques of Schlieren method and xenon emission filtering. 4,5) Simulations were also used to obtain numerical information in arc plasmas, which could hardly be measured. 8,9,[111][112][113][114][115][116][117][118] The molten metals produced during arc welding have also been investigated. The melting and deforming behaviors of weld pools and molten fluxes in SA welding were captured by an Xray transmission imaging technnique 119,120) to support the mechanism clarifications. 121,122) The deforming behaviors of molten electrodes were observed directly using high-speed video recording during a SA welding 120) and a GTA welding. 123,124) Especially in a GTA welding, the transport of emitter additives in a tungsten electrode is also important because it affects the stability of arc discharge and the lifetime of the electrode. The investigations were carried out experimentally 125,126) and computationally. 127) Additionally, the thermal energy transferred from an arc plasma to a base metal was measured for a GTA welding condition. 128) Despite those efforts, it is still difficult to measure the temperature and flow velocity in the molten metals.
Computational fluid mechanics for molten metals as well as arc plasmas have provided new insights into the invisible flow dynamics and heat transports. Section 4.2 describes Smoothed Particle Hydrodynamics (SPH) method which is expected to simulate complex molten metal flows with deformation and phase change in a way simpler than other computation methods. Section 4.3 presents a computational demonstration applying the SPH method to the melting and breaking-off process of an Flux-Cored (FC) wire.

SPH representation
Flow simulations often use computational grid systems; for arc welding simulation, complicated calculation procedures are required to capture movement of solid-liquid interfaces due to phase change and free surface deformation of the molten metal. To overcome such arduous works, we have been employing the SPH method which is a gridless Lagrangian-based method originally proposed to simulate fluid dynamics in astrophysics. 129,130) Therein, the flow is represented by motions of particles with mass and additional physical quantities. The quantity f at a certain position r a  is where T , C, and Q a  denote temperature, specific heat, and net heat generation, respectively. To simulate incompressible flows, the density homogenizing algorithm 132) is additionally implemented into the conventional SPH method. We have applied this incompressible SPH method to the thermofluid simulations of molten metal flows in several kinds of arc welding processes such as GTA weding, 133,134) GMA welding, [135][136][137][138][139] SA welding, 140,141) and Flux-Cored Arc (FCA) welding. 142) Because this incompressible SPH method offers high adaptability for other types of metal processing, it has also revealed the molten metal formation processes with convective and diffusive heat transports during dissimilar Resistance Spot (RS) welding 143,144) and Electroslag (ES) welding 145) collaborating with the experiments, 146,147) and dross formation during gas cutting. 148) 4.3. Simulation of FC wire droplet detachment FCA welding is a welding process which uses an FC wire as an electrode. An FC wire consists of a metal pipe and flux filled in the pipe. The FCA welding can be applied to all welding positions. During FCA welding, a flux column is formed at the tip of the wire under specific welding conditions, and then molten metal droplets are transferred to a weld pool along the flux column. It was reported that the flux column formation contributed to transfer of uniform droplets and reduction of defects. 149) However, the mechanism of such flux column formation is still poorly understood. This section presents a computational demonstration of a flux column formation and droplet transfer at a TiO 2 -based wire in FCA welding, using the incompressible SPH method. 142) The diameters of the whole wire and the flux region were set to be 1.2 mm and 0.54 mm, respectively. The wire feed rate was set at 16 m min −1 . It was supposed that the welding was conducted with a DC current of 300 A at electrode positive polarity. The shielding gas was pure argon supplied at 20 l min −1 . The diameter of each calculation particle was set to be 0.05 mm. The time-step width was chosen to be 0.05 ms. The details of the computational conditions are found in the literature. 142) Figures 8(a) and 8(b) depict the metal vapor fraction and temperature distributions in and around the arc plasma, which was obtained by a preliminary grid-based steady simulation in a 2D axisymmetric frame. Typical distributions for GMA welding were obtained. The metal vapor was generated mainly from the anode wire and conveyed downward with the flow. Because of the large radiation loss from the metal vapor, the temperature near the central axis was lower than the surrounding region. As a result, the highest temperature region around 15 000 K were produced not below but near the side of the electrode wire. Therefore, the electrode wire was heated most at its side surface by the arc plasma. Although the arc plasma actually fluctuates with the cathode wire deforming, those steady profiles of the arc plasma were imposed as the boundary conditions to the SPH simulation of an FC wire for simplicity. Figures 8(c) and 8(d) show the state and temperature distributions of the calculation particles representing the FC wire. The corresponding movies are found as Animation 5 in the supplementary data. The SPH method simulated a complex motion of the two-phase four-component medium which involved the phase changes, liquid flow, heat transport, large deformation, metal-flux mixing, two-component droplet formation, and eventual detachment. Because of the lower thermal conductivity and higher specific heat of TiO 2 , the thermal energy supplied from the arc plasma through the metal pipe surface was hardly conducted to the inside of the flux column; besides, the temperature hardly increased. Consequently, the flux remained a solid state and formed a column.
Figures 9(a)-9(c) display the flux column lengths, droplet diameters, and droplet detachment frequencies, respectively. The average length of the flux column was estimated by averaging 10 measurements. The distance between the left edge and the right edge of the droplet immediately after detachment was defined as the droplet diameter. The average The Japan Society of Applied Physics by IOP Publishing Ltd droplet size was also obtained by averaging 10 measurements. In the figures, the solid circles represent the average values. The error bars indicate the 95% confidence interval for the average values. For all the results, the computation and the experiment showed the close average values and the overlaps of their 95% confidence intervals, which supported the acceptable validity of this computational work. The discrepancies were caused by the model simplification and irregular dynamic behaviors of the melting FC wire and the arc plasma in the experiment. The more discussion with the detailed description of the experiment is found in the original literature. 142)

Concluding remarks and perspectives
This article has reviewed and discussed the recent progresses of the studies with the concept of "Computational plasma fluid mechanics." The inhouse simulation code PLASTIPC and prototypes have successfully captured hydrodynamic instabilities and reproduced plausible flow dynamics in the thermal plasma -nonionized gas coexisting systems. The unique computational method has made it feasible to study the collective formation and growth processes of binary alloy nanoparticles by means of numerical analysis. The SPH method with an incompressibility modification has achieved complex thermofluid behaviors of molten metal involving phase change, flow, heat transport, material mixing, and large deformation during arc welding.
As shown in the figures and animations, visual information aids in a person's first intuitive understanding. Moreover, computational simulation possesses numerical data of all quantities abiding by underlying physics, which cannot be obtained by experimental ways.
Plasma fluid systems are particularly difficult to measure the physical quantities. For instance, we do not have sufficient knowledge of the plasma fringe, where vortices mix the plasma and nonionized gas and/or nanoparticles are generated. The thermofluid dynamics at the plasma fringe and a consequently produced turbulent flow crucially affect the growth process of nanoparticles and cause ambient gas contamination in arc welding. Therefore, it is essential to study thermal plasma processes as a comprehensive fluid system in which hot plasma, cold nonionized gas, and materials coexist and interact with one another. For that, the viewpoint and approaches of fluid mechanics in addition to plasma physics will be indispensable. Especially, computational study will play a more important role in giving us new and deeper insights.

Acknowledgments
The works presented in this article were supported partly by Japan Power Academy (TOKUBETSU SUISHIN KENKYU 2019) and the Japan Society for the Promotion of Science Grant-in-Aid  Figure B·1 shows the temperature-depending kinematic viscosities of air, SF 6 , and Ar. The curves of air and SF 6 were obtained by converting the data, 36) and that of Ar was obtained by converting the data. 150) (a) (b) (c) where U u 1 º represents a velocity component in x 1 direction, f signifies a variable, and x 1 D represents a distance between the adjacent discrete points. h g stands for a weight parameter which is determined automatically in computation. 151) Subscript m denotes the mth discrete point.
PLASTIPC treats the transient terms selectively by Adams methods such as Adams-Bashforth, Adams-Moulton, and Adams-Bathforth-Moulton methods. 54) As presented in Sect. 2.2.1, the arc-plasma-discharge flow system in SF 6 gas was simulated well by the 3rd-order Adams-Moulton method described as where n ( ) denotes a time-step number, and asterisk * means an intermediate value to convergence.
In the implicit pressure-velocity coupling algorithms such as Semi-Implicit Method for Pressure-linked Equations (SIMPLE) and Pressure-Implicit with Splitting of Operators (PISO) including IPISO (Improved PISO), 152) solving the pressure correction equation works effectually for incompressible flows. To simulate compressible flows based on these algorithms, the density variation r¢ due to the pressure variation p¢ at compression/expansion is included into the pressure correction equation. 35 35) modified by 1st-order upwind scheme 35) TVD 153,154) with limited TCDF 155) SL0801-14 © 2023 The Author(s). Published on behalf of The Japan Society of Applied Physics by IOP Publishing Ltd (total variation diminishing) scheme 153,154) as follows It is noteworthy that r 0 1 ( )   y was required to avoid computational breakdown because this TVD scheme was used not for the velocity equation but for only the compressible part in the pressure correction equation. The similar manner is also taken at the second pressure correction step for p in the PISO algorithm.
Because PLASTIPC works as a totally implicit solver, additional convection-diffusion equations can easily be implemented and calculated stably even with high rates of chemical reactions and/or nanoparticle formations such as Eqs. (12), (13), and (14).