Safeguarding against Inactivation Temperatures during Plasma Treatment of Skin : Multiphysics Model and Phase Field Method

One of the most appealing applications of cold plasmas is medical treatment of the skin. An important concern is the capability to safeguard the non-targeted cells against inactivation temperatures during the plasma treatment. Unfortunately, it is problematic to experimentally determine the highest transient temperatures in these cells during the plasma treatment. In the present work, a complete multiphysics model was built based on finite element analysis using phase field method coupled with heat transfer and fluid dynamics to study the discharge phenomenon of cold plasma with helium carrier gas ejected out of a tube for skin treatment. In such plasmas with carrier gas, the fractions of plasma constituents are small compared to the carrier gas, so thermofluid analysis is needed for the carrier gas as the major contributor to the fluid and heat flow. The phase field method has been used to capture the moving helium gas in air, which has enabled us to compute fluid dynamics parameters for each phase individually. In addition to computational fluid dynamic analyses, we have also considered heat transfer in the fluids and to the skin using the Fourier law of heat conduction, which led to a multiphysics system. In the present paper, various flow velocities and tube-to-target distances (TTDs) have been considered to reveal the dependence of the fluid discharge output parameters on the flow and efficiency of heat transfer to the skin and the surrounding environment. The built model is a useful tool for future development of plasma treatment devices and to safeguard the non-targeted cells against inactivation temperatures.


Introduction
The progress in plasma medicine research has been fast and fruitful in recent years [1][2][3][4].One of its most appealing medical applications is treatment of living cells and tissues with cold plasmas (e.g., [5][6][7][8][9][10][11][12][13][14][15]).However, potential thermal damage of the non-targeted cells is an important concern for medical treatment using cold plasma jets.In particular, temperatures beyond 40 • C can cause detrimental effects in cells [16].In recent development of cold plasma treatment devices, higher temperature plumes are not well studied and most of the available literature has focused on a temperature lower than 30 • C [9,17].Plasma plumes with temperatures higher than the inactivation temperature of the cells and biological molecules need to be studied, as complete reliance on low temperature plasma may lead to loss in the maximal treatment capability.The species in the generated plasma would gain higher kinetic energies with increased plume temperature, which would lead to increased flux of the plasma constituents over the treatment zone.Unfortunately, it is problematic to experimentally determine the highest transient temperatures in the non-targeted cells during the plasma treatment because of the influence of the plasma itself.It is also unreliable to use the surface temperature of the treated tissue after the plasma treatment has stopped to deduce the highest transient temperature during the treatment because of the rapid change in the temperature due to the transfer of heat into the surrounding air or due to blood perfusion.As such, in order to facilitate future development of plasma treatment devices, and to safeguard the non-targeted cells against inactivation temperatures during plasma treatment, it is pertinent to develop accurate computational methods, which can reveal such information.Schröder et al. [18] was one of the pioneer groups who set up a parameterized model to simulate the species densities generated by a plasma needle and also the heat transfer to a skin layer, with a view to improve the design of medical plasma tools while, at the same time, ensuring safe temperatures in the skin of the patients.However, the model developed by Schröder et al. did not consider the plasma plume and its transfer to the treatment zone and fluid dynamic analysis is necessary when considering such systems and its effect on the skin layer surface.
In the present paper, we will formulate a different model to simulate the discharge of the carrier gas from the tube of a cold-plasma jet system into the air under atmospheric pressure before the gas comes into contact with the target.The modeled system consists of two fluids which are continuously interacting with and moving into each other; this is a reasonable assumption since there is a weak ionization rate in the cold plasma discharges [19] and therefore it is possible to model the cold plasma discharge using its carrier gas.Shimizu et al. [20] pointed out that cold/non-thermal plasma discharges had very low gas temperatures due to a weak ionization rate for electron collisions and reactions.Moreover, the fractions of plasma constituents were found to be small and therefore no significant additional heat was produced as a result of reactions in the discharge, so the bulk temperature of the carrier gas would be the dominating factor in determining the plume temperature [17].We consider that the carrier gas of the cold plasma is helium.Helium is a commonly employed carrier gas thanks to its low breakdown voltage and the uniformity of its generated plasmas [21].The phase field method will be used to capture the moving gas in air, which will enable us to compute fluid dynamics parameters for each phase individually.The phase field method was previously applied in other research fields, such as in the simulation of two-phase fluid systems, where predictions agreed excellently with experimental results [22].In addition to computational fluid dynamic analyses, we will also consider heat transfer in the fluids and to the skin using the Fourier law of heat conduction → q = −k∇T, where → q is the local heat flux density, k is the thermal conductivity of the skin and ∇T is the temperature gradient, which led to a multiphysics system.It is remarked that blood circulation and the volumetric heat flow in metabolism are neglected here, which might lead to differences from a realistic skin.In our current work, a variety of flow velocities and tube-to-target distances (TTDs) have been taken into account to analyze the dependence of the fluid discharge output parameters on the flow and efficiency of heat transfer to the skin and the surrounding environment.The consideration of air along with plasma carrier gas provides us with information on the forced gas flow cooling due to presence of air at the initial stage.Moreover, the use of phase field method is inspired from the previous work of Witterstein et al. [23] in which the cell differentiation process is modeled as phase change in the biological materials using the phase field method; this idea was indeed interesting and presented a great methodology.Furthermore, the system has been solved using the finite element method (FEM), which was previously used in simulation of plasma needle [24].
The significance of the present work is to provide a multiphase flow model which can track the evolution of the carrier gas in air and to provide the temperature change in the skin layer and the surrounding air.Although the present study focuses on the plasma treatment of the skin, the methodology is generic and it is relatively easy to adapt it for plasma treatment of other tissues.The model can be used to help design better equipment with usage of their maximum capacity for plasma treatment while the non-targeted cells are safeguarded against inactivation temperatures during the plasma treatment.In addition, apart from the difficulties in transient temperature measurements during plasma discharge, it is important to thoroughly investigate the potential possibilities of increased plume temperature using complete multiphysics computer simulation models.Cold plasma devices already licensed by the respective authorities for medical treatment of the infected wounds and infective skin diseases included (1) kINPen MED (neoplas tools GmbH, Greifswald, Germany) [25,26]; (2) MicroPlaSter plasma torch system developed and built by the Max Planck Institute for Extraterrestrial Physics in Garching/Germany and the company ADTEC Plasma Technology Co. Ltd., Hiroshima, Japan/London, UK [27] and (3) PlasmaDerm VU-2010 device (CINOGY GmbH, Duderstadt, Germany) [28].Our developed model could help to assess potential risks of thermal damages to the skin that might be tedious to study experimentally.Moreover, the feasibility of using phase field method in dynamic tracking of cold plasma discharge in air is investigated by exhaustively benchmarking the present model with previously performed experiments; the obtained theoretical results are in a great agreement with the experimental measurements in which it proves the accuracy, reliability and completeness of the present theoretical model.Moreover, it needs to be noted that the multiphysics simulation of cold plasma discharge for skin treatment coupled with heat transfer requires a continuous interface tracking method such as the phase field method as presented in this paper.

Materials and Methods
The phase field method has been used to capture the evolution of the He gas in air.Moreover, different carrier gas temperatures are considered in the presents study; the range of carrier gas temperatures that are considered in the present work embrace the currently used limits and even temperatures that are higher compared to the widely used limits.Furthermore, considering the initial state of the system, the helium gas is placed in a tube located above the air domain while the skin layer is placed below the air domain, so that the discharged gas can come into contact with the skin after passing through the air in the so-called "mixing region" that contains air.In addition, the plasma flow is considered to be continuous; this can be closely related to the experiments in which the plasma reactor operates throughout the experiment.The system setup is shown in Figure 1, and the modeled zone that consists of discharge tube, air and skin is marked on the schematic diagram.
Math.Comput.Appl.2017, 22, 24 3 of 14 wounds and infective skin diseases included (1) kINPen MED (neoplas tools GmbH, Greifswald, Germany) [25,26]; (2) MicroPlaSter plasma torch system developed and built by the Max Planck Institute for Extraterrestrial Physics in Garching/Germany and the company ADTEC Plasma Technology Co. Ltd., Hiroshima, Japan/London, UK [27] and (3) PlasmaDerm VU-2010 device (CINOGY GmbH, Duderstadt, Germany) [28].Our developed model could help to assess potential risks of thermal damages to the skin that might be tedious to study experimentally.Moreover, the feasibility of using phase field method in dynamic tracking of cold plasma discharge in air is investigated by exhaustively benchmarking the present model with previously performed experiments; the obtained theoretical results are in a great agreement with the experimental measurements in which it proves the accuracy, reliability and completeness of the present theoretical model.Moreover, it needs to be noted that the multiphysics simulation of cold plasma discharge for skin treatment coupled with heat transfer requires a continuous interface tracking method such as the phase field method as presented in this paper.

Materials and Methods
The phase field method has been used to capture the evolution of the He gas in air.Moreover, different carrier gas temperatures are considered in the presents study; the range of carrier gas temperatures that are considered in the present work embrace the currently used limits and even temperatures that are higher compared to the widely used limits.Furthermore, considering the initial state of the system, the helium gas is placed in a tube located above the air domain while the skin layer is placed below the air domain, so that the discharged gas can come into contact with the skin after passing through the air in the so-called "mixing region" that contains air.In addition, the plasma flow is considered to be continuous; this can be closely related to the experiments in which the plasma reactor operates throughout the experiment.The system setup is shown in Figure 1, and the modeled zone that consists of discharge tube, air and skin is marked on the schematic diagram.

Phase Field Theory Coupled with Heat Transfer
The phase field method is one of the widely used methods in interface tracking.The present theoretical model involves two-phase fluid systems in which the carrier gas (helium) discharges out of the tube with a constant velocity and it leads to formation of a velocity field in the mixing region.The phase field method is used to accomplish the interface tracking between these fluids that are located above the skin layer.The theoretical approach regarding this interface tracking method is shown in this section.The evolution of the phase field variable (Φ) is explained using the Cahn-Hilliard equation:

Phase Field Theory Coupled with Heat Transfer
The phase field method is one of the widely used methods in interface tracking.The present theoretical model involves two-phase fluid systems in which the carrier gas (helium) discharges out of the tube with a constant velocity and it leads to formation of a velocity field in the mixing region.The phase field method is used to accomplish the interface tracking between these fluids that are located above the skin layer.The theoretical approach regarding this interface tracking method is shown in this section.The evolution of the phase field variable (Φ) is explained using the Cahn-Hilliard equation: where γ represents the mobility, λ is the density of mixing energy and ε represents the interface thickness width.Moreover, when the term ∂φ ∂t + u • ∇φ is equal to zero, it will be converted to the well-known evolution function.The value of the phase field variable Φ would be either 1 or −1 in regions in which there is a fluid-fluid interface.Considering the present case where there is an interface between the plasma carrier gas and air, the phase field variable would take different values on different sides of this interface.Furthermore, the free energy density of the system can be expressed as a function of the dimensionless phase field parameter (Φ) as shown in Equation ( 2): where f tot is the total free energy density of the system that has the international system (SI) unit of (J/m 3 ).Moreover, the involvement of the free energy density of the system in the evolution of the phase field parameter is shown in Equation ( 3): On the right-hand side of Equation ( 3), the relaxation time that is controlled by the mobility (γ) aims to minimize the total free energy.The SI unit of the mobility in this case would be (m 3 •s/kg).
In the present two-phase fluid system, the free energy density of this mixture (plasma carrier gas and air) can be expressed as the sum of the elastic and the mixing energy.Furthermore, the consideration of the phase field parameter would be crucial due to the two-phase nature of the system.Therefore, the free energy density of the mixture as a function of Φ and ∇Φ are expressed as shown in Equation (4).
It is remarked here that Equation (4) assumes the Ginzburg-Landau form: The parameter λ represents the mixing energy density that has the SI unit of (N), and the parameter ε represents the interface thickness with SI unit of (m).The volume fractions of the fluid components in the system (plasma carrier gas and air) can be defined as (1 + Φ)/2 and (1 − Φ)/2, which are denoted as α and β, respectively.The surface tension considered in the present model can be defined using the mixing energy density (λ) and the thickness of the interface or in order words capillary width (ε): where σ represents the surface tension that has the SI unit of (N/m).The value of ε (interface thickness) was set to be half of the maximum mesh size, which, in our case, was 0.265 mm and the value of λ was set to be 0.2811 N, which led to better numerical convergence.The values of λ and ε depended on the grid/mesh resolution and on the computational resources available.Different values could be chosen for these parameters provided that numerical convergence was achieved.
Considering Equations ( 1) to (5), it is possible to re-write the Cahn-Hilliard equation that governs the phase field variable as: where G represents the chemical potential term with SI unit of (Pa): The Cahn-Hilliard equation shown in Equation ( 6) is a fourth-order partial differential equation (PDE) and is decomposed into two second-order PDEs for numerical analysis.Further derivation, decomposition and finite element formulation that are used in the present work are adopted from the previous work of Yue et al. [29].In the present study, a seamless triangular mesh was employed with the maximum element size of 0.53 mm and the minimum element size of 0.003 mm.The curvature factor of the mesh was set as 0.3 and linear discretization was employed in the present model.The mesh figures for the two configurations with TTDs of 1 and 5 mm are shown in Figure 2.  In addition, different mesh sizes were studied and no changes in the convergence efficiency and in the final obtained results were observed.Considering the fluid dynamics of the system, the Navier-Stokes momentum equation is: where is density, P is pressure, u is the velocity and F is given by = + + , ρg is the body force, Fst is the surface tension and Fvf represents the volume force.Moreover, the continuity equation that is solved alongside with the Navier-Stokes momentum equation is: In order to solve these equations numerically, a no-slip boundary condition is considered at the inner-boundary of the discharge tube and also at the upper surface of the skin that is exposed to the plasma plume; the no-slip boundary condition arises mainly as a result of the adhesion between the fluid and the solid surface.The regime was considered laminar since the largest Reynolds number was ~14, which was far below the number of ~4000 in a turbulent regime.Furthermore, the heat transfer in the mixing region (the domain above the skin surface) is explained using the heat equation, with the use of the temperature-dependent thermophysical properties of air and plasma carrier gas: where Cp is the specific heat capacity, represents the density and k is the thermal conductivity.Moreover, the Navier-Stokes momentum equation (Equation ( 8)) and the continuity equation (Equation ( 9)) coupled with the heat transfer (Equation ( 10)) have been used in our previous work [30], but the level-set interface tracking method was used there to explore the plasma evolution in air and water.Furthermore, in order to understand the coupling method used to determine the skin temperature during plasma exposure, one needs to examine (1) the temperature in the mixing region (Tmix) and ( 2) the temperature in the skin layer (Tskin).The temperature in the mixing region is mainly determined by the temperature (Tplasma) of the plasma plume being discharged out of the tube.Moreover, Tmix would determine the temperature of the skin's upper boundary from the view point of the mixing region, so, for simplicity in this case, we denote the upper skin boundary temperature as (Tupper → mix).Similarly, we denote the skin upper boundary temperature from the view point of skin In addition, different mesh sizes were studied and no changes in the convergence efficiency and in the final obtained results were observed.Considering the fluid dynamics of the system, the Navier-Stokes momentum equation is: where ρ is density, P is pressure, u is the velocity and F is given by F = ρg + F st + F v f , ρg is the body force, F st is the surface tension and F vf represents the volume force.Moreover, the continuity equation that is solved alongside with the Navier-Stokes momentum equation is: In order to solve these equations numerically, a no-slip boundary condition is considered at the inner-boundary of the discharge tube and also at the upper surface of the skin that is exposed to the plasma plume; the no-slip boundary condition arises mainly as a result of the adhesion between the fluid and the solid surface.The regime was considered laminar since the largest Reynolds number was ~14, which was far below the number of ~4000 in a turbulent regime.Furthermore, the heat transfer in the mixing region (the domain above the skin surface) is explained using the heat equation, with the use of the temperature-dependent thermophysical properties of air and plasma carrier gas: where C p is the specific heat capacity, ρ represents the density and k is the thermal conductivity.Moreover, the Navier-Stokes momentum equation (Equation ( 8)) and the continuity equation (Equation ( 9)) coupled with the heat transfer (Equation ( 10)) have been used in our previous work [30], but the level-set interface tracking method was used there to explore the plasma evolution in air and water.Furthermore, in order to understand the coupling method used to determine the skin temperature during plasma exposure, one needs to examine (1) the temperature in the mixing region (T mix ) and ( 2) the temperature in the skin layer (T skin ).The temperature in the mixing region is mainly determined by the temperature (T plasma ) of the plasma plume being discharged out of the tube.Moreover, T mix would determine the temperature of the skin's upper boundary from the view point of the mixing region, so, for simplicity in this case, we denote the upper skin boundary temperature as (T upper→mix ).Similarly, we denote the skin upper boundary temperature from the view point of skin layer as (T upper→skin ).As such, the heat transfer in the mixing region is coupled with the skin layer temperature using (T upper ) from these two different viewpoints.The present theory has been implemented in COMSOL MULTIPHYSICS commercial software (COMSOL, Inc.Los Angeles, CA, USA) and the system is solved numerically using the finite element method.The boundary conditions used to solve the PDE system are marked on Figure 3.The material properties used in the present model are summarized in Table 1.
Math.Comput.Appl.2017, 22, 24 6 of 14 USA) and the system is solved numerically using the finite element method.The boundary conditions used to solve the PDE system are marked on Figure 3.The material properties used in the present model are summarized in Table 1.Table 1.Summary of material properties used in the present work.Data for the skin were adopted from [18].

System Geometry and Setup
In order to simulate the system, we have used right circular cylinders to represent the tube, air and the skin layer.The system is considered axisymmetric about the vertical z--axis.The tube dimensions are 1 mm in radius and 5 mm in height with skin radius of 10 mm and thickness of 1.7 mm [18].The cylindrical air domain has a radius of 10 mm and varying heights corresponding to different TTDs.In the present work, we have considered two different heights, which are 1 and 5 mm, to cover different TTDs.The chosen range of TTDs is commensurate with those commonly adopted in experiments.The tube is assumed to have a continuous flow of helium gas.The initial velocity of the carrier gas in the tube has been set at different values which are 100 and 500 mm•s −1 .The flow out of the tube is considered continuous so helium will occupy the discharge zone by interacting with the air and pushing the latter out of the region of interest where the output parameters are scored.In addition, the initial temperature of the system has been set to be 23 °C, while the temperature of the gas flowing out of the nozzle is in turn set to be 30 °C, 50 °C and 70 °C.For a tube with a total length of 20 mm and for an operating voltage of 5 kV for the plasma reactor, the temperature of the plasma being transferred into the tube will be 30 °C.The studies for 50 °C and 70 °C will provide information on the reliance on the plasma temperature of the flow and efficiency of heat transfer to the skin and the surrounding environment.In the current study, a tube section with a length of 15 to 20 mm is simulated, in which the temperature variations are lesser [21].Table 1.Summary of material properties used in the present work.Data for the skin were adopted from [18].

System Geometry and Setup
In order to simulate the system, we have used right circular cylinders to represent the tube, air and the skin layer.The system is considered axisymmetric about the vertical z--axis.The tube dimensions are 1 mm in radius and 5 mm in height with skin radius of 10 mm and thickness of 1.7 mm [18].The cylindrical air domain has a radius of 10 mm and varying heights corresponding to different TTDs.In the present work, we have considered two different heights, which are 1 and 5 mm, to cover different TTDs.The chosen range of TTDs is commensurate with those commonly adopted in experiments.The tube is assumed to have a continuous flow of helium gas.The initial velocity of the carrier gas in the tube has been set at different values which are 100 and 500 mm•s −1 .The flow out of the tube is considered continuous so helium will occupy the discharge zone by interacting with the air and pushing the latter out of the region of interest where the output parameters are scored.In addition, the initial temperature of the system has been set to be 23 • C, while the temperature of the gas flowing out of the nozzle is in turn set to be 30 • C, 50 • C and 70 • C. For a tube with a total length of 20 mm and for an operating voltage of 5 kV for the plasma reactor, the temperature of the plasma being transferred into the tube will be 30 • C. The studies for 50 • C and 70 • C will provide information on the reliance on the plasma temperature of the flow and efficiency of heat transfer to the skin and the surrounding environment.In the current study, a tube section with a length of 15 to 20 mm is simulated, in which the temperature variations are lesser [21].

Benchmarking of the Model
In the present work, the developed theoretical model has been exhaustively benchmarked with a previously performed experiment [21].Moreover, two sets of results are used to accomplish the benchmarking of our theory, which are (1) the temperature variation of the plasma plume versus axial distance away from the outlet and (2) the temperature variation with respect to the gas flow rate out of the tube that is measured at 15 mm distance away from the outlet.Furthermore, a similar experimental setup of Akhlaghi et al. [21] has been modelled.Our benchmarking needed experimental data on temperature variations associated with mixing between the discharged plume and the ambient air.To the best of our understanding, the experimental data employed were the only ones available in the literature.The result the temperature variation versus the axial distance away from the outlet nozzle is shown in Figure 4.The experimental results (see Figure 3 in [21]) for the 6 kV reactor operating voltage are used by obtaining the temperature of the plume at the nozzle outlet through extrapolation.
Math.Comput.Appl.2017, 22, 24 7 of 14 benchmarking of our theory, which are (1) the temperature variation of the plasma plume versus axial distance away from the outlet and ( 2) the temperature variation with respect to the gas flow rate out of the tube that is measured at 15 mm distance away from the outlet.Furthermore, a similar experimental setup of Akhlaghi et al. [21] has been modelled.Our benchmarking needed experimental data on temperature variations associated with mixing between the discharged plume and the ambient air.To the best of our understanding, the experimental data employed were the only ones available in the literature.The result of the temperature variation versus the axial distance away from the outlet nozzle is shown in Figure 4.The experimental results (see Figure 3 in [21]) for the 6 kV reactor operating voltage are used by obtaining the temperature of the plume at the nozzle outlet through extrapolation.The results from the present model are in a good agreement with the previously performed experiments.The average temperature deviation (∆Tmean) was found to be ∼0.347K; this deviation is likely due to the unknown and uncontrollable external factors involved in the experiments.For distances below 10 mm, the temperature predictions are taken from the present model.In addition, the sensitivity of the model in tracking different plasma discharge speeds has been studied by benchmarking the temperature versus gas flow rate at a fixed distance away from the nozzle.The experimental results are also obtained from the previous work of Akhlaghi et al. [21], in which the reactor operating voltage was set to 4 kV.The comparison between the experimental and theoretical results are shown in Figure 5.
In Figure 5, the average temperature deviation (∆Tmean) is ∼0.04 K, which shows excellent agreement between the model predictions and experimental results.The results from the present model are in a good agreement with the previously performed experiments.The average temperature deviation (∆T mean ) was found to be ∼0.347K; this deviation is likely due to the unknown and uncontrollable external factors involved in the experiments.For distances below 10 mm, the temperature predictions are taken from the present model.In addition, the sensitivity of the model in tracking different plasma discharge speeds has been studied by benchmarking the temperature versus gas flow rate at a fixed distance away from the nozzle.The experimental results are also obtained from the previous work of Akhlaghi et al. [21], in which the reactor operating voltage was set to 4 kV.The comparison between the experimental and theoretical results are shown in Figure 5.
In Figure 5, the average temperature deviation (∆T mean ) is ∼0.04 K, which shows excellent agreement between the model predictions and experimental results.

Computed Results
The present computations were performed on a supercomputer with dual Intel Xeon E5-2630 v3 2.40 GHz processors (Intel, Santa Clara, CA, USA) (a total of 16 physical cores hyperthreaded to 32), and the average total computation time for each case presented in this paper was ~6.65 min.Sample results for velocity field streamlines in the discharge zone and temperature distribution in the skin layer for the plasma plume temperature of 70 °C are shown in Figure 6a-d for TTDs of 1 and 5 mm.In general, the temperature of skin under the tube outlet would be higher compared to other sections further away from the outlet.Moreover, the temperature in the skin layer will be higher for plasma plumes with higher velocities.For the larger plasma inlet velocity and TTD, formation of eddies can be observed near the tube outlet (see Figure 6d).These eddies further enhance the heat transfer in the discharge zone due to mixing.

Computed Results
The present computations were performed on a supercomputer with dual Intel Xeon E5-2630 v3 2.40 GHz processors (Intel, Santa Clara, CA, USA) (a total of 16 physical cores hyperthreaded to 32), and the average total computation time for each case presented in this paper was ~6.65 min.Sample results for velocity field streamlines in the discharge zone and temperature distribution in the skin layer for the plasma plume temperature of 70 • C are shown in Figure 6a-d for TTDs of 1 and 5 mm.In general, the temperature of skin under the tube outlet would be higher compared to other sections further away from the outlet.Moreover, the temperature in the skin layer will be higher for plasma plumes with higher velocities.For the larger plasma inlet velocity and TTD, formation of eddies can be observed near the tube outlet (see Figure 6d).These eddies further enhance the heat transfer in the discharge zone due to mixing.

Computed Results
The present computations were performed on a supercomputer with dual Intel Xeon E5-2630 v3 2.40 GHz processors (Intel, Santa Clara, CA, USA) (a total of 16 physical cores hyperthreaded to 32), and the average total computation time for each case presented in this paper was ~6.65 min.Sample results for velocity field streamlines in the discharge zone and temperature distribution in the skin layer for the plasma plume temperature of 70 °C are shown in Figure 6a-d for TTDs of 1 and 5 mm.In general, the temperature of skin under the tube outlet would be higher compared to other sections further away from the outlet.Moreover, the temperature in the skin layer will be higher for plasma plumes with higher velocities.For the larger plasma inlet velocity and TTD, formation of eddies can be observed near the tube outlet (see Figure 6d).These eddies further enhance the heat transfer in the discharge zone due to mixing.The phase field variable shown in Figure 7 indicates the mixing between the plasma plume and air in the discharge zone.The positive and negative Φ values correspond to the presence of air and plasma plume, respectively.Initially, the discharge zone is filled with air (Φ = +1).The fraction of air is reduced when the discharge zone is getting filled by the cold plasma (Φ = −1).The mixing between air and plasma can be quantified through Φ.
The results for the temperature in the discharge zone are shown in Figure 8a-f for TTDs of 1 and 5 mm.In general, the temperature of these two domains increases with time.The small dips in the temperature of the discharge zone are due to the filling up of the zone with the carrier gas.The time when the zone gets filled with helium depends on the velocity of the injection and the TTD, which has a direct influence on the discharge zone volume.Furthermore, for a constant TTD (i.e., 1 mm), the maximum average temperature within this domain would be higher for a larger velocity, which is mainly due to the increased heat transfer between the air and the plasma carrier gas with the flow velocity.Furthermore, when considering a constant velocity and plasma temperature, the configuration with a larger TTD will lead to a larger averaged increase in the discharge zone temperature.Furthermore, as described above, eddies can also be formed for larger TTDs, which can further enhance the heat transfer during continuous injection of plasma plume in the domain, and thus lead to larger averaged increase in the temperature.
(d) TTD = 5 mm, plasma inlet velocity = 500 mm•s −1 .The color bar represents the velocity distribution (in mm•s −1 ) and contour-plot labels represent the temperature of the skin layer (in K).
The phase field variable shown in Figure 7 indicates the mixing between the plasma plume and air in the discharge zone.The positive and negative Φ values correspond to the presence of air and plasma plume, respectively.Initially, the discharge zone is filled with air (Φ = +1).The fraction of air is reduced when the discharge zone is getting filled by the cold plasma (Φ = −1).The mixing between air and plasma can be quantified through Φ.
The results for the temperature in the discharge zone are shown in Figure 8a-f for TTDs of 1 and 5 mm.In general, the temperature of these two domains increases with time.The small dips in the temperature of the discharge zone are due to the filling up of the zone with the carrier gas.The time when the zone gets filled with helium depends on the velocity of the injection and the TTD, which has a direct influence on the discharge zone volume.Furthermore, for a constant TTD (i.e., 1 mm), the maximum average temperature within this domain would be higher for a larger velocity, which is mainly due to the increased heat transfer between the air and the plasma carrier gas with the flow velocity.Furthermore, when considering a constant velocity and plasma temperature, the configuration with a larger TTD will lead to a larger averaged increase in the discharge zone temperature.Furthermore, as described above, eddies can also be formed for larger TTDs, which can further enhance the heat transfer during continuous injection of plasma plume in the domain, and thus lead to larger averaged increase in the temperature.(d) TTD = 5 mm, plasma inlet velocity = 500 mm•s −1 .The color bar represents the velocity distribution (in mm•s −1 ) and contour-plot labels represent the temperature of the skin layer (in K).
The phase field variable shown in Figure 7 indicates the mixing between the plasma plume and air in the discharge zone.The positive and negative Φ values correspond to the presence of air and plasma plume, respectively.Initially, the discharge zone is filled with air (Φ = +1).The fraction of air is reduced when the discharge zone is getting filled by the cold plasma (Φ = −1).The mixing between air and plasma can be quantified through Φ.
The results for the temperature in the discharge zone are shown in Figure 8a-f for TTDs of 1 and 5 mm.In general, the temperature of these two domains increases with time.The small dips in the temperature of the discharge zone are due to the filling up of the zone with the carrier gas.The time when the zone gets filled with helium depends on the velocity of the injection and the TTD, which has a direct influence on the discharge zone volume.Furthermore, for a constant TTD (i.e., 1 mm), the maximum average temperature within this domain would be higher for a larger velocity, which is mainly due to the increased heat transfer between the air and the plasma carrier gas with the flow velocity.Furthermore, when considering a constant velocity and plasma temperature, the configuration with a larger TTD will lead to a larger averaged increase in the discharge zone temperature.Furthermore, as described above, eddies can also be formed for larger TTDs, which can further enhance the heat transfer during continuous injection of plasma plume in the domain, and thus lead to larger averaged increase in the temperature.The relationships among the skin-layer temperature, the plasma temperature, TTD and inlet velocity are shown in Figure 9.A higher plasma temperature will lead to a larger increase in the skin-layer temperature.Moreover, the skin-layer temperature increases more for smaller TTDs and for higher plasma inlet velocities.However, for treatment time up to 100 s, the skin-layer temperature will not increase above 40 °C even with a plasma temperature of 70 °C, showing that skin treatment with cold plasmas under these conditions is relatively safe.If a treatment time longer than 100 s is needed, the temperature curves can be further extrapolated.In comparison, the temperature of the discharge zone shown in Figure 8a-f increases more rapidly when compared to the skin layer due to presence of a gas-gas interface, which further enhances the heat transfer between air and the plasma carrier gas.The present model could determine the plasma effect (in terms of temperature) vs depth in the skin.We scored the temperature of the skin at three different depths (0.5, 1.0, and 1.7 mm) versus the length/radius of the simulated skin domain for the plume temperature of 70 °C at different TTDs of 1 and 5 mm (see Figure 10a,b).The temperatures were lower at larger depths in the skin, while the radial heat dissipation was more significant at smaller depths where the effect of the plume discharge was still dominant.
The licensed cold plasma device known as kINPen MED (neoplas tools GmbH, Greifswald, Germany) was already used in palliative treatment of head and neck tumors [26] and the present temperature model might further help risk assessment of plasma devices for tumor treatment.The relationships among the skin-layer temperature, the plasma temperature, TTD and inlet velocity are shown in Figure 9.A higher plasma temperature will lead to a larger increase in the skin-layer temperature.Moreover, the skin-layer temperature increases more for smaller TTDs and for higher plasma inlet velocities.However, for treatment time up to 100 s, the skin-layer temperature will not increase above 40 • C even with a plasma temperature of 70 • C, showing that skin treatment with cold plasmas under these conditions is relatively safe.If a treatment time longer than 100 s is needed, the temperature curves can be further extrapolated.In comparison, the temperature of the discharge zone shown in Figure 8a-f increases more rapidly when compared to the skin layer due to presence of a gas-gas interface, which further enhances the heat transfer between air and the plasma carrier gas.The present model could determine the plasma effect (in terms of temperature) vs depth in the skin.We scored the temperature of the skin at three different depths (0.5, 1.0, and 1.7 mm) versus the length/radius of the simulated skin domain for the plume temperature of 70 • C at different TTDs of 1 and 5 mm (see Figure 10a,b).The temperatures were lower at larger depths in the skin, while the radial heat dissipation was more significant at smaller depths where the effect of the plume discharge was still dominant.
The licensed cold plasma device known as kINPen MED (neoplas tools GmbH, Greifswald, Germany) was already used in palliative treatment of head and neck tumors [26] and the present temperature model might further help risk assessment of plasma devices for tumor treatment.

Conclusions
In the present work, a complete multiphysics model was built based on finite element analysis to investigate the discharge phenomenon of cold plasma with helium carrier gas ejected out of a tube for skin treatment.In such plasmas with carrier gas, the fractions of plasma constituents are small compared to the carrier gas, so thermofluid analysis is needed for the carrier gas as the major

Conclusions
In the present work, a complete multiphysics model was built based on finite element analysis to investigate the discharge phenomenon of cold plasma with helium carrier gas ejected out of a tube for skin treatment.In such plasmas with carrier gas, the fractions of plasma constituents are small compared to the carrier gas, so thermofluid analysis is needed for the carrier gas as the major

Conclusions
In the present work, a complete multiphysics model was built based on finite element analysis to investigate the discharge phenomenon of cold plasma with helium carrier gas ejected out of a tube for skin treatment.In such plasmas with carrier gas, the fractions of plasma constituents are small compared to the carrier gas, so thermofluid analysis is needed for the carrier gas as the major contributor to the fluid and heat flow.The built model will be a very useful tool for cold plasma medicine applications where thermofluid analysis is necessary.Moreover, the feasibility of increasing the plasma temperature above the conventional temperature limit (30 • C) has been studied; this would be useful in fully exploiting the potential effectiveness and the treatment efficiency of cold plasma devices.The present model can be used to simulate different plasma carrier gases.In the near future, we aim to publicly distribute our computer program together with a full library of different plasma carrier gases.

Figure 1 .
Figure 1.Schematic diagram of simulated sections when plasma is discharged out of the tube.Red arrows: plasma plume; yellow arrows: heat flux vectors in the skin layer.

Figure 1 .
Figure 1.Schematic diagram of simulated sections when plasma is discharged out of the tube.Red arrows: plasma plume; yellow arrows: heat flux vectors in the skin layer.
Math.Comput.Appl.2017, 22,24 5 of 14 the previous work of Yue et al.[29].In the present study, a seamless triangular mesh was employed with the maximum element size of 0.53 mm and the minimum element size of 0.003 mm.The curvature factor of the mesh was set as 0.3 and linear discretization was employed in the present model.The mesh figures for the two configurations with TTDs of 1 and 5 mm are shown in Figure2.

Figure 2 .
Figure 2. Triangular mesh employed for the two different discharge zone sizes.TTD: tube-to-target distance.

Figure 2 .
Figure 2.Triangular mesh employed for the two different discharge zone sizes.TTD: tube-to-target distance.

Figure 3 .
Figure 3. Boundary conditions used in the present work.Blue text: conditions used for fluid dynamics module; Red text: conditions used for heat transfer module.

Figure 3 .
Figure 3. Boundary conditions used in the present work.Blue text: conditions used for fluid dynamics module; Red text: conditions used for heat transfer module.

Figure 4 .
Figure 4. Temperature vs distance from the outlet, at 6 kV operating voltage.

Figure 4 .
Figure 4. Temperature vs distance from the outlet, at 6 kV operating voltage.

Figure 5 .
Figure 5. Temperature variations vs gas flow rate at 4 kV operating voltage, measured at 15 mm away from the outlet.SLPM: standard liter per minute.

Figure 5 .
Figure 5. Temperature variations vs gas flow rate at 4 kV operating voltage, measured at 15 mm away from the outlet.SLPM: standard liter per minute.

of 14 Figure 5 .
Figure 5. Temperature variations vs gas flow rate at 4 kV operating voltage, measured at 15 mm away from the outlet.SLPM: standard liter per minute.

Figure 7 .
Figure 7. Variation of the phase field variable with time in the discharge zone for (a) total simulation time of 100 s and (b) the expanded view for first 5 s.(Φ > 0 represents the presence of air in the zone and Φ < 0 represents the presence of the plasma plume).

Figure 7 .
Figure 7. Variation of the phase field variable with time in the discharge zone for (a) total simulation time of 100 s and (b) the expanded view for first 5 s.(Φ > 0 represents the presence of air in the zone and Φ < 0 represents the presence of the plasma plume).

Figure 7 .Figure 8 .Figure 8 .
Figure 7. Variation of the phase field variable with time in the discharge zone for (a) total simulation time of 100 s and (b) the expanded view for first 5 s.(Φ > 0 represents the presence of air in the zone and Φ < 0 represents the presence of the plasma plume).

Figure 9 .
Figure 9. Variations of the skin-layer temperature with the plasma temperature, TTD and inlet velocity.

Figure 10 .
Figure 10.Plasma effect (in terms of temperature) at different depths in the skin layer for TTD of (a) 1 mm and (b) 5 mm, for plume temperature = 70 °C and velocity = 500 mm•s −1 .Yellow arrows: plasma plume being discharged onto the skin.

Figure 9 .
Figure 9. Variations of the skin-layer temperature with the plasma temperature, TTD and inlet velocity.

14 Figure 9 .
Figure 9. Variations of the skin-layer temperature with the plasma temperature, TTD and inlet velocity.

Figure 10 .
Figure 10.Plasma effect (in terms of temperature) at different depths in the skin layer for TTD of (a) 1 mm and (b) 5 mm, for plume temperature = 70 °C and velocity = 500 mm•s −1 .Yellow arrows: plasma plume being discharged onto the skin.

Figure 10 .
Figure 10.Plasma effect (in terms of temperature) at different depths in the skin layer for TTD of (a) 1 mm and (b) 5 mm, for plume temperature = 70 • C and velocity = 500 mm•s −1 .Yellow arrows: plasma plume being discharged onto the skin.