+高级检索
网刊加载中。。。

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

Plasma Distribution in RF Plasma Oxidation Process Based on NS-DSMC Hybrid Method and Langmuir Probe Diagnostic  PDF

  • Liu Wansuo 1 ✉
  • Wang Zhibao 2,3
  • Wang Qing 1
  • Lin Zeng 1 ✉
1. School of Mechanical Engineering and Automation, Northeastern University, Shenyang 110819, China； 2. NeuMat (Taian) Surface Technology Co., Ltd, Taian 250111, China； 3. Taian Industrial Technology Research Institute of Vacuum Surface Engineering, Taian 250111, China

Updated：2024-04-23

DOI：10.12442/j.issn.1002-185X.E20230023.

  • Full Text
  • Figs & Tabs
  • References
  • Authors
  • About
CN CITE
OUTLINE

Abstract

The combined method of Navier-Stokes equation and direction simulation Monte-Carlo method was used to simulate the flow in the reaction chamber, and the relationship between gas flow velocity and plasma density was analyzed by the Langmuir probe detection equipment. Results show that the flow uniformity can significantly impact the plasma uniformity. Under the pressure of 6.5 Pa, the plasma density and flow velocity show a positive correlation, which is consistent with the experiment results in atmosphere pressure. This research verifies that improving the flow uniformity can enhance the plasma oxidation uniformity.

Plasma-enhanced chemical vapor deposition (PECVD) is one of the commonly used techniques for vacuum coating. This technique was firstly developed in the late 1970s and has been widely used in the production of consistent low-temperature coatings, particularly in the manufacture of semiconductor films. With the Langmuir probe diagnostic device and hybrid method combining the Navier-Stokes (NS) equation and the direct simulation Monte Carlo (DSMC) method, the plasma state can be diagnosed and analyzed, which ameliorates the plasma distribution to achieve greater uniformity and to improve the growth and quality of the film.

The study of plasma uniformity is typically conducted in the gas environment. Hao et al[

1] studied the effects of chamber radius and chamber height on the plasma parameters and uniformity in the inductively coupled plasma source with frequency of 2 MHz. The pressure ranges from 0.5 Pa to 10 Pa, and the results show that the uniformity of plasma distribution in the chamber gradually deteriorates with increasing the pressure. Gao et al[2] used the Langmuir probe and COMSOL simulation software to investigate the electron density of plasma at different pressures. Significant differences exist in the plasma distribution at various pressures. The experiment and simulation results are in good agreement at pressure of 2 Pa and low radio frequency (RF). Most researches focus on different pressure environments, particularly the relationship between the number of non-flowing gas particles and the amount of plasma. However, in the complex flow process in PECVD chamber, stronger correlation exists between the number of gas particles affected by gas flow velocity and the plasma quantity.

In 1980, Coburn et al[

3] used the spectral diagnostics to correlate gas flow rate with plasma density. Sun et al[4] used the two-dimensional self-consistent model to study the uniformity of plasma in the wafer chamber and demonstrated the effect of strong diffusion on plasma uniformity. Larijani et al[5] used the optical emission spectroscopy to investigate the plasma production through the activation of gas phase CH4+H2 by both hot filaments and plasma discharge in chemical vapor deposition. It is found that the flow rate of methane gas affects the plasma uniformity. Zhang et al[6] described the effect of O2 addition on the spatial uniformity of one-dimensional helium plasma jet arrays. With increasing the O2 flow rate, the addition of excessive O2 may result in the quenching of discharge in the plasma jet array. Liang et al[7] used the two-dimensional self-consistent electrostatic fluid model and found the radial uniformity of plasma in capacitively coupled nitrogen discharge. Although the effect of electrode gap on plasma uniformity was analyzed, the impact of flow rate changes resulting from reduction in electrode gap was neglected.

These results all demonstrate the effect of flow velocity on plasma and the influence of gas environment, but the flow field is rarely investigated. This omission may be attributed to the restrictions of simulation methods. In most cases, the environment of plasma chamber is complex and contains multiple flow regimes, leading to difficulty of current fluid simulation methods. Therefore, it is important to develop an effective method to visualize the relationship between gas flow velocity and plasma density.

In the environments with multiple flow states, such as PECVD, there are significant differences between the calculated results obtained by NS method based on the continuity assumption and the actual values in regions with large Knudsen numbers (Kn). DSMC method is based on the probability and statistics, which can accurately simulate the flow field, but it can only apply to regions with high Kn. Therefore, a hybrid NS-DSMC method which combines the efficiency of NS method and the accuracy of DSMC method is proposed. This method is commonly used to simulate the rarefied flows in outer space and can be applied to calculate complex flow environments in PECVD[

8]. In this research, the relationship between the velocity of the rarefied flow and the plasma density was studied.

1 Experiment

Fig.1 shows the schematic diagram of experiment device: the vacuum chamber had an inner diameter of 340 mm, the height was 350 mm, and the spacing was 60 mm between the upper and lower plates. The upper plate was made of metal mesh. The gas was injected into the vacuum chamber through three lines, and the flow was controlled by the mass flow controller. The upper plate was connected to RF power source (13.56 MHz) and generated plasma through capacitive coupling. The experiment parameters were set based on the optimal parameters for the oxidation process. The working vacuum was 6.5 Pa, and RF power was 200 W. The ratio of gas mass flow rate of Ar:O2=1:1, and the Ar/O2 gas flow rate was 5 mL/min.

Fig.1  Schematic diagram of PECVD reaction system

The Langmuir probe used the conventional tungsten probe which was adjusted in position by the drive system, eliminating the sawtooth waveform with the output of ±110 V. Because the probe tip might be easily contaminated and the accuracy of experiment result could be affected, the probe was sequentially cleaned by acetone and alcohol solution before the experiment. The Langmuir probe was located at the position of 20 mm away from the lower plate during the diagnosis. The needle tip was moved outward along the probe, starting from the center of the lower plate. The initial position of the probe was recorded as 0, representing the chamber center. The position 10 indicated the edge of plate, and its distance from initial position was 10 cm.

2 Numerical Method and Validation

2.1 NS method

Based on the open-source software OpenFOAM, a steady-state solver rhoSimpleFoam was used in the continuous flow regime. A three-dimensional model was used for the simulation. The calculation grid adopted the tetrahedral mesh with better adaptability over a wide range and locally used the hexahedral mesh with higher accuracy. The continuous regime simulation method used the density-based compressible flow model[

9]. The fundamental governing equations included the continuity equation, momentum equation, and energy equation, and they were discretized by the finite volume method. The laminar flow model was used in the turbulence model. Thermal radiation was neglected in the calculation. NS method is based on the assumption of continuity. However, with increasing the non-continuum effect, the kinetic theory suggests that NS equation is gradually invalid[10].

2.2 DSMC method

Transition regimes in the microchannel flow field were resolved by DSMC procedure in the dsmcFoam module of OpenFOAM software[

11]. DSMC method is based on the statistical methods instead of the solution of continuity assumption equations. DSMC method describes the system state by tracking the positions and velocities of a large number of particles. Therefore, this method can be used for various flow regimes. Each calculated particle in DSMC represented a specific number of real particles. Through grid division, the motion state of particles in each grid area could be calculated to represent the flow state of that specific area. The molecular free path (λ) was calculated as follows:

λ=25-2ω7-2ω15μρm2πkT1/2 (1)

where ω is the temperature coefficient of viscosity, m is the atomic mass, k is the Boltzmann constant, T is the temperature, μ is the gas dynamic viscosity, and ρ is the gas density. If the mesh size is a fraction of free path[

12], the time-step size is a fraction of the average collision time, and the average number of particles per cell is 25. A variable hard sphere collision model was used[13–14]. The Maxwell wall condition is suitable for rarefied flow interaction, and it was also used in the simulation[15]. The computational efficiency of DSMC method depends on the number of particles, which in turn is related to the free path. Therefore, DSMC method is mainly used to solve the non-continuous state of the flow field.

2.3 NS-DSMC method

NS-DSMC hybrid method was adopted to simulate the microchannel jet flow from transition flow state to continuous flow state in this research. Various techniques were used to solve the problem of algorithmic coupling.

The statistical scatter calculated by DSMC only caused numerical fluctuations of the boundary grids, exerting slight effect on the internal grids. Therefore, the overlapping grid method could be used to transfer the internal data for coupling[

16]. DSMC region was extended into NS region to create the overlapping region[17]. The coupling process must be conducted in the continuous domain for smooth program transition from the particle domain to the continuous domain. The size of the overlapping grid had slight effect on the results[18], which transferred the internally stable solution of DSMC and ensured the validity of NS method in the coupling region. Therefore, DSMC area was extended to NS area, as shown in Fig.2.

Fig.2  Schematic diagram of calculation area division

A state-based coupling approach was proposed[

19–21]. Boundary conditions were imposed on the rarefied gas molecular simulations by Chapman-Enskog velocity distributions[18]. Compared with the boundary conditions of the Maxwellian distribution, the Chapman-Enskog distribution is more suitable for the boundary conditions close to the equilibrium region[16]. The flow velocities of simulated molecules entering DSMC boundary cells were determined by the Chapman-Enskog velocity distribution function based on the macroscopic flow variables[22].

Furthermore, by adopting the sub-relaxation technique, the statistical errors of DSMC can be effectively reduced[

23]. Both sub-relaxation technique and overlapping grid were used to transfer the stable results from DSMC domain to NS domain. The overlapping grid and the Chapman-Enskog distribution were used to exchange boundary conditions between continuous domains and molecular domains[24].

The time decoupling of the steady-state coupling process was achieved by the Schwarz alternation method. By adopting the serial Schwarz alternation method, the exchange of steady-state solutions between subdomains was achieved, promoting the evolution of entire flow field. Using the Schwarz alternation method and Dirichlet boundary conditions, the location of the solver interface for the two domains was pre-specified to facilitate the data exchange. This method is suitable for the coupled low-speed steady-state problems, which uses the implicit iterative coupling method to circumvent the explicit integration to achieve the global steady-state[

18].

NS-DSMC results are consistent with DSMC results, as shown in Fig.3. The DSMC method has high computational accuracy[

25], but it requires a large number of computational resources in the continuous regime. The computational efficiency of coupling method should be primarily considered, rather than the computational accuracy. Therefore, in order to achieve the same level of accuracy with higher computational efficiency, the results of NS-DSMC method are compared with those of DSMC method.

Fig.3  Variation of NS-DSMC results and DSMC results

3 Results and Discussion

3.1 Langmuir probe detection

Under the same discharge parameters, Fig.4 shows the electron density and probe position of the plasma measured by the Langmuir probe. The plasma electron density gradually increases from the plate middle to the plate edge. In the central region of the plate, i.e., at position 0, the electron density of the plasma is slightly higher than that at position 1. The electron density reaches 2.1×1010 cm-3 at position 0, whereas it reaches 2.0×1010 cm-3 at position 1. At position 9, the plasma electron density reaches the maximum value of 2.6×1010 cm-3. However, with moving the probe tip further away from the plate center, the electron density of the plasma is rapidly decreased. When the distance between the probe tip and the plate center exceeds 9 cm, i.e., at position 10, the electron density decreases to 2.3×1010 cm-3. Because the plate has symmetrical structure, the plasma spatial distribution exhibits the symmetrical shape along the axial direction. Eq.(2) shows the general formula for calculation of plasma uniformity[

26–27]:

Uniformity= 1-Max-Min2×Average×100% (2)

Fig.4  Electron density distribution of plasma measured by Langmuir probe

where Max refers to the maximum plasma intensity, Min refers to the minimum plasma intensity, and Average refers to the average plasma intensity.

The device is operated at working pressure of 6.5 Pa and RF power of 200 W to activate the mixed gas of argon and oxygen (Ar:O2=1:1, Ar flow rate=5 mL/min). The plasma uni-formity between the plates is approximately 85%, which is good but still needs further improvement. Usually, the main reasons for the non-uniform distribution of plasma are the standing wave effect, skin effect, edge effect, and the uneven distribution of gas. However, if the skin depth is greater than half of the plate spacing, the skin effect can be ignored[

28]. The distance between the plates in this experiment was 60 mm, and the skin depth was about 100 mm. The standing wave effect normally exists in the process of large-area plasma discharge[29–30]. Using RF power supply with frequency of 13.56 MHz, the RF wavelength is significantly larger than the electrode diameter. Therefore, the impact of the standing wave effect in this experiment can be basically ignored. The effects of standing waves and skin effect can be ignored in this experiment. Therefore, the main reason for the non-uniform distribution of plasma in this experiment may be the edge effect and the uneven gas distribution.

3.2 Flow field simulation results

To further investigate the influence of edge effect and non-uniform gas distribution on the unevenness of plasma distribution, the flow distribution in the axial plane of the probe was discussed. The flow field was simulated by NS-DSMC method throughout the entire chamber. The simulation results show that when the gas enters the vacuum chamber, only a small amount of gas can reach the plate center. More gas molecules exist at the plate edge rather than in the middle. According to Fig.5, the flow characteristics of the gas movement are asymmetric: a large amount of gas exists on the left side, compared with that on the right side. This phenomenon results in the non-uniform distribution of plasma. The plasma density at the plate edge is higher than that at the center position.

Fig.5  Schematic diagram of cross-section of PECVD reaction sys-tem (a); two-dimensional streamline of velocity simulation (b)

The distribution of electron field in PECVD device was analyzed. Fig.6 shows the electric field distribution simulated by ANSYS software. It can be observed that the electric field is concentrated near the plate edge, whereas it is relatively uniform in the central region. The electron density and velocity of the flow field have similar variations, as shown in Fig.7. It can be seen that the flow is basically consistent with the electron density in the central region (area A) with uniform electron field. With moving the position to edge part (area B), the strength of the electric field is increased. The increase in the difference between electron density and flow rate can be explained by the combined effect of flow rate and electron field strength on electron density.

Fig.6  Simulation morphology of electric field

Fig.7  Electron density and velocity of flow field

In the central region, the electric field is uniform, and an approximately linear relationship exists between velocity and electron density. The relationship between plasma density and gas flow velocity is consistent with the results in Ref.[

31]. Yambe et al[31] established the relationship between plasma and flow velocity under atmosphere pressure. The plasma density is logarithmically related to the gas flow velocity. There is no obvious logarithmic relationship under the condition of 6 Pa, but the positive correlation between velocity and electron density is consistent with the results in Ref.[31].

3.3 Structural improvement

To enhance the plasma distribution of PECVD equipment and to improve the surface quality of workpiece, the structure of PECVD device should be optimized by considering the uniformity of gas distribution. The intake device is gradually developed into a showerhead, as shown in Fig.8.

Fig.8  Schematic diagram of ameliorated equipment

The new structure has the original microchannel with 1.2 mm in diameter coupled with the microchannel with 0.8 mm in diameter. The hole depth (microchannel length) is increased to enhance the flow resistance of the showerhead. This design achieves a uniform pressure distribution in the showerhead and minimizes the flow difference between the holes. In addition, the holes are uniformly spread along the horizontal x and y directions to maximize the uniformity of inflow. The results show that although a significant velocity difference exists between the center and the edge, the uniformity of the flow field in the center region is significantly improved.

The simulation results show that the amelioration of the air intake device can improve the gas balance in the vacuum chamber, as shown in Fig.9. Along the x direction, within the range of 0–0.08 m, the flow velocity ranges from 4×10-5 m/s to 4.8×10-5 m/s. The gas distribution between the plates significantly improves after the amelioration.

Fig.9  Flow field simulation of new structure

4 Conclusions

1) Non-uniformity and edge effects can lead to the decrease in the uniformity of plasma distribution.

2) The enhanced electric field caused by edge effects increases the plasma density. Under a uniform electric field in the rarefied environment of 6.5 Pa, the gas flow velocity is positively correlated with the plasma density. The simulation results are consistent with the experimental ones under atmosphere pressure.

3) Hybrid NS-DSMC method can accurately obtain the flow field values and can be effectively used to analyze the plasma uniformity.

4) By modifying the inlet structure, the gas uniformity can be improved, thus enhancing the plasma uniformity. This research provides support for the investigation about the influence of rarefied flow uniformity on plasma distribution, presenting wide application potential.

References

1

Hao Z Y, Hua Y, Song J et al. Physics of Plasmas[J], 2020, 27(4): 043502 [Baidu Scholar] 

2

Gao F, Zhang Y R, Li H et al. Physics of Plasmas[J], 2017, 24(7): 073508 [Baidu Scholar] 

3

Coburn J W, Chen M. Journal of Applied Physics[J], 1980, [Baidu Scholar] 

51(6): 3134 [Baidu Scholar] 

4

Sun X Y, Zhang Y R, Chai S et al. Physics of Plasmas[J], 2019, 26(4): 043503 [Baidu Scholar] 

5

Larijani M M, Le Normand F, Crégut O. Applied Surface Sci-ence[J], 2007, 253(8): 4051 [Baidu Scholar] 

6

Zhang C, Shao T, Zhou Y et al. Applied Physics Letters[J], 2014, 105(4): 044102 [Baidu Scholar] 

7

Liang Y S, Liu Y X, Zhang Y R et al. Journal of Applied Phys-ics[J], 2015, 117(8): 083301 [Baidu Scholar] 

8

Raman S K, Kexin W, Kim T H et al. International Journal of Mechanical Sciences[J], 2020, 176: 105396 [Baidu Scholar] 

9

Gadel-Hak M. J Fluids Eng[J], 1999, 121(1): 5 [Baidu Scholar] 

10

Fan J, Shen C. Journal of Computational Physics[J], 2001, [Baidu Scholar] 

167(2): 393 [Baidu Scholar] 

11

Palharini R C, White C, Scanlon T Jet al. Computers & [Baidu Scholar] 

Fluids[J], 2015, 120: 140 [Baidu Scholar] 

12

Cai G B, Liu L H, He B J et al. Aerospace[J], 2022, 9(11): 706 [Baidu Scholar] 

13

Rafi K M M, Deepu M, Rajesh G. International Journal of Thermal Sciences[J], 2019, 146: 106063 [Baidu Scholar] 

14

Moss J N, Price J M. Journal of Thermophysics and Heat Transfer[J], 1997, 11(3): 321 [Baidu Scholar] 

15

Virgile C, Albert A, Julien L. Acta Astronautica[J], 2022, 195: 295 [Baidu Scholar] 

16

Schwartzentruber T E, Boyd I D. Journal of Computational Physics[J], 2006, 215(2): 402 [Baidu Scholar] 

17

Schwartzentruber T, Scalabrin L, Boyd I. 45th AIAA Aerospace Sciences Meeting and Exhibit[C]. Reno: AIAA, 2007: 613 [Baidu Scholar] 

18

Aktas O, Aluru N R. Journal of Computational Physics[J], 2002, 178(2): 342 [Baidu Scholar] 

19

Roveda R, Goldstein D B, Varghese P L. Journal of Spacecraft and Rockets[J], 1998, 35(3): 258 [Baidu Scholar] 

20

Roveda R, Goldstein D B, Varghese P L. Journal of Spacecraft and Rockets[J], 2000, 37(6): 753 [Baidu Scholar] 

21

Farber K, Farber P, Gräbel J et al. Applied Mathematics and Computation[J], 2016, 272(3): 648 [Baidu Scholar] 

22

Li Z, Dang L, Yang Y et al. AIP Conference Proceedings[J], 2018, 2027(1): 030175 [Baidu Scholar] 

23

Sun Q, Boyd I D. Journal of Thermophysics and Heat Trans- [Baidu Scholar] 

fer[J], 2005, 19(3): 329 [Baidu Scholar] 

24

John B, Damodaran M. IEEE Transactions on Magnetics[J], 2009, 45(11): 4929 [Baidu Scholar] 

25

Wagner W. Journal of Statistical Physics[J], 1992, 66: 1011 [Baidu Scholar] 

26

Kim S S, Chang H Y, Chang C S et al. Applied Physics [Baidu Scholar] 

Letters[J], 2000, 77(4): 492 [Baidu Scholar] 

27

Bliokh Y P, Brodsky Y L, Chashka K B et al. Plasma Sources Science & Technology[J], 2009, 18(1): 015009 [Baidu Scholar] 

28

Kolobov V I, Economou D J. Plasma Sources Science & Technology[J], 1997, 6(2): R1 [Baidu Scholar] 

29

Lee I, Graves D B, Lieberman M A. Plasma Sources Science & Technology[J], 2008, 7(1): 136 [Baidu Scholar] 

30

Ahn S K, Chang H Y. Applied Physics Letters[J], 2008, 93(3): 2091 [Baidu Scholar] 

31

Yambe K, Taka S, Ogura K. Physics of Plasmas[J], 2014, 21(4): 043511 [Baidu Scholar] 

copyright © 2018-2020

总访问量

地址：西安市未央区未央路96号 邮政编码：710016 联系电话：029－86231117

E-mail：rmme@c-nin.com; rmme0626@aliyun.com

版权所有:稀有金属材料与工程 ® 2025 版权所有 技术支持:北京勤云科技发展有限公司 ICP:陕ICP备05006818号-3