Numerical simulation of pulsed planar evaporation into background gas based on direct Monte Carlo simulation and solution of the BGK model kinetic equation

A numerical study of the planar gas expansion under pulsed evaporation into the background gas is carried out. The chosen conditions are typical for nanosecond laser deposition of thin films and nanostructure synthesis, with the saturated gas pressure at the surface of 5.4 MPa and the background pressure of 50 and 500 Pa. The problem is solved based on the direct simulation Monte Carlo method and direct numerical solution of the BGK model kinetic equation. A generally good agreement was obtained for all computed macroscopic quantities, with the exception of the higher density peak in the compressed layer and a wider shock front in the background gas for the BGK model.


Introduction
Nanosecond pulsed laser ablation in a background gas is widely used for deposition of thin multicomponent films, synthesis of nanoparticles, elemental analysis, processing, cleaning, and surface structuring. Investigation of the laser-induced plume dynamics are important for optimizing the processes related to the various applications. The problem has been the subject of numerous theoretical and numerical studies. Analytical models of gas expansion into the background gas have been developed based on the conservation laws to estimate the limiting plume length [1] and calculate the temporal evolution of the contact surface and internal and external shock waves [2,3]. A blast wave model and an empirical drag force model were shown to well predict shock front evolution [4].
The first numerical studies were one-dimensional. The oscillatory behaviour of plume expansion dynamics was shown on the basis of the solution of Euler equations in spherical coordinates [5]. Planar direct simulation Monte Carlo (DSMC) calculations have been performed for analysis of plume oscillations [6,7] and the dynamics of velocity splitting of ablated particles [8]. Later different twodimensional fluid models have been applied for analysis of laser-induced plume expansion for the continuum regime [9]. Influence of the ambient atmosphere on the initial stage of vapor expansion was investigated on the basis of solutions of the gas dynamics equations [10]. Vortex structures near the evaporation surface were investigated by solving the system of Navier-Stokes equations for a gas mixture [11]. A combined approach combining the DSMC calculation of the initial stage of the plume expansion and the method of test particles for further diffusion of evaporated particles through the background gas was employed to study the ablation of a two-component substance [12] and to analyze the deposition rate and stoichiometry of the deposited films [13]. A combined approach combining the method of large particles for the initial stage and the DSMC for the subsequent expansion stage, taking into account the collective motion of the laser-ablated species and the background-gas particles, was applied to analyze the expansion dynamics for the background pressure of 0.7 to 70 Pa [14]. Recently, the entire flow was calculated by the DSMC method for the case of atmospheric background pressure [15,16].
In spite of the large number of works in this area, there has been no systematic study of the process. Numerical studies were carried out for various special cases without mutual verification, and most of the used numerical approaches were approximate. To accurately solve this problem, we propose to use different kinetic approaches: DSMC method and the solution of the Boltzmann kinetic equation with a model collision integral. Previously, these approaches were employed to calculate the pulsed evaporation into vacuum [17]. The problem is solved in a one-dimensional planar formulation, which corresponds to the initial stage of the expansion of the laser plume.
The main goal of the analysis is to cross-verify the physical models as well as numerical methods implemented in our computer codes. Computationally the problem is very demanding as it involves enormous gradients of all calculated quantities due to the huge difference between parameters of the background and evaporated gas as well as rapidly changing degrees of rarefaction from nearly continuum to nearly free-molecular flow. Moreover, the numerical method to solve the kinetic equation must be able to handle discontinuities of the distribution function propagating from x = 0.

Formulation of the problem
A one-dimensional planar problem of pulsed evaporation of molecules into a background gas is considered. The laser-induced plume is proposed to be neutral and the Clausius-Clapeyron equation is used to describe relation between the surface temperature and the saturated gas pressure (in other words, the classical mechanism of evaporation is assumed).
Molecules are evaporated with the energy corresponding to a surface temperature T0. It is assumed that during time interval  particle flux  is constant and equal to is the most probable thermal speed, k is the Boltzmann constant, m is the molecular mass. The molecules are treated in the frames of the hard sphere model. During evaporation, the condition of complete absorption of backscattering particles is set on the evaporation surface. After completion of evaporation, two variants of the boundary condition are considered: the absorption of all returning particles or their specular reflection.
The initial parameters were set in such a way as to correspond to the typical conditions of pulsed laser deposition. The temperature of the evaporating surface is set to T0 = 5000 K, the temperature of the background gas is Tb = 300 K. The duration of evaporation is  = 10 ns. For evaporated and background particles, we set the same values of the mass m and the diameter of the collision cross section d, corresponding to argon (m = 40 a.m.u., d = 3.7 Å). Assuming that the total number of evaporated particles is 10 15 and evaporation occurs from a spot with a radius of 1 mm, the saturated vapor pressure on the surface is p0 = 5.4 MPa and density n0 = 7.827·10 25 m -3 . These conditions correspond to the evaporation of 34 monolayers of the matter. The background gas pressure was set to 50 Pa (which corresponds to typical conditions for film deposition) and 500 Pa (which corresponds to the conditions for the synthesis of nanoparticles). The ratio of saturated vapor density to background density is 6480 and 648, respectively. For such conditions, the gas flow is characterized by a sufficiently large scatter in the rarefaction degree. The local Knudsen number changes in space from 0.01 to 10, while the overall Knudsen number can be estimated as 0.01.

Description of kinetic approaches
The traditional DSMC scheme [18] is used. The elaborated code has been previously applied for analysis of pulsed expansion under evaporation into vacuum [17,19,20]. The physical domain is divided into computational cells, and the solution is advanced in discrete time steps. The modeling process is divided into two parts: first, the collisionless movement of particles in space and, second, the simulation of collisions between particles. The temporal evolution of the gas cloud is constructed by following the model molecules. The interparticle collisions are simulated in accordance with the "no-time-counter" scheme. The more detail description can be found in [17]. The deterministic calculation is based on the solution of the Bhatnagar-Gross-Krook (BGK) model kinetic equation [21] by means of a high-order conservative discrete-velocity method developed by the second author in a sequence of works [22][23][24] and implemented into "Nesvetay" parallel code. The code solves the kinetic equation in the general three-dimensional formulation; the one-dimensional calculations can be carried out using a special arrangement of the spatial mesh and making some modifications to the numerical flux computation procedure. To improve computational efficiency and reduce the required computational time for the present study the numerical method was extended to deforming spatial meshes based on the arbitrary Lagrangian-Eulerian (ALE) approach from [25]. Complete details on the construction of the ALE method and its verification will be reported in a separate publication.
The overall calculation procedure is as follows. At t = 0 the spatial domain of 0.1 mm length is non-uniformly divided in 400 spatial cells so that the cell size grows from 1.25·10 -4 mm near the surface to 4.4·10 -4 mm at the right boundary. The three-dimensional velocity mesh contains 140x100x100 cells for pb = 50 Pa and 50x36x36 cells for pb = 500 Pa. The initial stage of the flow t ≤  is computed on the fixed mesh. Then, the movement of the spatial mesh is prescribed from for t >  so that the spatial domain is extended in the positive direction. At the final flow time t = 100 the spatial domain extends up to x = 4.15 mm with the cell size varying from 0.005 mm near the surface up to 0.0183 mm near the right boundary.
Combination of coarse spatial cells and generally larger cell sizes during the latter stages of the flow results in an order of magnitude reduction of the required computational time and computer memory as compared to the conventional method on a fixed spatial mesh.

Results and discussion
The calculations have been carried out up to time t = 100 using both approaches. The initial stage of the process is qualitatively similar to the one-dimensional studies reported in the earlier works [26][27][28].
Here for the purpose of the comparison of the approaches we are more interested in later moments of the flow once the evaporation is stopped. Figures 1 and 2 show the profiles of density, velocity, temperature, pressure, and the Mach number during evaporation into the background gas for time moments t = 10, 30, 50 and 100 . For the output times under consideration, almost the entire region occupied by evaporated particles is supersonic. It is seen that a shock wave is formed which propagates into the background gas. Its structure is more visible for the larger of the two considered background pressures, which is expected. Behind it a complicated flow pattern is formed with a wellpronounced peak in density, which decays rapidly with time.
Overall, a reasonably good agreement between DSMC and BGK solutions is observed. Firstly, positions of all waves match quite well. Secondly, the amplitudes of all quantities agree perfectly except density, for which the BGK solution predicts somewhat high peak values. The propagating shock wave structure is different for the lower of two pressures (50 Pa), which is expected due to the use of the molecular-velocity-independent collision frequency in the BGK model. However, for 500 Pa the differences on the shock wave front are almost invisible.
We also performed calculations with specular reflection of particles from the evaporation surface. Figure 3 depicts comparison with the complete condensation boundary condition for t = 50  and background pressure 500 Pa. One can see that the positions and amplitude of all waves agree quite well. The differences occur near the boundary x = 0 and have little effect on the external flow pattern.  Figure 1. Profiles of density, velocity, temperature, and pressure at times of 10, 30, 50, 100  for the background pressure 50 Pa (left) and 500 Pa (right) and complete condensation of particles at the evaporation surface.   During time  < t < 50, only 12% of the evaporated particles recondense on the surface, which manifests itself in the density decrease near the surface. This value of the post-pulse particle backflow is in good agreement with the known data on the back flux for evaporation into vacuum [29].

Conclusion
We have studied numerically one-dimensional gas expansion under pulsed evaporation into background gas for two different background pressures. The considered flow is time-dependent with varying rarefaction regimes, which makes it computationally difficult. The results of two kinetic approaches (DSMC and BGK model kinetic equation) are compared against each other. A generally good agreement for all computed macroscopic quantities is shown. This allows us to conclude that both approaches are well suited for the present physical studies. Further work will include solving the problem in the axisymmetric formulation.