Theoretical analysis

Equations (1.2-9) and (1.2-10) allow one to perform analysis of flow through the orifice (Lch/hch=0) connecting two tanks filled by gas (Figure 2.3-1). From the engineering point of view, it is very difficult to create infinitely thin rigid wall with one hole. But the theoretical expressions shown below allow one judge whether numerical model is able to reproduce physical processes related to rarefied gas flow or not.

Following two cases will be used to check the correctness of current numerical simulation. The first case: total number of enclosed atoms is constant (N=N1+N2=const);

tanks are of equal size; flow is isothermal T1=T2.

Rate of change of number of gas atoms in left tank (Figure 2.3-1)

1 2 1

Integration of eq. (2.3-1) gives:

26

Constant C1 can be determined using initial condition: N1(t=0)=N10

( )

Number of gas atoms in left tank versus time:

( ) (

¹⁰

)

^{2 1 10 2}

If the flow is isothermal and tanks are of equal size and all gas was initially stored in the left tank, then the eq. (2.3-6) becomes:

( )

Equations (2.3-6) and (2.3-7) allows one to estimate the rate of change of number of gas atoms in tanks. Expressions mentioned above are valid in case of rarefied gas (δ<1),

27

while in case of continuum regime parameters of flow become sensitive to rarefaction degree and shape of orifice as well [15]. Authors of [46] have proposed and equation that relates dimensionless flow rate with pressure drop and rarefaction parameter δ practically for the whole range of gas rarefaction:

2

where the interpolating coefficients are given in Table 2.3-1.

The second case: initially pressures in both tanks were equal; tanks are of equal size; flow is not isothermal T1<T2.

Rate of change of number of gas atoms in left tank:

2 1

Initial number of gas atoms in the left tank: N1(t=0)=N10

1

Number of gas atoms in left tank versus time:

28

Relations stated by eq. (2.3-13) and (2.3-7) allow one to judge whether the numerical solution is valid or not. Similarly to eq. (2.3-7) the last equation has been derived with assumption of highly rarefied flow. Sharipov [44] have discussed non-isothermal flow of gas through thin slit within the wide range of rarefaction, the dimensionless flow rates cause by pressure WP and temperature WT gradients are shown in Table 2.3-2. It should be noted that flow rates have been calculated with assumption of small temperature and pressure gradients, respectively:

2 1 1

In contrast to the relations mentioned above, work [39] deals with rarefied gas flow through long rectangular channel. According to results shown in that paper the mass flow rate can be estimated using following relations in case of pressure and temperature driven flows, respectively:

where the coefficient QP versus δ is given in Table 2.3-3.

( )

¹

( )

^{2 3} 2 1 0

29

where the coefficient QT versus δ is given in Table 2.3-4.

It should be noted that eq. (2.3-15a) and (2.3-15b) with following assumptions:

channel is long (it is not specified “how long is long”); pressure and temperature gradients are small; walls of the channel scatter atoms diffusively. Equations (2.3-15) give very important estimations:

• Mass flow is in direct proportion to: average pressure in the capillary tube;

temperature and pressure gradients along the channel; cross-sectional area of a channel.

• Mass flow is in inverse proportion to length of the tune and average temperature of gas in the channel.

These estimations can be used for further qualitative comparison, described early, analytical and numerical results, obtained in current studying.

30

Chapter 3: Description of mathematical model

3.1 Introduction to molecular dynamic method

Molecular dynamics (MD) is a form of computer simulation method in which atoms and molecules are allowed to interact for a period of time by approximations of known physics, giving a view of the motion of the atoms. Because molecular systems generally consist of a vast number of particles, it is impossible to find the properties of such complex systems analytically; MD simulation solves this problem by using numerical methods. Molecular dynamics is a multidisciplinary method. Its laws and theories stem from mathematics, physics, and chemistry, and it employs algorithms from computer science and information theory. It was originally conceived within theoretical physics in the late 1950s[65].

Molecular dynamics simulations compute the motions of individual molecules in models of solids, liquids, and gases. The key idea here is motion, which describes how positions, velocities, and orientations change with time. In fact, molecular dynamics constitutes a motion picture that follows molecules as they dart to and fro, twisting, turning, colliding with one another, and perhaps, colliding with their container.

This usage is not unique: molecular dynamics may also refer to the motions of real molecules when studied primarily by molecular beam [66] or spectroscopic [67]

technique. Molecular dynamics simulation is the modern realization of an old, essentially old-fashioned, idea in science; namely, the behavior of a system can be computed if one have, for the system’s parts, a set of initial conditions plus forces of interaction.

One should realize that the aim of an MD simulation is not to predict precisely what will happen to a system that has been prepared in a precisely known initial

31

condition. In this respect, MD simulations differ fundamentally from numerical schemes for predicting the trajectory of satellites through space: in the latter case, true trajectory must be predicted. It is impossible to launch an ensemble of satellites and make statistical predictions about their destination. However, in MD simulations, statistical predictions are good enough.

Highly simplified the molecular dynamics simulation algorithm can be described like the simulation which proceeds iteratively by alternatively calculating forces and solving the equations of motion based on the accelerations obtained from the new forces. In practice, almost all MD codes use much more complicated versions of the algorithm, including two steps (predictor and corrector) in solving the equations of motion and many additional steps for e.g. temperature and pressure control, analysis and output.

In common case MD method can be described by the following steps:

1. Give atoms initial positions and choose short time step;

2. Calculate forces and accelerations of atoms;

3. Move atoms;

4. Move time forward;

5. Repeat from the second step as long as you need.

Beginning in theoretical physics, MD method gained popularity in materials science and since 1970s also in biochemistry and biophysics. In chemistry, MD serves as an important tool in protein structure determination and refinement using experimental tools such as X-ray crystallography and NMR. It has also been applied with limited success as a method of refining protein structure predictions. In physics, MD is used to examine dynamics of atomic-level phenomena that cannot be observed

32

directly, such as thin film growth and ion-subplantation. It is also used to examine the physical properties of nano devices that have not or cannot yet be created.

A lot of software exists, both commercial packages and free shareware from universities, for MD simulations. Table (3.1-1) lists some of the available software packages.

Almost all existing MD software packages are highly tailored programs (studying of material properties in nanoscale or biopolymers studying), they also have a lot of limitations (closed program source code, fixed boundary, initial conditions or other conditions).

In this thesis, special MD program code are developed for thermal transpiration effect investigation because such approach gives freedom of choice for all parameters of a model. Fortran programming language was chosen for the program source code. The choice of Fortran language can be explained by properties such like: its simplicity of the command syntax, wide applications in solving engineering problems, and high performance capabilities.

3.2 Initialization

To start MD simulation, initial positions and velocities to all particles in the system must be assigned. As the equilibrium properties of the system do not (or, at least, should not) depend on the choice of initial conditions, all reasonable initial conditions are acceptable. Since one wish to simulate the gas state of a particle model system, it is logical to prepare the system with random atom positions of gas atoms in space. Since a particle cannot approach infinitely close to another, it means that it is necessary to check the distance between particles. Otherwise potential energy will grow to inadequate value.

33

Initial velocities may be randomly assigned or from a previous simulation. In current simulations, in order to save computational time, velocities are taken randomly according to Maxwell velocities distribution for temperature T. For three dimensional case, probability density function for velocity magnitude is described as:

( )

where V - particle’s velocity magnitude; kB - Boltzmann’s constant; m0- mass of the particle.

And the cumulative distribution function is:

( )

^{4 F2} ₂ ^exp

(

^{F 2}

)

₄¹

(

^F

)

Error function was calculated by series expansion:

( ) ( ) ² 2 ²

As a result initial conditions definition procedure can be represented by the following list of steps:

34

1. Define atoms’ initial positions by using random function or arrange them into lattice;

2. Check the distance between particles. If distance does not satisfy required conditions procedure starts from beginning;

3. Define velocity vector magnitude for each particle, by using random number generator and cumulative distribution function;

4. Define angles between velocity vectors and coordinate axes. Angles are uniformly distributed.

3.3 Boundary conditions

Interaction between gas atoms and inner surface of the channel is a very important issue that significantly influences on flow behaviors [34, 68-70]. In description of interaction between gas molecules and solid surface, there are two simple models proposed by Maxwell [19] in 1879: the specular reflection model and the diffuse reflection model. The word 'specular' originated from specularis in Latin which means 'mirror'. The word 'diffuse' originated from diffusus which in Latin means 'spread abroad', 'scatter'.

It should be noted, that both models mentioned above are limiting simplifications of the real gas-solid boundary dealing with atomic structure of solid.

This is so-called realist boundary condition or wall with atomic structure. The last approach allows us to carry out precise investigation of gas flow, but in the same time it supposes use of huge amount of computational power. Consequently it would be beneficial to couple merit of the first group of BC (gas atoms’ parameters after collision with solids can be found analytically) with high precision of the realistic BC. Such

35

combination allowed us to create New statistical Boundary Conditions for argon-tungsten interactions, which will be discussed below.

3.3.1 Specular model

The specular reflection model (Figure 3.3-1a) assumes that the incident molecules reflect on the body surface as the elastic spheres reflect on the entirely elastic surface, i.e., the normal to the surface component of the relative velocity reverses its direction while the parallel to the surface components remain unchanged. Thus the normal pressure originated from the reflected molecules equals to that originated from the incident molecules; the shear stress subjected by the surface from the reflected molecules has the opposite sign to that from the incident molecules and the net shear stress is zero; the total energy exchange with the surface is zero. Implementation of specular reflection is simple. When particle reaches surface, normal velocity component changes sign.

3.3.2 Diffuse model

The diffuse reflection model (Figure 3.3-1b) assumes that molecules leaving surface scatter with an equilibrium Maxwell distribution. The condition of equilibrium is the equality of the surface temperature, the temperature in the Maxwell distribution and the static temperature of the incident flow. In diffuse reflection, velocity of each molecule after reflection is independent of its incident velocity. However, velocities of reflected molecules are randomly distributed within half-range Maxwell distribution.

Equilibrium diffuse reflection requires that both the surface temperature and the temperature associated with the reflected gas according to Maxwellian distribution must be equal to the gas temperature [14].

In current studying diffuse reflection model was used to make sure that created computational model is able to reproduce the theoretical results based on diffusive

36

boundary conditions. Implementation of such model is more sophisticated than in case of specular reflection. Maxwellian cumulative velocity distribution of scattered molecules corresponding to the surface’s temperature Tsurf (Eq. 3.2-2) is used for this reflection. According to the sampling method based on inversion of the cumulative distribution [71], by letting F(V^s) equal to a random function ranf uniformly distributed between 0 and 1, solving the equality, the velocity of scattered atom is obtained. The next step is velocity vector direction definition, for 3D case two angles must be defined (Figure 3.3-2 shows local and global 3D coordinate systems for a scattered particle), first angle α^s (in surface plane) is uniformly distributed between -π and π:

(^ranf) ^{/ 2}

αs =π −π (3.3-1)

The second and angle β^s (in plane perpendicular to the surface) is uniformly distributed between-π and π:

(^ranf) ^{/ 2}

βs =π −π (3.3-2)

The described above reflection model is called Maxwellian type boundary condition. Diffusive model of BC has been extensively used in numerical and analytical studies of rarefied gas flow through micro- channels [15, 47, 48, 72]. Unfortunately, those simple models are fail in case of temperature driven flow [2], moreover, recently it has been shown, that interactions of gas with clean metal surface cannot be precisely described by using diffuse BC[61].

It should be noted that there is a special specular-diffusive type of boundary conditions where the degree of diffusivity or secularity is measured in terms of accommodation coefficient [73-75]. This coefficient shows how much energy and/or momentum gas atom loses while collision with solids. Implementation of

37

accommodation coefficients partially fixes demerits of diffusive model, but this approach requires determination of accommodation coefficients for all flow regimes because their values are found to be in strong dependence with thermodynamic parameters of gas and surface [76, 77].

3.3.3 Wall constructed with atomic structures

The third model of BC that has been used in current study is “Real channel”, i.e.

the channel with walls represented by a set of tungsten atoms. “Real channel” was created by removing of excess atoms of the tungsten bar (atoms of the bar are arranged according to BCC structure of tungsten). This process is illustrated in Figure 3.3-3a.

Each wall of the channel consists of six atomic layers and atoms that comprised the outer layers are held fixed in their equilibrium lattice positions, while the atoms in other layers were permitted to move according to the appropriate classical equations of motion. In order to demonstrate the influence of channel’s irregularities on the flow we have considered two kinds of real channel: with smooth and wavy walls, see left and right pictures, respectively, in the Figure 3.3-3b. Initial velocities of the tungsten atoms were determined from the Maxwell distribution function corresponding to the desired temperature of the channel’s surface Tsurf.

This model of BC is the most accurate and comprehensive one, nevertheless it requires a lot of computational resources and increases the number of simulated molecules significantly. Consequently the time of computations rises as well.

3.3.4 New statistical Boundary Conditions

The energy transfer and other processes accompanying the scattering of rarefied gases from solid surfaces have been the subject of a series of studies. Weinberg and Merrill [78] determined angular distributions for gas atoms scattered by a single-crystal W(110) surface. The experimental results of Janda, Hurst and co-workers [102] allowed

38

the researchers to relate the average kinetic energy of scattered argon atoms to the surface temperature, as well as to the incident kinetic energy. The theoretical explanations for argon atoms scattering from a self-assembled monolayer on Ag(111) have been proposed recently by Fan and Manson [79]. Furthermore, Gibson et al. [80]

conducted a detailed study of Ar scattering from an ordered 1-decanethiol–Au(111) monolayer. Recently, Chase et al. [81] have conducted experimental and molecular dynamics studies of argon scattering from liquid indium. They have shown how the angular and energy distributions of scattered atoms depend on incident energy.

Inapplicability of the simple hard-sphere model (this model assumes that atoms don’t have any potential energy on interaction, i.e. atoms influence on each other while collision only) for the description of gas-surface interactions is also presented.

While there are many publications related to gas-surface interactions, their results are still insufficient to define boundary conditions that can describe the gas flow in micro/nano systems. The best way to study natural processes is to conduct an experiment, but one should realize that accurate measurements of gas-surface interactions on a microscopic level are very expensive and time consuming. Fortunately, the recent achievements in computer science and numerical methods made it possible to investigate such processes using MD simulation method.

In this thesis, MD method was applied to study the argon gas scattering processes on a W(110) surface. This approach made it possible to precisely describe the interaction between argon gas and tungsten. The aims of this part of the work were to study effects of argon scattering on the tungsten surface and to propose boundary conditions describing correlations between the parameters of incident and scattered atoms. The method applied in the present thesis can be simply expressed as the bombardment of a tungsten surface with argon atoms, where further analysis of the

39

scattered atoms’ trajectories was conducted. Analysis of both angular distributions and distributions of velocities of scattered atoms were performed using mean values and root mean square deviations (RMSDs). The combinations of these parameters provide complete information about process of gas atoms scattering process. It is shown that results of current study are in line with experimental and theoretical results obtained by the other researches. The information obtained in simulations was statistically analyzed and represented by polynomial functions of incident energy and angle of incidence. All the functions that state relationship between parameters of impinging gas atoms and scattered atoms have been obtained using the Least Squares Method (LSqM). As a conclusive step of the work, an algorithm describing an implementation of the relations mentioned above will be applied to a real gas flow with tungsten boundary in this study.

These relations can be used to specify boundary conditions for argon-tungsten interactions. This model of boundary condition as well as results of its implementation is discussed in next chapter.

3.4 Potential function and interatomic interaction force

In order to simulate practical molecules, the determination of the suitable potential function is very important. A pair of neutral atoms or molecules is subject to two distinct forces in the limit of large separation and small separation: an attractive force at long ranges (van der Waals force, or dispersion force) and a repulsive force at short ranges (the result of overlapping electron orbitals, referred to as Pauli repulsion from Pauli exclusion principle). The interactions among tungsten atoms were taken as being sums of pairwise Morse potential:

( ) ( )

40

where the potential’s parameters [82] were: DW=0.9906eV, βW=14.116nm^-1, RW=0.3032nm. It should be noted that Morse potential can be used for simulation of metallic structure just in case of pure metal with a temperature significantly below the melting point.

To describe the interaction between both argon-tungsten and argon-argon atoms, Lennard-Jones potential functions are applied and depicted as following:

12 6

where the parameters values used for the case of argon-tungsten and argon-argon, respectively, were: ε^WAr/kB=25.17K, RWAr=2.93Å and ε^Ar/kB=119.18K, RAr=3.4Å [83].

The intermolecular potential of inert monatomic molecules such as Ar is known to be reasonably expressed by this function with the parameters listed above. Moreover, many computational and statistical mechanical studies have been performed with this potential as the model pair potential.

The Lennard-Jones potential is an approximation. The form of the repulsion term has no theoretical justification; the repulsion force should depend exponentially on the distance, but the repulsion term of the L-J formula is more convenient due to the ease and efficiency of computing r¹² as the square of r⁶. Its physical origin is related to Pauli principle: when electronic clouds surrounding atoms start to overlap, energy of the

在文檔中 氬氣原子與鎢原子表面交互作用之新型邊界條件發展及其應用於奈米管流之研究 (頁 42-0)