Chapter 2 Numerical Methods
2.3. Solution Procedure
A general implicit discretized time-marching scheme for the transport equations is employed to solve the discretized equations. It can be written as:
1 1
where the superscripts n and n+1 mean old value (at time t) and new value (at time t+dt) of the variables, respectively. The high order differencing terms and cross diffusion terms are treated using known quantities and retained in the source term and updated explicitly.
In an extended SIMPLE [Chen (1989), Shang et al. (1995), Shang and Chen (1997), and Zhang et al. (2001)] family pressure-correction algorithm, the pressure correction
equation for all-speed flows is formulated using the perturbed equation of state, momentum and continuity equations. The simplified formulations can be written as
p and Du is the pressure-velocity coupling coefficient. Cconsidering
1 1
k n k k k n k n
, the following all-speed
pressure correction equation is obtained,
where the superscript k represents the last iterative value.
The numerical algorithm for solving unsteady, compressible, viscous and heat-conducting gas flows consists the following steps:
1. Initialize the u, v, P, T and other parameters. Set time step t. 2. Set the initial conditions for the calculated time step.
3. Extended SIMPLE algorithm
3.1 Solve the momentum equations at the predictor step.
3.2 Solve the pressure-correction equation.
3.3 Update the velocity, and pressure.
3.4 Solve the secondary pressure-correction equation.
3.5 Update the velocity, and pressure again.
3.6 Solve the energy and species conservation equations.
3.7 Check for convergence of the iteration process. Repeat step 3 until a converged solution is obtained.
4. March to next time step and return to step 2 until a steady-state solution is obtain.
A basic description of the simulation processes is available in Figure 2-2. In addition, parallel computing is implemented and tested on distributed-memory machines using spatial domain decomposition.
Chapter 3
Verifications and Parallel Performance of the Gas Flow Model
In order to validate the feasibility of the developed gas flow model, several simulations from incompressible, low-speed to compressible, high-speed flow and problems with considering conjugate heat transfer were compared with previous simulations wherever possible. The weakly compressible lid-driven cavity flow has been frequently used as the benchmark problem for validation. It was chosen to test the capability of gas flow model in solving the incompressible flow problems. A two-dimensional square cavity with a heat conducting square cylinder at different locations is considered for validating that the gas flow model is suitable for flow problems considering the conjugate heat transfer. The validation of the gas flow model applied in a supersonic gas flow past a confined square in a micro-channel is to demonstrate the ability of the gas flow model in solving high-speed flow problems.
3.1. Lid-driven Cavity Flow
3.1.1. Simulation Conditions
The problem considers the fluid in a two-dimensional square cavity with an upper wall moving in the x-direction at a velocity u as shown in Figure 3-1. The other walls are stationary with the no-slip boundary condition. The pressure and temperature at the boundary are given a neumann boundary condition, meaning that the normal gradients of pressure and temperature are zero. A wide range of Reynolds numbers has been
studied. Solutions are obtained for configurations with non-uniform meshes consisting of 128 by 128 cells.
3.1.2. Validation with Previous Simulations
The streamlines for the driven cavity flow with increasing Reynolds number from 100 to 10,000 along with those by Ghia et al. [1982] are shown in Figure 3-2 and 3-3.
As is well known, the center of primary vortex is offset near the top right corner at Re
=100. It moves towards the geometric center of the cavity with increasing Re. The simulation results are good agreement with those published in the literature. It is clearly that the current gas flow model is capable of reproducing the flow fields as Ghia et al.
[1982] at near-incompressible flow limit in the wide range of Reynolds numbers.
3.2. Low-Speed Gas Flow with Conjugate Heat Transfer
3.2.1. Simulation Conditions
To demonstrate the capability of the solver to simulate gas flow with considering the conjugate heat transfer, we have chosen the conjugate heat transfer problem simulated by Rahman et al. [2008], as shown in Figure 3-4.
The Richardson number ( Ri g
ThT L ui
/ i2 ), a dimensionless parameter represents the importance of natural convection relative to the forced convection, for this investigation is set as 0 to 5, where g is the gravitational acceleration, is the thermal expansion coefficient, L is the length of the square cavity, ui is the inlet velocity, and Th and Ti are the temperature of heated wall and inlet, respectively. L and are assumed as the characteristic length, and velocity, and the corresponding Reynolds number is keeping equal to 100. The inlet width and the square block width are equal to 0.1L and 0.2L, respectively . The solid fluid thermal conductivity ratios f
k
k is equal to 5, where ks and kf are the heat conductivities of the solid and gas, respectively.
The boundary conditions are given as follows:
At the inlet (BC):
At the solid-fluid interfaces of the block:
0, 0, f s ,
3.2.2. Validation with Previous Simulations
Flow over a solid cylinder in a square cavity is considered. Air is used as the working fluid in the cavity. 100 100 uniform computational cells are used for simulations throughout the study. Conjugate heat transfer is considered by solving a steady-state heat conduction equation within the square block and by enforcing the heat flux continuity at the interfaces between gas and solid.
Figure 3-5 to 3-10 show the stream lines and isotherms at Ri 0, Ri 1, and 5
Ri for four configurations of the blocks in the cavity along with the data of Rahman et al. [2008]. The comparisons show good agreement between the two
simulations while the temperature difference (ThTi) is small.
At smaller temperature gradient, the stream lines and isotherms are almost the same as those obtained by Rahman et al.; however, as the temperature difference becomes large, the flow and thermal patterns shown in Figure 3-11 and 3-12 deviate greatly from those by Rahman et al. [2008]. The results indicate that the Boussinesq approximation as often assumed by most of the simulations for mixed convection problems; especially at high Richardson number is highly questionable. For this type of flow, a compressible viscous gas flow model is necessary.
3.3. Micro-scale High-Speed Gas Flow with Slip Boundary Conditions
3.3.1. Simulation Conditions
A 2D compressible laminar flow around a square cylinder with size a (a=1.4 μm) confined in a micro-channel (height Hch 10a, length Lch 50a) is simulated to demonstrate the capability of handling supersonic flow with slip boundary conditions as shown in Figure 3-13.
Velocity-slip and temperature-jump boundary conditions are implemented in the study.
The blockage ratio a Hch is equal to 0.1 and the distance between square and the channel inlet (La) is equal to 5.5a.
The boundary conditions applied in this problem are as follows:
At the inlet:
At the solid-fluid interfaces of the square:
0, 0, f s ,
3.3.2. Validation with Previous Simulations
The distribution of flow properties in the channel is shown in Figure 3-14. The results show that there is a strong bow shock, also called a detached shock, forming in front of the square. Downstream of the bow shock, the pressure, density, and temperature of the flow rise rapidly; the velocity of the fluid drops from "supersonic"
to "subsonic".
The comparison of the normalization distribution of horizontal velocity and temperature obtained by Shterev and Stefanov [2010] (steps
0.00625, 8000 1600
cells), and the present (2000 400 cells) are shown in Figure 3-15 and 3-16. The results obtained that the agreement between gas flow model and SIMPLE-TS data is very good.
Figures 3-17 to 3-20 show the comparisons of horizontal velocity and temperature along the center line of the channel (y H ch 2) with the simulation data [Shterev and Stefanov, 2010] for different spatial steps. The obtained by the gas flow model with
much lower resolution (500 100 - 2000 400 cells) is in a good agreement with those calculated by the SIMPLE-TS method using much higher resolution (2000 400 - 8000 1600 cells). It clearly shows that to obtain a receivable result by SIMPLE-TS method needs to carry out very long calculations compared to those of gas flow model.
3.4. Parallel Performance Study
3.4.1. Test Conditions
The two-dimensional micro-scale supersonic gas flow as shown earlier is used for parallel performance study. 2,000 by 400 uniform computation cells are considered in the study. The parallel performance of gas flow model is investigated by using two type of Krylov subspace method, GMRES and BiCGStab, where the subdomain problems were solved by ILU. The test is operated on the V’ger cluster system (Xeon 3GHz dual-core, dual-CPU) at the Center for Computation Geophysics, National Central University, Taiwan.
3.4.2. Results and Discussion
Figure 3-21 shows the variation of the parallel performance and computation times per time step as a function of the number of processors. The parallel performance is tested by applying GMRES or BiCGStab as the KSP solver and ILU as the sub-domain solver. Result shows speedup is about 46 times as 64 processors are used for the test case with GMRES as the KSP solver. The parallel speedup using GMRES-ILU is slightly better than using BiCGStab-ILU, and the absolute runtime of case using GMRES-ILU is relatively short.
3.5. Summary
In this chapter, the validations of the gas flow model and parallel performance were presented. The validation is found to be very good agreement with previous simulations. Parallel efficiency is reasonably good up to 64 processors with parallel efficiency of ~70%.
Chapter 4
Applications in Low-Temperature Plasma Discharge
This chapter presents the gas flow model coupling with the plasma fluid modeling applied in low-temperature plasma discharge. Figure 4-1 to 4-4 exhibit the plasma properties include electron number density, electron temperature, and ion number densities in an atmospheric-pressure plasma jet w/ and w/o considering the neutral gas flow model. These results demonstrate that influence of gas flow model in the low-temperature plasma discharge should not be underestimated.
4.1. Simulation of Silane/Hydrogen Gas Discharge in a Plasma Enhanced Chemical Vapor Deposition (PECVD) Chamber
This phase demonstrates a large-scale realistic PECVD using the mixture Silane/Hydrogen gas. In order to get a uniform thickness of coating, optimization of the flow streamlines with respect to the influence of deposition parameters by using numerical simulation is necessary. Figure 4-5 show the schematic diagram of the PECVD chamber. Due to the symmetry of the plasma chamber, only half domain is considered in this simulation.
4.1.1. Simulation Conditions
Simulation conditions in this study include: (1) chamber pressure (600 mtorr); (2) a square glass plate (20 Å ~ 20 cm); (3) substrate temperature (250C); (4) gap distance between shower-head and substrate (14 mm) and (5) inflow rate ratio of silane to hydrogen (50 : 80 sccm).
The boundary conditions used in the gas flow model are defined as At the gas inlet:
0, i, i, extrapolation, i
u v v T T P Y Y (4-1)
where v is calculated depending on the gas flow rate, the cross section of the shower i head, and background pressure using a temperature of 300 K.
At the outlet:
The schematic of the computational cells used in this simulation is shown in Figure 4-6.
4.1.2. Results and Discussion
The steady-state neutral flow field including the distributions of species number density, temperature and mass-averaged velocities are simulated by using the developed gas flow model. The applied frequency was 25 kHz and the total gas pressure was 600 mTorr. These properties were then used in the plasma fluid modeling as the background gas properties. Figure 4-7 to 4-8 show the distributions of H2 and SiH4, and gas temperature. In this case, the Reynolds number is small with inflow velocity of silane/hydrogen and gap distance as the characteristic velocity and length, respectively. A small Reynolds number leads to a small Peclet number (Pe Re Pr ), which means that the conduction is dominated. This phenomenon is clearly seen in
Figure 4-7, where the temperature distribution between the electrode and subtract is almost linear. In addition, the density near the substrate surface decreases greatly because of the heated substrate at elevated temperature (250C). Non-uniform background density is important in determining the ionization rate during the simulation. The detailed flow structure is given by the streamlines in Figure 4-9.
4.2. Simulation of a Helium Micro-Cell Plasma
4.2.1. Simulation Conditions
Figure 4-10 shows the schematic diagram of a micro-cell plasma investigated in this study. The micro-cell consists of two ring-shaped electrodes made of aluminum separated by an insulator. The powered electrode is connected to an RF power source (f=13.56 MHz) with amplitude of 300 Volt, and the other one is grounded. Helium gas is applied as the working gas under the atmospheric-pressure condition. This study is numerically solved in a cylindrical coordinate in the right region in Figure 4-10.
40 60 uniform computational cells is used in x- and y-direction, respectively. Table 4 lists the substance properties used in this study. A gravitational field g is considered in the negative y-direction.
The boundary conditions for helium micro-cell plasma include:
At domain boundaries:
4.2.2. Results and Discussion
The present study numerically investigates the characteristics of flow field and heat transfer in a helium micro-cell plasma at atmospheric pressure. Figure 4-11 and 4-12 show the force due to ion-molecule collisions generated by the plasma fluid model for x- and y-momentum equations, respectively. And the energy source duo to ion Joule heating and electron elastic collision is shown in Figure 4-13.
Figure 4-14 shows the temperature distribution without considering the plasma momentum sources. The helium gas is heated and the maximum temperature in the plasma region is about 307 K. The increasing temperature produces an increase in the flow field due to the buoyancy effect. As expected, the buoyancy effect leads a weak clockwise flow field shown in Figure 4-15. However, the situation is totally different when the plasma momentum sources are considered. The clockwise flow field is reversed into a counter-clockwise field shown in Figure 4-16. There is a large component of the force (Fplasma SU2 SV2 ) directed toward the powered electrode 1 2 shown in Figure 4-17. The intensity of the force per unit volume is extremely larger than the buoyancy force, which leads an opposite flow field. The flow is accelerated to a speed of 0.8 m/s. Figure 4-18 shows the distributions of temperature with considering the plasma momentum sources. two-dimensional parallelized plasma fluid modeling code developed by another group
member [Lin, 2010]. The converged steady-state results of flow field, temperature, and other neutral gas properties are present.
4.3.1. Simulation Conditions
The schematic of the helium dielectric barrier discharge atmospheric-pressure plasma jet used in this study is illustrated in Figure 4-20. This system consists of two parallel electrodes made of copper and each of the electrodes is cuboid of 25 50 8 mm3. The lower electrode is connected to an AC power source (f=25 kHz) with amplitude of 250 Volt, and the upper one is grounded. Each electrode is covered with a 35 70 1 mm3 ceramic plate as the dielectric. The gap spacing between the two dielectric plates is kept at 1 mm throughout the study. Helium gas with nitrogen impurity (100ppm) is applied under the atmospheric-pressure condition. The total gas flow rate passing with a cross section of 1 50 mm2 is fixed at 20 slm.
In addition, the initial background gas temperature is assumed to be 300 K, and the surroundings are filled with air (78% N2 and 22% O2). The substrate surface is specified under two kinds of boundary conditions: a) an adiabatic wall (Neumann boundary) and b) an isothermal wall (dirichlet boundary).
The boundary conditions for the helium dielectric barrier atmospheric-pressure plasma jet include:
At the gas inlet:
, 0, , extrapolation,
i i i
u u v T T P Y Y (4-7)
where u is calculated depending on the gas flow rate, the cross section of inlet, and i background pressure using a temperature of 300 K.
At the outlet:
If the gas flows out of the domain,
, , ,
where the subscript “amb” means the state of the surrounding environment.
At the substrate surface:
4.3.2. Results and Discussion
A two-dimensional helium dielectric barrier discharge atmospheric-pressure plasma jet is investigated by solving the governing conservation equations with the boundary conditions. According to the gap between the two dielectric layers, the Reynolds number is estimated to be about 60. The gas flow is assumed to be laminar.
The Knudsen number related to the Mach number and the Reynolds number is less than 0.001, and the gas is treated as continuum. For this reason, the no-slip condition is imposed at the boundary. The 160 160 non-uniform computational cells shown in the Figure 4-21 are employed in the simulation. The parallelized neutral gas flow model uses 8 processors in x-direction and 5 processors in y-direction, respectively.
Figure 4-22 shows the time-history of the velocity and temperature profiles at the export of helium dielectric barrier discharge at atmospheric-pressure with gas flow rate of 20 slm under two kinds of boundary conditions of substrate surface. The computations are continued until it is obtained that the fluid properties have reached a statistically stationary state. It is seen the gas velocity at the exit of the plasma jet increase achieved a maximum velocity of about 14 m/s. The result in an adiabatic substrate surface is quite similar to that in an isothermal substrate surface. The steady-state solutions of the maximum temperatures at the exit of the plasma jet are 321 K and 322.5 K for the adiabatic and isothermal substrate surfaces, respectively.
The two-dimensional spatial distributions of pressure and over-all density are shown in Figure 4-23. The inlet pressure is approximately 760.5 torrs. In Figure 4-24, the spatial distributions of temperature for d 1 mm, H d 10, Re 60 and gas flow rate of 20 slmare shown for two different thermal boundaries. The gas is heating due to ion Joule heating and electron-neutral elastic collision. The maximum shown temperature in the plasma channel is about 323 K.
The velocity components in x- and y-direction are presented in Figure 4-25. Figure 4-26 shows the spatial distributions of the mean-speed, stream lines, and velocity vector.
In the stagnation region, the fluid velocity is zero and the surrounding flow is turning into the wall direction. Results illustrate that there are several pairs of vortices formed in the region between the jet exit and substrate and shown almost symmetric pattern.
Figure 1 shows the steady-state velocity vectors and streamlines at Re 60 , 10
H d , and gas flow rate of 20 slm. The locations of primary and secondary vortexes are depicted in Figure 4-27. A counter-clockwise primary vortex is formed near to the jet exit owing to low-pressure formation near the jet exit. When the momentum of the jet is unable to overcome retarding effect of the primary vortex and
the opposing frictional force of the impingement plate, the clockwise secondary vortex is formed. Figure 4-28 show the numerical simulated stream lines at different time level for d 1 mm, H d 10, Re 60 and gas flow rate of 20 slm. The results denote that the primary and secondary vortices glow and move toward the outlet as time increases. Finally, the flow converges to its steady state.
Figure 4-29 shows the horizontal velocity profiles in the plasma channel at different x position. The entrance length, a length in the channel until the flow velocity profile is fully developed, correlation with the Reynolds Number for laminar flow can be expressed as 0.06Re. It is clear that the flow is nearly fully developed with a parabolic velocity profile. The temperature profile in the plasma channel at various x positions is illustrated in Figure 4-30. The temperature difference between an adiabatic and an isothermal substrate surface is about 3 K.
Figure 4-31 shows the horizontal velocity and temperature profiles along the center line of the plasma jet. The gas flow velocity increases with increasing the distance from the inlet, and reaches a nearly constant value about 11.8 m/s in the plasma channel. A pronounced reduction in the velocity of the helium gas can be observed out of the plasma channel. The gas temperature in the plasma region is continuously increasing due to the electron-neutral elastic collision and ion Joule heating. The maximum values in the plasma channel, respectively, are about 324.5 K and 322.8 K for adiabatic and for isothermal boundaries. The gas temperature outside
Figure 4-31 shows the horizontal velocity and temperature profiles along the center line of the plasma jet. The gas flow velocity increases with increasing the distance from the inlet, and reaches a nearly constant value about 11.8 m/s in the plasma channel. A pronounced reduction in the velocity of the helium gas can be observed out of the plasma channel. The gas temperature in the plasma region is continuously increasing due to the electron-neutral elastic collision and ion Joule heating. The maximum values in the plasma channel, respectively, are about 324.5 K and 322.8 K for adiabatic and for isothermal boundaries. The gas temperature outside