Investigation of the Contamination Control in a Cleaning Room with a Moving AGV by 3D Large-Scale Simulation

The motions of the airflow induced by the movement of an automatic guided vehicle (AGV) in a cleanroom are numerically studied by large-scale simulation. For this purpose, numerical experiments scheme based on domain decomposition method is designed. Compared with the related past research, the high Reynolds number is treated by large-scale computation in this work. A domain decomposition Lagrange-Galerkin method is employed to approximate the Navier-Stokes equations and the convection diffusion equation; the stiffnessmatrix is symmetric and an incomplete balancing preconditioned conjugate gradient (PCG)method is employed to solve the linear algebra system iteratively. The end wall effects are readily viewed, and the necessity of the extension to 3 dimensions is confirmed. The effect of the high efficiency particular air (HEPA) filter on contamination control is studied and the proper setting of the speed of the clean air flow is also investigated. More details of the recirculation zones are revealed by the 3D large-scale simulation.


Introduction
Automatic guided vehicles (AGVs) play a more and more important role in the material handling system of modern computer integrated manufacturing systems (CIMS).Based on the demands of just-in-time delivery, AGVs are now used more and more widely in modern hospitals, medical centers, port logistics, airports, and semiconductor industry.With the increase in both size and weight of the wafers, AGVs are used to carry the wafers instead of operators to save labor and to decrease wafer damage and contamination.
Cleanrooms are commonly used to provide a contamination control environment to produce high quality and precision products in modern semiconductor industry.The external clean airflow in the cleanroom plays a role of removing microcontaminants and hazardous gas, which are usually generated in manufacturing and transportation process.A vertical external airflow is usually necessary and efficient to keep the cleanness of the cleanroom; however, some microcontaminants that circulate within recirculation zones and deposit on the products and equipment are extremely difficult to be removed.Therefore, simulation of the motion of the airflow in the cleanroom shows its importance.
However, existing techniques available in the literature show a variety of concerns about the characteristics [1][2][3][4][5][6], system control [7][8][9][10], and path planning [11][12][13][14] of AGVs; there are few reports devoted to the simulation of the motion of the airflow in the cleanroom; see [15][16][17][18][19][20][21].By using an upwind scheme, Kanayama et al. analyzed a twodimensional flow around an AGV in a cleanroom [18].The results indicated that the wafer might not be contaminated by the particles from the floor.But unfortunately, the airflow looks rather different from what actually occurs because of the limitation of dimensions and mesh size in computation model; Suh-Jenq Yang proposed a two-dimensional model; without the grating zones to simulation of an AGV moving in cleanroom, only the recirculation zones in a two-dimensional plane were displayed, and the end wall effects were totally neglected.
The current study is to improve the simulation of airflow in the cleanroom and to investigate the contaminant control inside.It is expected to have a better view about the recirculation zones by using large-scale computation based on a domain decomposition method.To handle the problem caused by the nonlinear convective terms of flow problems, which result in the nonsymmetry of the stiffness matrix, an adapted Lagrange-Galerkin method for the domain decomposition method is used.Compared with the classical methods, which employ product-type methods such as GPBiCG or BiCGSTAB as the iteration solver [22], the symmetry of the matrix enables preconditioned conjugate gradient (PCG) method to be employed to solve the interface problem of the domain decomposition system.By using an incomplete balancing domain decomposition preconditioner, problems with up to 30 million degrees of freedom (DOF) can be solved [23] on a cluster with 20 Intel(R) Core(TM) i7 920 processors.In order to investigate the effects of end walls, two types of boundary conditions are considered.The necessity of the vertical external airflow, which plays an important role in removing the microcontaminants, is also to be proved by numerical experiments in this work.
The remaining sections are arranged as follows.Section 2 gives a brief description about the physical model.Section 3 describes the formulations, scheme, and the domain decomposition implementation.Numerical results and discussions are present in Section 4, and finally, Section 5 gives concluding remarks in this research and points out the directions for the future study of airflows around a moving AGV in cleanrooms.

Modeling
To maintain the cleanliness of the cleanroom, the external airflow entering the cleanroom should be filtered by the high efficiency particular air (HEPA) filter.The clean airflow plays a role of removing microcontaminants and hazardous gas from the cleanroom.A grating zone at the bottom is to let out the airflow and keep a normal pressure.To facilitate the analysis, the moving AGVs (from right to left) are supposed to be still while the flow around the AGV is set to be moving from left to right, as is shown in Figure 1.
The height and width of the cleanroom are set to 2.1 m and 9.2 m, respectively, and a wafer cassette is above the top surface of the AGV, which is consistent with the experiment in [18].The airflow enters the cleanroom from the ceiling with a uniform velocity.Initially, the flow around is stationary, and the velocities of flow at the left side and right side are set to be the same.As the computation starts, the effects of the boundary setting on two sides appear, simulating the AGV moves toward the right side with a constant velocity.To avoid the drawbacks of such a simulation, the pressure conditions on both sides are set to be free.
In order to investigate the distribution of micro contaminant, the tyres are assumed to be the source of the micro containments and the concentration at the bottom of AGV is set as constant.

Governing Equations.
Let Ω be the boundary of a threedimensional polyhedral domain Ω.  1 (Ω) is the Sobolev space of the first order and  2 0 (Ω) is the subspace of  2 (Ω) functions with zero mean value.Under the assumption that the flow field is incompressible, viscous, and laminar, solving the model can be summarized as finding (, ) ∈  1 (Ω) 3 ×  2 (Ω) such that for any  ∈ (0, ), the following set of equations hold: where with the Kronecker delta   and the viscosity  [kg/ms].Let Γ 2 ⊂ Ω; the governing equation of microcontaminants diffusion is to find  ∈  1 (Ω) such that where  is the concentration;  is the diffusion coefficient [m 2 /s], and  is the source term [K/s].An adapted Lagrange-Galerkin method is applied to the nonlinear terms in (1) and ( 4), which works as follows: let Δ be the time increment and   ≡ [/Δ] be the total step number, thus   ≡ Δ, where  is an integer.The trajectory of fluid particle is defined by where (, ) : Ω × (0, ) → R 3 is a smooth function denoting the velocity and () : (0, ) → R 3 is a solution of this ODE, representing the particle's position at .With definition in (5), a particle's position at  −1 can be approximated by  1 ((, ), Δ) : Ω → R 3 with the formula where denotes a first-order approximation; see [24]; therefore, the material derivative term in (1) can be written as follows: Here, the notation ∘ designates the composition of functions.Similarly, (4) can also be replaced in an analogy way.

Stabilized FEM Scheme.
Let the subscript ℎ be the representative length of the triangulation, a triangulation I ℎ of Ω consisting of tetrahedral elements, where The finite element spaces used in this work are as follows: The weak form of the Navier-Stokes equations (1) after finite element discretization can thus be written in oneequation as follows: where (⋅, ⋅) means the inner product and the bilinear forms (⋅, ⋅) and (⋅, ⋅) are given by As is shown in ( 10) and ( 11), 1 element is used for both velocity and pressure, which does not satisfy the so-called LBB compatibility condition [25], a penalty stabilization method [25] is employed, and a stabilization term is added to the left hand side of (11), where   = min{Δ/2, ℎ 2  /24 } is the stabilization parameter.The scheme is as follows:  STEP 2: find It should be noted that the element information term ( −1 ) may cause some difficulty in solving governing equations (1) and ( 4), especially in this work when all the elements information is distributed on different processor elements (PEs) and the searching is done in parallel.Fortunately, in scheme of ( 15) and ( 16), the same  −1 ℎ is used for both equations; therefore, the searching results can be shared by them in one nonstationary step.

Boundary Conditions.
The model is extended to three dimensions in this work and the width in  2 direction is set to 0.5 m, while the AGV's width in  2 direction is 0.4 m.A diagram of the three-dimensional AGV model is given by Figure 2 To compare the numerical results with numerical results and experimental results obtained by Kanayama et al. [18], in which the AGV was moving at a constant speed, a similar boundary setting was used in this research.The first aim is to compare the results with [18] to verify the consistence and to reduce the effects of the end walls; the velocity in  2 direction is set to be zero ( 2 = 0), as is shown in (17).
In order to investigate the end wall effects produced by the three-dimensional model, the restriction in  2 direction is removed and the boundary conditions are given as follows: In this simulation, the tyres are assumed to be the source of the micro containments and the concentration is set to be 0.01 mg/L at the bottom of AGV.

Parallel Implementation and Preconditioning.
In the domain decomposition system used in this work, computation models are firstly divided into several parts before the domain decomposition computation; parts are further decomposed into subdomains, and for each parts, domain decomposition is performed by the current processor element (PE), see Figure 3.
During the computation of PCG loop, PEs work independently almost all the time and only the information of the elements that belong to the current parts is available; however, when computing the residual, the PEs need to work synchronously and the collection of local residuals from all the PEs is performed.
As a hybrid Schwarz type preconditioner, the classical balancing domain decomposition (BDD) preconditioner [26] shows good efficiency in solving the interface problem of domain decomposition system.It can be written as where   and   are the coarse level preconditioner and the local level preconditioner, respectively, and  is the Schur complement of the linear system.In order to obtain   , the inverse of the coarse matrix is required [26].As the dimension of coarse matrix is usually very small, a parallel Cholesky factorization is performed; the coarse matrix is distributed to each PE and thus the computation cost is reduced.During the update of one row of the coarse matrix, each processor is idle until the final updated elements of the same row belonging to the previous processor have been completed.It is clear that in the row-oriented factorization of such sparse coarse matrix,      expensive.A similar phenomena in structural analysis was reported by Ogino et al. [27].
A strategy to neglect the fill-in at the beginning of each distributed portion of the row in some PEs is adopted.The  coarse matrix is modified to keep the sparsity, and thus the preconditioner is changed to

Numerical Results
The criteria of the stationary stage is set by an element based where The convergence judgment of each CG loop is made by the Euclidean norm, quasispectral method and the flow looks more diffusive than others; they all show good agreements within the tendency of flows [28], which convince us the success of the current scheme.
The efficiency of the incomplete balancing domain decomposition preconditioner is also tested by comparing it with other preconditioners and results are shown in Figure 5.
Figure 5(a) shows that comparing with BDD preconditioner and no domain decomposition preconditioners, the IBDD preconditioner is also as fast as the BDD preconditioner, except the fact that the inverse of the coarse matrix is modified, and therefore   is changed to Q .Figure 5(b) demonstrates the time elapsed in solving the coarse problem (indicated by " Coarse Prob.") and the total time elapsed in a CG loop (indicated by " Total").With the increase in interface DOF, the IBBD preconditioner uses less time, and this convinces us its capability in large-scale simulations.

Contamination Control in a Cleaning
Room.The ADVENTURE CAD and ADVENTURE Metis [29] were used to create the geometry and mesh, and the local density of mesh around the AGV was increased to have a better simulation about the flow field around the AGV.In order to have a clear image of the meshing, a coarse mesh is displayed in Figure 6.As is mentioned in Section 3.1, the three-dimensional model is decomposed in several parts by domain decomposition method.By using 9 single-core CPUs (9 PEs in total), a nonoverlapped domain decomposition result of the threedimensional AGV model is demonstrated by Figure 7.
A finer mesh of 920 × 210 × 50 divisions was then created for large-scale computation to improve the accuracy of the simulation and to overcome the numerical difficulty caused by high Reynolds number.A comparison of the velocity vectors around the AGV got by the current solution and Kanayama et al. 's 2D results is shown in Figure 8.
Aiming at getting a stationary stage of the ADV at a constant speed,  1 was fixed in this research.Using the boundary setting in (18), pseudo 2D results in Figure 8(a) give an image about the velocity vectors around the AGV after the coordinate transformation, which agrees with the numerical results of a constant speed AGV reported by Kanayama et al. [18] in Figure 8(b) very well; the upflow from the grating zone does not reach the wafer cassettes.Compared with numerical results reported by Kanayama et al. [18] and Yang et al. [15,16], more details of the flow pattern around the AGV are given by the current large-scale simulation: the eddy behind the wafer cassettes and the AGV is displayed more clearly by Figure 8(a) and the eddy in front of the wafer cassettes, which is caused by the angular pediment of the AGV, is also demonstrated by Figure 8(a).
The experiment in [18] was done in three dimensions, as is shown in Figure 9, where the flow pattern around the AGV when deceleration is displayed, the frog shows threedimensional characters although the nozzles are arranged in one line.However, the numerical simulations were done in two dimensions.This surprising phenomenon motived us to simulate the 3D model and make the comparison of results in different dimensions.With the essential boundary conditions given in (18), true 3D model results are illustrated by Figure 10, which shows the velocity vectors of two planes:  2 = 0.2 (a front plane) and  2 = 0.3 (the central plane).
Compared with the pseudo 2D results in Figure 8(a), it is not difficult to confirm the necessity of using threedimensional modeling in simulations to an AGV moving in a cleanroom.
(1) The eddy behind the AGV and the wafer cassette becomes smaller, which reflects the end-wall effects of Γ side .(2) Similarly, the eddy in the plane  2 = 0.2 is smaller than that of the plane  2 = 0.3, as it is closer to Γ side and more influence is received from the threedimensional boundary settings.
It is worth pointing out that in spite of the fact that the computation domain is not simply-connected, the Lagrange-Galerkin does not encounter any difficulties at these "holes"; both Figures 9 and 10 show good continuity at the corners, convincing us of the solvability of complex problems.
To compare the numerical results of Kanayama et al., the space between the grating zone and the bottom of the AGV is supposed to be the source of pollution first, which generates 1 contaminant per second.The comparison of the 0 0.001 0.002 0.003 0.004 0.006 0.007 0.008 0.009 0.01  The necessity of the vertical external airflow was further proved by numerical experiments with boundaries described in Section 2.3.As is shown in Figure 12, when no external clean air exists, the cassette is polluted by the microcontaminants surrounding the AGV; while with the external vertical clean air, the concentration of micro contaminates on the surface of the cassette is close to zero.
In this work, the concern about the proper setting of the speed of external vertical clean air under boundary conditions described in Section 2.3 had been throughly studied.A series of concentration values of microcontaminants on surface of the cassette were obtained under various settings of  3 , and the relationship between them is displayed in Figure 13.
From Figure 11, it is easy to see that with the increase of  3 , the decreasing speed of contaminant concentration is slowing down; at  3 = 0.7, the concentration is below 1 × 10 −4 mg/L and the trend is turning to be stable; therefore, it can be regarded as a suitable setting of the speed of the external clean air.
The Raynolds number in this model is about 33,000.In this work, by using domain decomposition, an AGV model of over 30 million DOF is solved.Isolines of the vorticity are illustrated in Figure 14.As can be seen, by large-scale computation, more details of the flow induced by the AGV are revealed.The recirculation zones can be more clearly viewed, which will provide more accurate information to remove the microcontaminants.
The model was divided into 920 × 210 × 50 elements (mesh size ℎ = 0.01) and analyzed parallelly by the scheme described in Section 2.2.The total DOF was 39,643,524 and IBDD [23] was employed as the preconditioner to accelerate the CG iteration.A Linux PC cluster with 20 PCs was utilized to carry out the computation, and for each PC the configuration of hardware was as follows: CPU: Intel(R) Core(TM) i7 920@2.67GHz, Memory: 12 [GB].Due to the unconditional stability of the adapted Lagrange-Galerkin method, Δ was set to 0.1 s.The computation took about 20 hours on this small PC cluster.

Conclusions
A moving AGV in a cleanroom is numerically simulated in three dimensions by large-scale computation in this work.The main conclusions can be summarized as follows.
(1) By using large-scale simulation, the results are consistent with the conventional simulations; moreover, more details of the computational models are revealed, which is important for contaminant control in practice.(2) To use three-dimensional modeling is necessary to improve the numerical simulation of moving AGVs in cleanrooms.(3) The proper setting of the speed of external vertical clean air is found, which can be applied to the design and optimization of cleaning room.
(4) The adapted Lagrange-Galerkin method shows good computation accuracy in solving complex problems.
(5) By using incomplete balanced domain decomposition preconditioner, the scheme has the solvability for large-scale problems with up to 30 million of DOF.
The information of recirculation zones is extremely valuable for modern health centers and semiconductor industry to optimize and improve the facilities and to make a clearer cleanroom.In spite of the pervasive applications of AGV systems, the simulation of airflows induced by AGVs still calls for much attention in the future.

Figure 1 :
Figure 1: A moving AGV in a cleanroom.

Figure 6 :
Figure 6: The meshing of the 3D AGV model.
) and in this work,  1 = 1.0 × 10 −4 and  2 = 1.0 × 10 −6 .4.1.Validation and Efficiency.The solver is tested by a 3D lid-driven cavity problem.The current results of Re = 400 (Figure 4(a)) and Re = 1,000 (Figure 4(b)) are also compared with some available 3-dimensional benchmarks.The current results show a better agreement with Tabata et al. and Kato et al., which use finite element method; Ku et al. used