In this study we have successfully proceeded numerical investigations of noble gas flow through nano- channels with various boundary conditions. Cases of isothermal and non-isothermal flow are considered. As a part of this study we have developed a new statistical boundary conditions and created a parallel program code for implementation on Graphic Processing Units.
Results discussed above allow us to conclude that gaming GPUs have good potential to become a foundation base for “home” or desktop supercomputers. Such GPUs have lower prices than professional systems, but the main disadvantage of gaming devices is the relatively low amount of installed RAM. Another important conclusion is that optimal choice of GPU for computations depends on which algorithms will be implemented on a specific GPU. Proceeded optimizations of CUDA program made it possible to create an effective computational model which allowed us to provide studies of the argon nano-scale flow in a wide range of such parameters like pressure, temperature and Knudsen number.
Rarefied gas flow between two tanks has been investigated in the study and our results are in a good agreement with theoretical and experimental studies of other researches. On the other hand, results of MD simulations have shown that diffuse boundary conditions represent rough rather than the smooth wall. It was found that the known theoretical prediction of flow rate through the channel with diffuse BC is valid only if the length-to-height ration is greater than 10, while in case of short channel, dependence between flow rate and length is nonlinear. Analysis of the obtained results leads us to conclusion that neither rarefaction parameter nor pressure difference
74
influence on the curvature of correlation between flow rate and length-to-height ration in case of pressure driven flow.
The discrepancy of 15% between the values of hydrodynamic conductance measured in simulations of gas flow through the smooth tungsten channel and the conductivity obtained from semi-empirical equation has been observed. However, it was found that in case of rough surface this gap becomes less than 5%. The last fact means that appearance of surface’s roughness turns ideal smooth surface into the more realistic one. On the other we have found hydrodynamic conductance of channel with diffuse BC is close to the conductance of the channel with rough walls. Thus, the last two issues allow one to conclude that flow inside the channel produced with high level of precision cannot be modeled using the assumption of diffuse BC and a more complex model of BC is required.
In this work we developed new boundary conditions which are able to reproduce the behavior of argon atoms collided with tungsten surface. The new proposed boundary conditions are based on Gaussian, uniform and Maxwellian probability distribution functions conducted together by using a set of polynomial functions. Such approach made this model of boundary condition more flexible and realistic compared with the simple specular-diffusive one. In the previous section we have proposed an algorithm to describe how the polynomial and probability distribution functions mentioned above could be used as boundary conditions. Generally, this works explains how to implement normalization and further conversion of the experimental or simulation results into boundary conditions.
We have observed that the distribution of the velocities of the scattered atoms becomes more Maxwell-like in the case of a rough surface, while the mean kinetic
75
energy of scattered atoms is almost independent of surface roughness. Furthermore, it was shown that even small surface’s irregularities transform gas-surface interactions to more diffusive-like.
Application of the new proposed boundary conditions for Ar-W interactions is proven to be able to significantly speed up computations. The computation time by use of this algorithm is only about one-twelfth of computation time of MD simulation involving the real tungsten substrate. This is implemented by replacement of atoms on solid surfaces with surfaces following the gas scattering law described by the new boundary conditions.
76
Tables
Table 2.3-1 Coefficients for interpolation formula (2.3-8).
P2/P1 A B C D E
0 0.182 0.006 0.0011 0.257 0.0020 0.1 0.218 0.006 0.0052 0.235 0.0072 0.5 0.276 0.026 0.0241 0.076 0.0137 0.7 0.297 0.049 0.0253 0.017 0.0090
Table 2.3-2 Reduced flow rates vs δ
δ WP WT δ WP WT
Table 2.3-3 Pressure driven flow coefficient QP vs δ in case of square cross-section
δ QP δ QP δ QP
77
Table 2.3-4 Temperature driven flow coefficient QT vs δ in case of square cross-section
δ QT δ QT δ QT
0 0.4193 0.2 0.3390 5 0.1366 0.001 0.4181 0.4 0.3071 8 0.1017 0.01 0.4110 0.5 0.2953 10 0.0868 0.02 0.4037 0.8 0.2684 15 0.0632 0.04 0.3912 1 0.2545 20 0.0495 0.05 0.3857 1.5 0.2275 30 0.0344 0.08 0.3716 2 0.2070 40 0.0263
0.1 0.3637 4 0.1539
78
Table 3.1- 1 Software packages for molecular dynamic simulations
Name Manufacturer/Creator Description
AMBER Peter Kollman's group at the University of California
Provides a set of programs for applying the AMBER forcefields to simulations of biomolecules. It is written in Fortran
90 and C with support for most major Unix-like systems and compilers.
CHARMM Martin Karplus, Accelrys
Allows generation and analysis of a wide range of molecular simulations.
The most basic kinds of simulation are minimization of a given structure and
production runs of a molecular dynamics trajectory.
NAMD University of Illinois at Urbana-Champaign.
It is a parallel molecular dynamics code designed for high-performance simulation of large biomolecular
systems
Materials
Studio Accelrys
It provides tools for modeling crystal structure and crystallization processes,
for the study of polymer properties, catalysis, and the study of
structure-activity relationships.
TINKER Jay Ponder
It is a computer software application for molecular dynamics simulation with a
complete and general package for molecular mechanics and molecular dynamics, with some special features
for biopolymers.
79
Table 3.5-1 Properties of graphic processing units used in current work [101]
PARAMETER GTS250 GTX580 GTX 460 GTX590
CUDA Cores 128 512 336 1024
Graphics Clock (MHz) 738 772 675 607
Processor Clock (MHz) 1836 1544 1350 1215
Memory Clock (MHz) 1242 2004 1800 1900
Standard Memory Configuration 1 GB GDDR3 Memory Interface Width 256-bit 384-bit 256-bit 384-bit Memory Bandwidth (GB/sec) 70.4 192.4 115.2 327.7
Compute Capability 1.1 2.0 2.1 2.0
Table 3.5-2 Time spent for different procedures (in ms) by CPU (Intel Q9400) and by GPU (NVIDIA GeForce GTS250).
PROCEDURE OpenMP(4 cores) CUDA (GTS 250) Acceleration
Output 700.99 0.07% 656 0.11% 1.07
VerletList 295702 31.36% 0 0.00% N/A
Forces 559109 59.30% 359.9986 0.06% 1553
Integration 87209.1 9.25% 691.99 0.11% 126.03
Temperature 76.99 0.01% 200.999 0.03% 0.383
Exchanges 0 0.00% 611998 99.69% 0
Total Time 942798.1 613906.99
80
Table 3.5-3 Simulation time (in sec) of copper crystal on various GPUs versus number of atoms. System was simulated over 20000 time steps.
Number of atoms
Code without 3rd Newton’s law Code with 3rd Newton’s law
GTX460 GTX580 C2075 GTX460 GTX580 C2075
14896 20.81 19.32 20.75 12.06 8.36 11.79
34461 54.89 52.29 55.59 31.84 30.68 32.54
66326 119.38 121.89 120.89 67.59 76.79 72.90
113491 236.81 222.88 228.23 137.14 133.57 136.46
178956 451.57 410.36 408.16 302.22 256.97 260.29
265721 816.64 679.69 701.08 454.46 414.75 415.74
Table 3.5-4 Normalized time of execution of MD program on various computational systems.
Number of atoms
CPU (intel i7-950) GTX590 Sequential OpenMP
CUDA+OpenMP 1 CPU 4 CPUs
14894 52.20 13.7 1.00
34461 155.78 40.35 2.98
66326 416.34 111.05 5.17
178956 2258.29 590.49 19.67 265721 4703.86 1205.8 33.89
81
Table 4.1-1. Correlation coefficients between <Ei> and <Es> for different angle of incidence βi
βi 0° 5° 10° 15° 20°
Correlation
coefficient 0.999437 0.999431 0.999467 0.999439 0.999435
βi 25° 30° 35° 40° 45°
Correlation
coefficient 0.999452 0.999453 0.999389 0.999429 0.999468
βi 50° 55° 60° 65° 70°
Correlation
coefficient 0.999478 0.999499 0.999486 0.999549 0.999626
Table 4.1-2. Coefficients A0 and B0 of the approximation polynomial (4.1-3)
k 0 1 2 3 4 5
A0k 0.7657 -8.5557*10-5 5.4918*10-5 -4*10-7 - - B0k 0.14 1.183*10-4 3.1536*10-5 -1.963*10-6 4.44*10-8 -3.1069*10-10
Table 4.1-3. Coefficients A1n,k of the approximation polynomial (4.1-4) n
k 0 1 2 3 4 5
0 4.0808E-01 -1.8143E-04 -6.2076E-05 4.0415E-06 -8.7349E-08 6.1889E-10 1 -3.1798E-01 -3.2165E-03 6.3264E-04 -3.1700E-05 5.8999E-07 -3.7420E-09 2 2.2149E-01 5.0261E-03 -9.3310E-04 4.4018E-05 -8.0211E-07 5.0384E-09 3 -9.4382E-02 -2.8516E-03 5.3339E-04 -2.4904E-05 4.5289E-07 -2.8462E-09 4 2.3375E-02 8.1348E-04 -1.5349E-04 7.1564E-06 -1.3024E-07 8.1957E-10 5 -3.2682E-03 -1.2525E-04 2.3666E-05 -1.1032E-06 2.0087E-08 -1.2649E-10 6 2.3874E-04 9.9072E-06 -1.8628E-06 8.6781E-08 -1.5799E-09 9.9487E-12 7 -7.0729E-06 -3.1471E-07 5.8677E-08 -2.7304E-09 4.9678E-11 -3.1269E-13
82
Table 4.1-4 Coefficients A2n,k of the approximation polynomial (4.1-6) k
n 0 1 2 3
0 1.0712E+00 7.3238E-02 -8.8442E-02 5.2060E-02 1 -2.0793E-03 1.0793E-02 -2.0503E-02 1.0220E-02 2 8.6448E-04 -8.6467E-03 1.0127E-02 -5.0983E-03 3 -9.8449E-05 7.6928E-04 -9.0643E-04 4.6986E-04 4 4.5169E-06 -3.1373E-05 3.7164E-05 -1.9772E-05 5 -1.0050E-07 6.5031E-07 -7.7981E-07 4.2601E-07 6 1.0868E-09 -6.6825E-09 8.1503E-09 -4.5652E-09 7 -4.5898E-12 2.7098E-11 -3.3664E-11 1.9270E-11
k
n 4 5 6 7
0 -1.6027E-02 2.6247E-03 -2.1649E-04 7.0685E-06 1 -2.2900E-03 2.5026E-04 -1.2243E-05 1.7791E-07 2 1.3124E-03 -1.8107E-04 1.2747E-05 -3.5960E-07 3 -1.2515E-04 1.7885E-05 -1.3046E-06 3.8140E-08 4 5.4091E-06 -7.9247E-07 5.9136E-08 -1.7658E-09 5 -1.1940E-07 1.7857E-08 -1.3563E-09 4.1131E-11 6 1.3062E-09 -1.9861E-10 1.5291E-11 -4.6908E-13 7 -5.6071E-12 8.6358E-13 -6.7172E-14 2.0781E-15
Table 4.1-5. Coefficients A3n,k of the approximation polynomial (4.1-7) n
k 0 1 2 3 4 5
0 7.6831E-01 -1.0661E-03 1.0522E-04 -3.4014E-06 4.5134E-08 -1.9840E-10 1 -5.5235E-01 -2.6447E-03 3.7658E-05 -7.3489E-06 1.7438E-07 -1.1959E-09 2 3.5233E-01 4.3213E-03 -2.0959E-04 1.4560E-05 -3.0241E-07 1.9698E-09 3 -1.3976E-01 -1.8036E-03 1.1142E-04 -7.7702E-06 1.6259E-07 -1.0674E-09 4 3.3855E-02 2.6667E-04 -2.1592E-05 1.8641E-06 -4.1163E-08 2.7775E-10 5 -4.7690E-03 -1.7706E-06 1.2882E-06 -2.1876E-07 5.3403E-09 -3.7634E-11 6 3.5569E-04 -2.7154E-06 6.5481E-08 1.1619E-08 -3.3836E-10 2.5370E-12 7 -1.0811E-05 1.6567E-07 -7.1913E-09 -1.9147E-10 8.0756E-12 -6.6354E-14
83
Table 4.1-6. Coefficients A4n,k of the approximation polynomial (4.1-8) k
n 0 1 2
0 5.0233E-01 9.3434E-02 1.4669E-03 1 -8.7760E-01 5.1610E+00 -4.5275E-02 7 -9.4203E-01 3.5421E-01 -3.8039E-03 8 1.6723E-01 -5.5725E-02 5.9857E-04 9 -1.9580E-02 5.9177E-03 -6.3514E-05 10 1.4527E-03 -4.0511E-04 4.3416E-06 11 -6.1898E-05 1.6135E-05 -1.7258E-07 12 1.1538E-06 -2.8400E-07 3.0305E-09
Table 4.1-7 Properties of wavy surfaces. A and B are the surface wave amplitude and wavelength, respectively.
Table 4.2-1 Flow rate through the channel versus the length-to-height ratio. Wreal is the flow rate through the real channel with smooth walls; Wdiff is the flow rate through the
channels with diffusive BCs. P2/P1=0.5; T1=T2=350K; δ1=0.2
84
Table 4.2-2 Hydrodynamic conductance of the channel with various BCs on the walls.
Parameters of simulations Lch/hch=0.55; T1=T2=350K; δ1=0.2 Semi Empirical eq. by Sreekanth [31] 0.3372 Specular BC 0.4955
Table 4.2-3 Gas flow rate through the channels with various BCs on the walls, length and pressure ratios.
Table 4.2-4 Flow rate through various channels versus the length-to-height ratio.
Parameters of simulations: P2/P1=0.5; T1=T2=350K; δ1=0.2
85
Table 4.3-1 Gas flow rate through the channels with various BCs on the walls.
Parameters of simulations: P2/P1=1; T2/T1=500/350; δ=0.2
Table 4.3-2 Gas flow rate through the channels with various BCs on the walls.
Parameters of simulations: P2/P1=1; T2/T1=500/350; δ=0.2 Lch/hch Diff. Type A Type D New BC
W W ε. % W ε. % W ε. % 0.55 0.115 0.129 11.96% 0.109 5.27% 0.126 9.19%
1.08 0.090 0.115 28.48% 0.095 5.83% 0.112 24.90%
2.68 0.050 0.072 42.24% 0.055 9.59% 0.07 39.13%
4.80 0.032 0.036 13.03% 0.030 6.29% 0.034 5.97%
9.57 0.017 0.020 15.1% 0.019 8.1% 0.019 8.1%
13.91 0.012 0.014 14% 0.013 6.78% 0.014 14%
86
Figures
Figure 2.1-1 One stage of Knudsen compressor [64]
Figure 2.2-1 Sketch of the system used for simulations. H and L – height and length of the tanks; hch and Lch – height and length of the channel;
87
Figure 2.3-1 Sketch of the 3D system used for theoretical analysis. Volumes of tanks are equal; lateral size of orifice is much smaller than the lateral size of tank.
Figure 3.3-1a Sketch of specular reflection model Left tank Right tank
T1, N
1
T2, N
2
88
Figure 3.3-1b Sketch of diffusive reflection model
Figure 3.3-2. Coordinate system, where αi and βi – incident azimuthal and polar angles, respectively; αS and βS are scattered beam’s azimuthal and polar angles, respectively;
α=αS-αi – azimuthal angle change caused by the scattering of the beam; Vni and Vτi – normal and tangential velocity components of the incident beam; VnS and VτS – normal and tangential velocity components of the scattered beam.
βs
βi
αi
αs
α Vi
Vs
Vτs s
Vn
i
Vn
Vτi
y x
z
αi
89
Figure 3.3-3a. Sketch of tungsten channel used for simulations. The tungsten channel (right picture) was created by removal of excess atoms of the tungsten bar (left picture).
Figure 3.3-3b. Side view of the channel with smooth and wavy walls. AW – wave amplitude BW – is the wavelength; hch – height of the channel.
Figure 3.4-1. Lennard-Jones potential and force functions for argon
Tungsten bar Tungsten channel
smooth channel wavy channel 2B hch
y(x)=AWsin(2πx/BW)
-2 0 2 4 6 8 10
2 3 4 5 6 7 8
r, 10-10m
Lennard-Jones potential
Interatomic interaction force
90
Figure 3.4-2. The Verlet list: a particle i interacts with those particles within the cutoff radius rc the Verlet list contains all particles within a sphere with radius rv-rc.
Figure 3.5-1. Cubic unit of copper.
rv
rc
i
91
Figure 3.5-2 Pseudo code of OpenMP program. N is the number of simulated atoms;
NnMx is the maximum number of neighbors; Nthr is the number of threads sharing the job, while ThrID – id number of current thread; jarlst(NnMx*N) stores indexes of neighboring atoms; iarlst(N) stores the number of neighbors; x(N), Vx(N) and fx(N)
store coordinates, velocities and forces, respectively.
01. Data initialization !initial velocities and positions 02. !$omp parallel private(ThrID,i,j,ii,rx) !begin Verlet lists 03. dev=omp_get_thread_num()
11. jarlst((i-1)*NnMx+iarlst(i))=j 12. !$omp end critical
13. endif
14. enddo
15. enddo
16. !$omp end parallel !end Verlet lists
17. !$omp parallel private(ThrID,i,j,ii,jm,rx,dflj) !begin Forces 18. dev=omp_get_thread_num()
19. do ii=1,N, Nthr
20. i=ii+ ThrID; j=1
21. do while (jarlst((i-1)*NnMx+j).gt.0) 22. jm=jarlst((i-1)*NnMx+j)
34. !$omp end parallel !end Forces
35. !$omp parallel private(ThrID,i,ii) !begin Velocity 36. do ii=1,N, Nthr
37. i=ii+ThrID
38. Vx(i)=Vx(i)+Fx(i)*dt/am0 39. x(i)=x(i)+ Vx(i)*dt 40. enddo
41. !$omp end parallel !end Velocity 42. do k=1,ntau !opening the loop in time 43. call Forces !lines 17-34
44. call Velocity !lines 35-40
45. If (maximum displacement of any atom is greater than 0.5(rc-rn)) then 46. call VerletList !lines 02-16
47. endif
48. if(mod(k,iPrint).eq.0) then !printing of atoms coordinates if “k” is multiple of “iPrint”
49. call Temp(am0,Vx,Temp) !calculation of temperature
50. Print to files (x,Vx,Temp) !printing atom’s coordinates and velocities into files 51. endif
52. enddo !closing the loop in time
92
Figure 3.5-3 Acceleration of MD program versus number of threads. Crystal of 178956 copper atoms was simulated for 20000 time steps. Simulations were executed on Intel
Q9400 quad core CPU.
Figure 3.5-4 Acceleration of “Verlet list” and “Forces” subroutines of MD program versus number of threads. Crystal of 66326 copper atoms was simulated for 1 time step.
Simulations were executed on Intel Q9400 quad core CPU.
0
93
Figure 3.5-5 Pseudo code of CUDA program. N is the number of simulated atoms;
NnMx is the maximum number of neighbors; BlkSz is the number of threads in each block, while N/BlkSz+1 blocks are run on the device; jarlst(NnMx*N) stores indexes of
neighboring atoms; x(N), Vx(N) and fx(N) store coordinates, velocities and forces, respectively.
Figure 3.5-6 Pseudo code of GPU’s subroutine for integration of equation of motion. N is the number of simulated atoms; threadidx%x, blockidx%x and blockdim%x are predefined variables that store information about current thread, current block and size of a block, respectively; x(N), Vx(N) and Fx(N) store coordinates, velocities and forces,
respectively.
01. Data initialization !initial positions and velocities 02. Upload coordinates and velocities to GPU
03. call VerletKernel<<<N/BlkSz+1,BlkSz>>>(N,NnMx,x,jarlst,…) !Verlet lists
04. call ForceKernel<<<N/BlkSz+1,BlkSz>>>(N,NnMx,x,jarlst,Fx) !calculation of forces 05. call VelKernel<<<N/BlkSz+1,BlkSz>>>(N,x,Fx,Vx,…) !integration of equation of motion 06. do k=1,ntau !opening the loop in time
07. Set forces to zero
08. call ForceKernel<<<N/BlkSz+1,BlkSz>>>(N,NnMx,x,jarlst,Fx,…) 09. call VelKernel<<<N/BlkSz+1,BlkSz>>>(N,x,Fx,Vx,…)
10. Download displacements of atoms from GPU
11. If (maximum displacement of any atom is greater than 0.5(rc-rn)) then 12. Set jarlst and displacements to zero
13. call NeighKernel<<<N/BlkSz+1,BlkSz>>>(N,NnMx,x,jarlst,…) 14. endif
15. if(mod(k,iPrint).eq.0) then !printing of atoms coordinates if “k” is multiple of “iPrint”
16. Download velocities and coordinates of atoms from GPU 17. call Temp(am0,Vx,Temp) !calculation of temperature
18. Print to files (x,Vx,Temp) !printing to atom’s coordinates and velocities to files 19. endif
20. enddo !closing the loop in time
01. attributes(global) subroutine VelKernel(N,NnMx,x,Fx,Vx,dl2)
02. i = (blockidx%x-1)*blockdim%x+threadidx%x !computation of ”i” based on grid size 03. if (i.le.N) then
04. Vx(i)=Vx(i)+Fx(i)*dt/am0 !velocity of the i-th atom 05. dx=Vx(i)*dt !instant displacement of the i-th atom 06. x(i)=x(i)+dx ! new coordinate of the i-th atom 07. dlx(i)=dlx(i)+dx !total displacement of the i-th atom 08. endif
09. end subroutine VelKernel
94
Figure 3.5-7a Pseudo code of GPU’s subroutine for construction of neighbors lists.
This algorithm does not allow implementation Newton’s third law. N is the number of simulated atoms; threadidx%x, blockidx%x and blockdim%x are predefined variables
that store information about current thread, current block and size of a block, respectively; jarlst(NnMx*N) array of indexes of neighboring atoms; x(N) stores
coordinates of atoms.
Figure 3.5-7b Pseudo code of GPU’s subroutine for construction of neighbors lists in case of Newton’s third law implementation. N is the number of simulated atoms;
threadidx%x, blockidx%x and blockdim%x are predefined variables that store information about current thread, current block and size of a block, respectively;
jarlst(NnMx*N) array of indexes of neighboring atoms; x(N) stores coordinates of atoms.
01. attributes(global) subroutine VerletKernel(N,NnMx,rnCh2,x,jarlst)
02. i = (blockidx%x-1)*blockdim%x+threadidx%x !computation of ”i” based on grid size 03. if (i.le.N) then
04. njlst=0 05. do j=1,N
06. rx=x(i)-x(j) !distance between atoms i and j
07. rr=rx*rx
08. if (rr.le.rnCh2.and.i.ne.j) then !check if the distance is less than rn
09. njlst=njlst+1 !number of neighbors of i-th atom
02. i = (blockidx%x-1)*blockdim%x+threadidx%x !computation of ”i” based on grid size 03. if (i.le.N-1) then
04. njlst=0 05. do j=i+1,N
06. rx=x(i)-x(j) !distance between atoms i and j
07. rr=rx*rx
08. if (rr.le.rnCh2.and.i.ne.j) then !check if the distance is less than rn
09. njlst=njlst+1 !number of neighbors of i-th atom
95
Figure 3.5-8a Pseudo code of GPU’s subroutine for Forces computations. This algorithm does not use Newton’s third law. N is the number of simulated atoms;
threadidx%x, blockidx%x and blockdim%x are predefined variables that store information about current thread, current block and size of a block, respectively;
jarlst(NnMx*N) array of indexes of neighboring atoms; x(N) and Fx(N) are arrays of coordinates and forces, respectively.
Figure 3.5-8b Pseudo code of GPU’s subroutine for Forces computations. This algorithm uses Newton’s third law. N is the number of simulated atoms; threadidx%x,
blockidx%x and blockdim%x are predefined variables that store information about current thread, current block and size of a block, respectively; jarlst(NnMx*N) array of
indexes of neighboring atoms; x(N) and Fx(N) are arrays of coordinates and forces, respectively.
01. attributes(global) subroutine ForceKernel(N,NnMx,rcCh,X,,jarlst,Fx) 02. i = (blockidx%x-1)*blockdim%x+threadidx%x
03. if (i.le.N) then
04. j=1
05. jm=jarlst((i-1)*NnMx+j) 06. do while (jm.gt.0)
07. rx=x(i)-x(jm)
08. rr=sqrt(rx*rx) !absolute value of distance between atoms 09. dflj=0.E0; dfx=0.0E0;
10. if (rr.gt.rcCh) goto 111 !check if the distance is less than rc
11. dflj=bt*De*(dexp1ij-dexp2ij)/rr !compute acting force using potential
12. dfx=dflj*rx 02. i = (blockidx%x-1)*blockdim%x+threadidx%x
03. if (i.le.N) then
04. j=1
05. jm=jarlst((i-1)*NnMx+j) 06. do while (jm.gt.0)
07. rx=x(i)-x(jm)
08. rr=sqrt(rx*rx) !absolute value of distance between atoms 09. dflj=0.E0; dfx=0.0E0;
10. if (rr.gt.rcCh) goto 111 !check if the distance is less than rc
11. dflj=bt*De*(dexp1ij-dexp2ij)/rr !compute acting force using potential
12. dfx=dflj*rx
13. istat=atomicadd(Fx(i),dfx) !compute the total force acting on atom i 14. istat=atomicadd(Fx(jm),-dfx) !compute the total force acting on atom jm
15. 111 continue
96
Figure 4.1-1 Sketch of argon beam colliding with smooth tungsten surface
Figure 4.1-2 Mean energy of scattered atoms versus energy of impinging atoms at βi=45°. Solid line - result from Janda experiments [102]; black cross
symbols represent MD results obtained in the current study.
Smooth tungsten surface
Incident atoms Scattered atoms
97
Figure 4.1-3 Angular distribution of normalized probability density function in case of Ar scattering from W(110) at TG=295K and incident angle βi=45°. +, Tsurf=500K; ○,
Tsurf=350K.
Figure 4.1-4 Mean energy of scattered atoms versus energy of impinging Ar atoms.
+,βi=0°; ,βi=20°; ○,βi=40°; Δ,βi=60°
98
Figure 4.1-5a Root mean square deviation of scattered atom’s velocities σV over mean velocity of scattered atoms μV versus kinetic energy impinging atoms. +,βi=0°;
,βi=15°; □,βi=30°; ○,βi=45°; Δ,βi=60°; ♦,βi=70°
Figure 4.1-5b Root mean square deviation of scattered atom’s velocities σV over mean velocity of scattered atoms μV versus kinetic energy impinging atoms.
99
Figure 4.1-6 Density of probability of Vs. Tsurf=350K, +, βi=20°, Vi=100m/s; , βi=20°, Vi=1500m/s; solid line – Maxwell distribution; dashed line – normal distribution.
Figure 4.1-7a Normalized root mean square deviation σα of the azimuthal angle αs of scattered atoms versus kinetic energy of impinging atoms and incident polar angle.
+,βi=0°; ,βi=15°; □,βi=30°; ○,βi=45°; Δ,βi=60°; ♦,βi=70°
100
Figure 4.1-7b Normalized root mean square deviation σα of the azimuthal angle αs of scattered atoms versus kinetic energy of impinging atoms and incident polar angle.
Figure 4.1-8a Normalized root mean square deviation σβ of the polar angle βs of scattered atoms versus kinetic energy of impinging atoms and incident polar angle.
+,βi=0°; ,βi=15°; □,βi=30°; ○,βi=45°; Δ,βi=60°; ♦,βi=70°
101
Figure 4.1-8b Normalized root mean square deviation σβ of the polar angle βs of scattered atoms versus kinetic energy of impinging atoms and incident polar angle.
Figure 4.1-9 Fraction of backscattered atoms versus energy of impinging Ar atoms.
+,βi=10°; ,βi=20°; ○,βi=40°; Δ,βi=60°
102
Figure 4.1-10a Average polar angle βs of scattered atoms versus kinetic energy of impinging atoms and incident polar angle. +,βi=0°; ,βi=15°; □,βi=30°; ○,βi=45°;
Δ,βi=60°; ♦,βi=70°
Figure 4.1-10b Average polar angle βs of scattered atoms versus kinetic energy of impinging atoms and incident polar angle.
103
Figure 4.1-11 Correlation coefficients between βi and βs versus energy of incident beam
Figure 4.1-12 Density of probability of βsand αs. Tsurf=350K, βi=20°, Vi=100m/s. +, distribution of αs; , distribution of βs; solid line – normal distribution of βs.
104
Figure 4.1-13 Density of probability of βsand αs. Tsurf=350K, βi=20°, Vi=1500m/s. +, distribution of αs; , distribution of βs; solid line – normal distribution of βs; dashed
line - normal distribution of αs.
Figure 4.1-14 Sketch of wavy surface
105
Figure 4.1-15 Fraction of atoms scattered by tungsten surface in backward direction versus energy of impinging Ar atoms. □,βi=10°; Δ,βi=20°; +,βi=40°; ○,βi=60°. Black symbols correspond to smooth surface. Green, red, blue and yellow symbols correspond
to surface of Type B, C, D and E, respectively.
Figure 4.1-16 Normalized root mean square deviation σα of the azimuthal angle αs of scattered atoms versus kinetic energy of impinging atoms and incident polar angle.
□,βi=10°; Δ,βi=20°; +,βi=40°; ○,βi=60°. Black symbols correspond to smooth surface.
Green, red, blue and yellow symbols correspond to surface of Type B, C, D and E, respectively.
106
Figure 4.1-17 Normalized root mean square deviation σβ of the polar angle βs of scattered atoms versus kinetic energy of impinging atoms and incident polar angle.
□,βi=10°; Δ,βi=20°; +,βi=40°; ○,βi=60°. Black symbols correspond to smooth surface.
Green, red, blue and yellow symbols correspond to surface of Type B, C, D and E, respectively.
Figure 4.1-18a Density of probability of βsand αs. Tsurf=350K, βi=20°, Vi=100m/s. +, distribution of αs (surface Type A); , distribution of βs(surface Type A); solid line –
normal distribution of βs(surface Type A); ○, distribution of αs (surface Type B); ●, distribution of βs(surface Type B); ∆, distribution of αs (surface Type C); ▲, distribution of βs(surface Type C); □, distribution of αs (surface Type D); ■, distribution
of βs(surface Type D); ☆, distribution of αs (surface Type E); ★, distribution of βs(surface Type E).
107
Figure 4.1-18b. Density of probability of βsand αs. Tsurf=350K, βi=20°, Vi=1500m/s. +, distribution of αs (surface Type A); , distribution of βs(surface Type A); solid line –
normal distribution of βs(surface Type A); ○, distribution of αs (surface Type B); ●, distribution of βs(surface Type B); ∆, distribution of αs (surface Type C); ▲, distribution of βs(surface Type C); □, distribution of αs (surface Type D); ■, distribution
of βs(surface Type D); ☆, distribution of αs (surface Type E); ★, distribution of βs(surface Type E).
Figure 4.1-19 Mean energy of scattered atoms versus energy of impinging atoms at βi=45°. □, surface Type A; ○, surface Type B; Δ, surface Type C; ◊, surface Type D; +,
surface Type E.
108
Figure 4.1-20 Root mean square deviation of scattered atom’s velocities σV over mean velocity of scattered atoms μV versus kinetic energy impinging atoms. □,βi=10°;
Δ,βi=20°; +,βi=40°; ○,βi=60°. Black symbols correspond to smooth surface. Green, red, blue and yellow symbols correspond to surface of Type B, C, D and E,
respectively.
Figure 4.2-1 Reduced flow rate versus relative size of the tank.
109
Figure 4.2-2 Density of probability of velocity distribution of gas atoms in tank.
T1=350K, ○, H/hch=1; □, H/hch=4; solid line – Maxwell distribution.
Figure 4.2-3a Normalized flow rate versus channel length-to-height ratio. Parameters of simulations: P2/P1=0.5; diffusive BC. +,δ=0.1 (current study); ◊,δ=0.1 (Sharipov[39]); □,δ=0.2 (current study); ○,δ=0.2 (Sharipov [39]); Δ,δ=0.4 (current
study); ♦,δ=0.4 (Sharipov [39]).
110
Figure 4.2-3b Normalized flow rate versus channel length-to-height ratio in a log-log scale. Parameters of simulations: P2/P1=0.5; diffusive BC. +,δ=0.1 (current study);
◊,δ=0.1 (Sharipov [39]); □,δ=0.2 (current study); ○,δ=0.2 (Sharipov [39]); Δ,δ=0.4 (current study); ♦,δ=0.4 (Sharipov [39]).
Figure 4.2-4a Normalized flow rate versus channel length-to-height ratio. Parameters of simulations: P2/P1=0.5 (black symbols); P2/P1=0.9 (red symbols); P2/P1=0.1 (blue
symbols); diffusive BC. □,δ=0.2 (current study); ○,δ=0.2 (Sharipov [39]).
111
Figure 4.2-4b Normalized flow rate versus channel length-to-height ratio in log-log scale. Parameters of simulations: P2/P1=0.5 (black symbols); P2/P1=0.9 (red symbols);
P2/P1=0.1 (blue symbols); diffusive BC δ=0.2;. □, current study; ○, Sharipov [39].
Figure 4.2-5 Normalized flow rate versus pressure difference. ■, Real Channel (Type A); ◊, Real Channel (Type B); □, Real Channel (Type C); +, Real Channel (Type E); Δ,
Diffuse BC; ●, Semi empirical eq. (1.4-1) [31]; ×, Specular BC.
F
113
Figure 4.3-2b Normalized flow rate versus channel length-to-height ratio in log-log scale. Parameters of simulations: P2/P1=1; T2/T1=500/350; diffusive BC.
+,δ=0.2 - current study; +,δ=0.2 – eq. (2.3-15b) derived by Sharipov [39];
◊,δ=0.1 - current study; ◊,δ=0.1 - eq. (2.3-15b) derived by Sharipov [39];
●,δ=0.4 - current study; ●,δ=0.4 - eq. (2.3-15b) derived by Sharipov [39].
114
References
[1] I.A. Waitz, G. Gauba, Y.-S. Tzeng, Combustors for Micro-Gas Turbine Engines, Journal of Fluids Engineering, 120 (1998) 109-117.
[2] F. Sharipov, V. Seleznev, Data on Internal Rarefied Gas Flows, Journal of Physical and Chemical Reference Data, 27 (1998) 657-706.
[3] R.M. Young, Analysis of a micromachine based vacuum pump on a chip actuated by the thermal transpiration effect, Journal of Vacuum Science & Technology B:
Microelectronics and Nanometer Structures, 17 (1999) 280-287.
[4] H. Yen-Lin, Thermal-Creep-Driven Flows in Knudsen Compressors and Related
[4] H. Yen-Lin, Thermal-Creep-Driven Flows in Knudsen Compressors and Related