Multiscale Multilevel Approach to Solution of Nanotechnology Problems

The paper is devoted to a multiscale multilevel approach for the solution of nanotechnology problems on supercomputer systems. The approach uses the combination of continuum mechanics models and the Newton dynamics for individual particles. This combination includes three scale levels: macroscopic, mesoscopic and microscopic. For gas–metal technical systems the following models are used. The quasihydrodynamic system of equations is used as a mathematical model at the macrolevel for gas and solid states. The system of Newton equations is used as a mathematical model at the mesoand microlevels; it is written for nanoparticles of the medium and larger particles moving in the medium. The numerical implementation of the approach is based on the method of splitting into physical processes. The quasihydrodynamic equations are solved by the finite volume method on grids of different types. The Newton equations of motion are solved by Verlet integration in each cell of the grid independently or in groups of connected cells. In the framework of the general methodology, four classes of algorithms and methods of their parallelization are provided. The parallelization uses the principles of geometric parallelism and the efficient partitioning of the computational domain. A special dynamic algorithm is used for load balancing the solvers. The testing of the developed approach was made by the example of the nitrogen outflow from a balloon with high pressure to a vacuum chamber through a micronozzle and a microchannel. The obtained results confirm the high efficiency of the developed methodology.


Introduction
The modern computing allows modeling very large and complex systems and processes.Computer modeling has become one of the most effective tools in many branches of science and production.The present work is devoted to the development and application of computer methods in the field of studying the complex gasdynamic processes in technical micro-and nanosystems.A feature of mathematical problems in this field is the simultaneous study of processes at many scales, including microand nanoscale.One of the modern and actively developing approaches to solving such problems is the multiscale approach which combines the methods of the continuum mechanics (MCM) and the Newton dynamics for individual particles.This allows efficient computer simulation of expensive and difficult to implement physical experiments.e-mail: polyakov@imamod.rue-mail: pvictoria@list.ruIn the present paper, an aspect of modeling technical microsystems is considered, which is related to computing the parameters of gas microflows under technical vacuum conditions.For a correct description of such processes, it is necessary to know in detail the properties of the real gases and reproduce them in a numerical experiment.Many properties of the gases are well studied experimentally in certain ranges of temperatures and pressures and are described in the literature.However, there are properties which can only be predicted theoretically on the basis of the kinetic theory of gases [1].In this case, the obtained theoretical data correspond to a very limited range of temperatures and pressures and can differ substantially from the real properties of the gas.One of the ways to obtain the missing information about the properties of the gas medium in a given range of temperatures and pressures is the molecular dynamic simulation [2,3].
A multiscale two-level approach to calculating the real gas flows in the microchannels of technical systems, applied in a wide range of Knudsen numbers, was presented in [4][5][6].At the molecular level it allows to determine the parameters of the equation of state of a real gas [7], to calculate the kinetic properties of the gas [8], to determine the type of the boundary conditions at the walls of the microchannel [5,6].This approach is based on a combination of MCM and Newton dynamics for individual particles and includes two levels of modeling: macroscopic and microscopic.As a model at the macrolevel, the system of quasigasdynamic (QGD) [9] equations is used, as a model at the microlevel, the method of molecular dynamics (MD) [2,3] is used.The implementation of this approach is based on splitting the problem into physical processes.At the macrolevel of detail the description of multicomponent gas media flows takes place.At the microlevel interactions are calculated for: 1) gas molecules between themselves (forming the equation of state of the mixture, determining the transport coefficients and realizing the mixing of the components); 2) gas molecules and atoms of solid surfaces (describing phenomena in boundary layers).
In what follows, a multiscale multilevel approach to solving nanotechnology problems with the help of supercomputer computing systems is described in some detail.This unification covers three scale levels: macroscopic, mesoscopic and microscopic.To realize it at the mesoscopic level the calculation of the motions of large particles in the gas flow is added, and at the microlevel the interactions of the gas molecules with the particle surface are calculated.The consideration of the mesoscopic level allows to track the detailed behavior of the individual large particles, for example, during their deposition on a substrate in order to create a given spatial nanostructure.

Mathematical formulation
The mathematical model includes three components corresponding to the three scale levels.
At the macroscopic level the QGD equations are used.In the case of a gas mixture, the system of QGD equations is written for each gas separately [9].These equations are written below in the three-dimensional case for a mixture of gases in invariant form with respect to the coordinate system together with the equations of constraints and the equations of state: Here all variables with the index l relate to a gas of the type l, each component has its own concentration n l , mass density ρ l = m l n l (m l is the mass of a molecule of the gas l).Each gas is also characterized by its temperature T l and macroscopic velocity u l .Other parameters of the mixture components: p l are partial pressures of gases in the mixture; E l , H l and ε l are total energy densities, enthalpies, and internal energies of the mixture components; µ l = µ l (T l ), D l = D l (T l ) and χ l = χ l (T l ) are kinetic coefficients of mixture components, namely, coefficients of dynamic viscosity, diffusion and thermal conductivity.
l coincide, up to a sign, with the fluxes of mass density, the momentum density of the corresponding components, and energy density.The exchange terms S (ρu k ) l and S (E) l take into account the redistribution of momentum and energy between the components of the mixture.
The system of equations ( 1)-( 4) is closed by the corresponding initial and boundary conditions.The initial conditions are taken in accordance with the equilibrium state of the gas medium in the absence of interaction with external factors.The boundary conditions on solid surfaces are given in the form of mass density, momentum density and energy density fluxes across the boundary, calculated from empirical formulas or calculated using MD methods [4][5][6].On the free surfaces of the computational domain, the so-called "soft" boundary conditions [9] are enforced.
At the microlevel, the evolution of the system under investigation is described by Newton's equations for a system of particles.The system of the equations of motion for particles of type l, which can be a gas, or a metal, has the following form where i is particle index, l is particle type, N l is total number of particles of type l.A particle of type l with index i has its own mass m l , position vector r l,i = r x,l,i , r y,l,i , r z,l,i , velocity vector v l,i = v x,l,i , v y,l,i , v z,l,i and total force F l,i = F x,l,i , F y,l,i , F z,l,i acting on it.The forces acting on the i-th particle are calculated as a sum of its interactions with the surrounding particles, depending on the potential energy and the component responsible for the external action.The potential energy of the system is represented as a sum of partial energies, the calculation of which takes place according to the formula of the chosen interaction potential.
where U is total potential energy, U ll is interaction potential for particles of type l with particles of type l , F ext l,i is the force of interaction with the external environment.The choice of the interaction potential is based on comparing the mechanical properties of the computer potential model and the real material.
The initial conditions at the microlevel are determined by the equilibrium or quasiequilibrium thermodynamic state of the particle system at a given temperature, pressure, and average momentum.The boundary conditions at the molecular level depend on the simulated situation.
At the mesoscopic level, in the presence of finely dispersed solid impurities in the gas mixture of the same chemical composition, a diffusion approximation is used in which the principal unknown quantity is the particle concentration C. For it, the following nonlinear convection-diffusion equation is used: Here D C is the interdiffusion coefficient of the gas suspension that can be calculated in various ways [10,11].The coefficient µ C describes the mobility of the impurity particles; in the general case it is determined in a complex manner and depends on the shape of the mesoparticles, the temperature and the viscosity of the medium.In this paper the well-known Stokes formula was used to determine the mobility coefficient, which is valid for spherical particles.The term F C is the total force acting on the particles and calculated at the molecular level, u is average velocity of the gas mixture flow.The initial conditions for (7) consist in assigning a constant and, as a rule, a very small concentration of particles in the entire free space.The boundary conditions at solid walls are determined by the chosen interaction model (adhesion, reflection, sliding, etc.).The boundary conditions at the entrance to the medium consist in the assignment of the particle flux.The "soft" boundary conditions are used at the output.
If necessary, to trace the movements of the individual large particles (substantially exceeding the dimensions of the molecules of the gaseous medium), equations similar to (5), and ( 6) are used together with the equation ( 7).

Numerical approach
The numerical approach to the solution is based on the problem splitting into physical processes.In this case, the equations of the quasihydrodynamics (1)-( 3) and the equation ( 7) are solved by the method of finite volumes on grids of various types.Newton dynamics equations are solved according to the Verlet integration [2] either in each grid cell independently or in groups of coupled cells.
The calculation of the macroparameters according to the QGD equations ( 1)-( 3) and the equation ( 7) is carried out over a temporal grid using an explicit numerical algorithm which is based on the finite volume method on grids of arbitrary type.For the convenience of solving the problem in areas of complex geometry, hybrid block meshes consisting of cells of various shapes and sizes can be used.In the two-dimensional case it is proposed to use quadrangles and triangles, in the three-dimensional case -polyhedrons with a number of vertices from four to eight.All parameters of the gas components (density, pressure, temperature, velocity vector components, etc.) refer to the mass centers of the cells.Flux variables refer to the centers of the faces of the cells.Spatial approximations of the basic terms of the QGD equations are performed by standard methods.The chosen computing scheme on time is explicit and two-stage (predictor-corrector).
The system of MD equations is used in additional calculations independently, or as a subgrid algorithm applied within each control volume.At the microscopic level, at every step of the simulation, a system of ordinary differential equations is solved, corresponding to the Newton second law and describing the motion of particles of the molecular dynamic system.To integrate the equations of motion, the Verlet integration [2] is used.
To carry out a correct calculation of the QGD, the model is supplemented by real gas equations of state, transport coefficients and other accompanying parameters (enthalpies, average mean free paths, etc.), as well as real boundary conditions.In the case of a mixture of gases, it is necessary to add the corresponding exchange terms to the equations for the momentum and energy of each component.The calculation of these dependencies, coefficients, and conditions is performed using the MD methods.
The stability condition of the splitting method by physical processes includes standard conditions for the stability of an explicit finite-difference scheme (∆t 1 ≈ h 2 , where ∆t 1 is time step, h is space step) at the macrolevel and the conditions for the stability of the Verlet schemes (∆t To verify the proposed numerical methodology, a verification procedure was used, based on performing calculations with grinding the mesh at the macrolevel, and on performing calculations for systems with different number of particles and increasing the number of particles at the microlevel.

Algorithms and program implementation
The problem modeling on the basis of the multiscale approach under consideration with three levels of detail is carried out with the help of special algorithms that in general, depending on the degree of microlevel use, are divided into four classes.
Algorithms of class 1 are of interest for the study of the properties of gas media and the properties of solid surfaces with which the gas medium contacts in technical applications.As a numerical implementation of the approach in this case, the Velocity Verlet scheme acts.With the help of 1 class algorithms, a database of molecular dynamic calculations (DMDC) is accumulated for the properties of gases and solid materials, which can be used in the framework of other algorithms.
Algorithms of class 2 assume the solution of problems only at macro-and mesolevels based on QGD equations ( 1)-( 3) and the convection-diffusion equation (7).In this case, the properties of the gas mixture components (the equations of state by pressure and energy, the kinetic coefficientsviscosity and thermal conductivity, the exchange terms in the equations for momentum and energy, the parameters of the boundary conditions) and the kinetic coefficients for equation (7) are determined from the above-mentioned DMDC accumulated in advance for the desired temperature and pressure range.
Algorithms of class 3 imply the simultaneous use of equations ( 1)-( 3), equations (7) and Newton equations of mechanics ( 5), (6) in the calculations of QGD.Algorithms of class 3 are realized in the framework of the method of splitting into physical processes.It is assumed here that in the gaseous medium and at its boundaries it is possible to confine ourselves to a local consideration of the processes of interacting the gases of the mixture with each other, with the surfaces of large particles and with solid walls.
Algorithms of class 4 also assume the simultaneous use in the calculations of equations ( 1)-( 3), ( 7) and ( 5), (6).The difference of this case from the previous one is that in some areas of the environment (usually at solid surfaces and in zones of a strong drop of the gas parameters), molecular dynamic calculations are carried out continuously without going to the macrolevel.In the same areas, the principle of locality of molecular interactions is not used, that is, in the general case, the algorithms of class 4 are non-local at the molecular level.
The parallelization is implemented based on the principles of geometric parallelism and rational partitioning of the computational domain.To balance the task loading, a dynamic algorithm is used.The code realization is made using MPI and OpenMP programming.

Calculation results
The code validation is carried out on the example of the problem of the nitrogen outflow from a balloon with high pressure to a vacuum chamber through a micronozzle and a microchannel.The computational geometry is shown in Fig. 1.A symmetric domain is selected, but the expected solution may not be symmetric.
The cylindrical micronozzle has a cross section of rectangular shape, the diameter D 0 ≈ 6.2 µm and length L 0 = 6D 0 ≈ 37.2 µm.The computational domain sizes of the balloon equate D 1 = 6D 0 and L 1 = 10L 0 .The diameter and length of the microchannel equate D 2 = 6D 0 and L 2 = 50L 0 .The nozzle and the microchannel connect the balloon with nitrogen and the open space of the vacuum Copper clusters of a regular cubic form with a linear dimension of 10.83 nm (30 × 30 × 30 elementary cells with the edge length of the FCC lattice 0.361 nm) were considered as large particles.
At the initial moment the gas is at rest: u 1 = u 2 = 0.In this case, in the balloon the gas is under standard normal conditions: T 1 = 295.15K, p 1 = 101325 Pa; in the nozzle, the microchannel and the vacuum chamber, the gas is at the same temperature, but at a lower pressure: The nozzle is blocked on the left by a partition which opens instantly at the beginning of the calculation.The inner surface of the nozzle is considered to be perfectly smooth and thermally insulated.
The calculations were performed using an algorithm of class 2 and were aimed at studying the flow parameters at the initial stage of flow acceleration.They showed that the internal energy of the gas in the balloon is, for a certain period of time, gradually converted into kinetic energy of the gas in the nozzle and in the channel.At the same time, a thermal shock wave passes through the gas medium.Subsequently, the gas temperature returns to the normal and the profiles of all gas macroparameters are aligned.
As an illustration of the above the Fig. 2 shows the distributions of the density, concentration of molecules, temperature, and local Mach number in the gas along the symmetry axis of the computational domain for a number of instants of time.An abbreviation a.u.means arbitrary units, here the spatial coordinates were normalized to the diameter of the micronozzle and the gas macroparameters were normalized to their values under normal conditions.Two-dimensional distributions of the concentration of gas molecules and the velocity modulus are shown in Fig. 3.They show the general picture of the flow near the nozzle and the entrance to the microchannel at the time t = 13.26µs, when the gas acceleration has already stopped, and a regular flow has been established.
In Fig. 4 the trajectories of copper clusters calculated in the case of triggering of single particles are shown.Clusters were placed in the middle of the calculated area of the balloon at zero speed.The flow of gas carried through the nozzle into the chamber accelerated them to a maximum speed equal to the average flow velocity.The calculated trajectories as a whole correspond to theoretical estimates.However, gravity was not taken into account.An even more accurate picture in the future will be obtained when this force is taken into account as well as the surface structure of the clusters and walls of the nozzle and microchannel.

Conclusion
A multiscale multilevel approach to the solution of nanotechnology problems by supercomputer systems is presented.The solution combines continuum mechanics models and Newton dynamics for individual particles at three scale levels: macroscopic, mesoscopic and microscopic.The description of the gas-metal technical systems as a mathematical model at the macrolevel is done using the quasihydrodynamic equations for gas and solid states.As a mathematical model at the meso-and   microlevels, the systems of Newton equations written for nanoparticles of the medium and large particles moving in the medium are used.The numerical implementation of the solution involves the problem splitting into physical processes, the method of grids at the macrolevel and the Verlet integration at the meso-and microlevels.Four classes of algorithms and their parallel implementation are discussed.The implementation uses the principles of geometric parallelism, rational partitioning the computational domain, load balancing of the tasks.A special dynamic algorithm is used for load balancing the solvers.The algorithm validation is done on the example of the problem of the nitrogen outflow containing clusters of copper into vacuum.The obtained preliminary results correspond to theoretical estimates and confirm the high efficiency of the developed methodology.

Figure 3 .
Figure 3. Two-dimensional distributions of gas concentration and velocity modulus near nozzle area at the time moment 13.26 µs

Figure 4 .
Figure 4. Illustration of coarse particles moving: streamlines of particles compressibility coefficients, adiabatic indices, specific heat capacities and individual gas constants of the mixture components (k B is the Boltzmann constant); Pr l , Sc l , Ma l and Re l denote Prandtl, Schmidt, Mach, and Reynolds numbers for the mixture components; λ l are mean free paths; vectors 2 ≈ 1/ max |F 2 |, ∆t 3 ≈ 1/ max |F 3 |, where ∆t 2 , ∆t 3 are time steps, max |F 2 |, max |F 3 | are maximum absolute values of the forces of interaction between the particles at a considered step in time) at the meso-and microlevels.