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

Authors: Qing-He Yao; Qing-Yong Zhu
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

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 stiffness matrix 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.


---

## Body

## 1. 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–6], system control [7–10], and path planning [11–14] of AGVs; there are few reports devoted to the simulation of the motion of the airflow in the cleanroom; see [15–21]. By using an upwind scheme, Kanayama et al. analyzed a two-dimensional 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. Section2 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.


## 2. Modeling

To maintain the cleanliness of the cleanroom, the external airflow entering the cleanroom should be filtered by thehigh 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.Figure 1
A moving AGV in a cleanroom.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.


## 3. Formulation and Implementation

### 3.1. Governing Equations

Let∂Ω be the boundary of a three-dimensional polyhedral domainΩ. H1(Ω) is the Sobolev space of the first order and L02(Ω) is the subspace of L2(Ω) 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 (u,p)∈H1(Ω)3×L2(Ω) such that for any t∈(0,T), the following set of equations hold:
(1)∂u∂t+(u·∇)u-1ρ∇·σ(u,p)=1ρfinΩ×(0,T),∇·u=0inΩ×(0,T),u=u^onΓ1×(0,T),∑j=03σijnj=0on∂Ω∖Γ1×(0,T),u=u0inΩ,att=0,
where
(2)Γ1⊂∂Ω,u=(u1,u2,u3) is the velocity [m/s]; p is the pressure [N/m2]; ρ is the density (const.) [kg/m3]; f=(f1,f2,f3) is the body force [N/m3]; u^=(u^1,u^2,u^3) is the boundary velocity [m/s], and σ(u,p) is the stress tensor [N/m2] defined by
(3)σij(u,p)≡-pδij+2μDij(u),Dij(u)≡12(∂ui∂xj+∂uj∂xi),i,j=1,2,3,
with the Kronecker delta δij and the viscosity μ [kg/ms].LetΓ2⊂∂Ω; the governing equation of microcontaminants diffusion is to find C∈H1(Ω) such that
(4)∂C∂t+u·∇C-aΔC=SinΩ×[0,Τ],C=C^onΓ2×[0,Τ],a∂C∂n=0on∂Ω∖Γ2×[0,Τ],C=C0inΩ,att=0,
where C is the concentration; a is the diffusion coefficient [m2/s], and S 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 Δt be the time increment and NT≡[T/Δt] be the total step number, thus tn≡nΔt, where n is an integer. The trajectory of fluid particle is defined by
(5)X′(t)=u(X,t),in(tn-1,tn),X(tn)=x,
where u(X,t):Ω×(0,T)→ℝ3 is a smooth function denoting the velocity and X(t):(0,T)→ℝ3 is a solution of this ODE, representing the particle’s position at t.With definition in (5), a particle’s position at tn-1 can be approximated by X1(u(X,t),Δt):Ω→ℝ3 with the formula
(6)X(tn-1)=X1(u(X,tn-1),Δt)+O(Δt2),
where
(7)X1(u(X,tn-1),Δt)≡x-u(X,tn-1)Δt
denotes a first-order approximation; see [24]; therefore, the material derivative term in (1) can be written as follows:
(8)∂u∂t+(u·∇)u=un-un-1∘X1(un-1,Δt)Δt(x)+O(Δt).Here, the notation∘ designates the composition of functions. Similarly, (4) can also be replaced in an analogy way.


### 3.2. Stabilized FEM Scheme

Let the subscripth be the representative length of the triangulation, a triangulation ℑh of Ω consisting of tetrahedral elements, where
(9)Ω=⋃K∈ℑhK.The finite element spaces used in this work are as follows:(10)Xh≡{vh∈C0(Ω¯)3;vh|K∈P1(K)3,∀K∈ℑh},Mh≡{qh∈C0(Ω¯);qh|K∈P1(K),∀K∈ℑh},Vh(b→)≡{vh∈Xh;vh(P)=b→(P),∀P∈Γ1},Θh(b)≡{θh∈Mh;θh(P)=b(P),∀P∈Γ2},Vh≡Vh(0),Θh≡Θh(0),Qh=Mh.The weak form of the Navier-Stokes equations (1) after finite element discretization can thus be written in one-equation as follows:(11)find{(uhn,phn)}n=1NT∈Vh(u^)×Qh,suchthatforn=1,…,NT,(uhn-uhn-1∘X1(uhn-1,Δt)Δt,vh)+a(uhn,vh)+b(vh,phn)+b(uhn,qh)=(1ρfn,vh),∀(vh,qh)∈Vh×Qh,
where (·,·) means the inner product and the bilinear forms a(·,·) and b(·,·) are given by
(12)a(u,v)≡∫Ω2μρD(u):D(v)dx,(13)b(v,q)≡-∫Ω1ρq∇·vdx,
respectively.As is shown in (10) and (11), P1 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
(14)∑K∈ThτK(1ρ∇phn,-1ρ∇qh)K
is added to the left hand side of (11), where τK=min{Δt/2,ρhk2/24μ} is the stabilization parameter.The scheme is as follows:forn=1,…NT, ∀θh∈Θh and ∀(vh,qh)∈Vh×Qh,STEP 1:findChn∈Θh(C^)by(15)(Chn-Chn-1∘X1(uhn-1,Δt)Δt,θh)+(a∇Chn,∇θh)=(S,θh).STEP 2:find(uhn,phn)∈V(u^)×Qby(16)(uhn-uhn-1∘X1(uhn-1,Δt)Δt,vh)+a0(uhn,vh)+b(vh,phn)+b(uhn,qh)+∑K∈ℑhτK(∇phn,-∇qh)K=(1ρfn,vh).It should be noted that the element information termX(tn-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 uhn-1 is used for both equations; therefore, the searching results can be shared by them in one nonstationary step.


### 3.3. Boundary Conditions

The model is extended to three dimensions in this work and the width inx2 direction is set to 0.5 m, while the AGV’s width in x2 direction is 0.4 m. A diagram of the three-dimensional AGV model is given by Figure 2(17)u1=u2=u3=0onΓAGV,u1=-U1u2=0onΓup+Γdown+Γgrating,u1=-U1u2=0u3=-0.264onΓHEPA,u2=0onΓside.Figure 2
A 3D AGV.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 x2 direction is set to be zero (u2=0), as is shown in (17).In order to investigate the end wall effects produced by the three-dimensional model, the restriction inx2 direction is removed and the boundary conditions are given as follows:
(18)u1=u2=u3=0onΓAGV+Γside,u1=-U1u2=0onΓgrating,u1=-U1u2=0u3=-0.264onΓHEPA,u1=-U1onΓup+Γdown.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.


### 3.4. 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 Figure3.Figure 3
A domain decomposition system.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
(19)MBDD-1=QC+(I-QCS)QL(I-QCS),
where QC  and  QL are the coarse level preconditioner and the local level preconditioner, respectively, and S is the Schur complement of the linear system. In order to obtain QC, 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, fill-in occurs in the position of some zero elements due to the presence of any previous nonzero element in the same row. As a result, the factors of the coarse matrix become less sparse than the original coarse matrix and these make the solution expensive. A similar phenomena in structural analysis was reported by Ogino et al. [27].A strategy to neglect thefill-inat 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
(20)MIBDD-1=QC~+(I-QC~S)QL(I-SQC~).


## 3.1. Governing Equations

Let∂Ω be the boundary of a three-dimensional polyhedral domainΩ. H1(Ω) is the Sobolev space of the first order and L02(Ω) is the subspace of L2(Ω) 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 (u,p)∈H1(Ω)3×L2(Ω) such that for any t∈(0,T), the following set of equations hold:
(1)∂u∂t+(u·∇)u-1ρ∇·σ(u,p)=1ρfinΩ×(0,T),∇·u=0inΩ×(0,T),u=u^onΓ1×(0,T),∑j=03σijnj=0on∂Ω∖Γ1×(0,T),u=u0inΩ,att=0,
where
(2)Γ1⊂∂Ω,u=(u1,u2,u3) is the velocity [m/s]; p is the pressure [N/m2]; ρ is the density (const.) [kg/m3]; f=(f1,f2,f3) is the body force [N/m3]; u^=(u^1,u^2,u^3) is the boundary velocity [m/s], and σ(u,p) is the stress tensor [N/m2] defined by
(3)σij(u,p)≡-pδij+2μDij(u),Dij(u)≡12(∂ui∂xj+∂uj∂xi),i,j=1,2,3,
with the Kronecker delta δij and the viscosity μ [kg/ms].LetΓ2⊂∂Ω; the governing equation of microcontaminants diffusion is to find C∈H1(Ω) such that
(4)∂C∂t+u·∇C-aΔC=SinΩ×[0,Τ],C=C^onΓ2×[0,Τ],a∂C∂n=0on∂Ω∖Γ2×[0,Τ],C=C0inΩ,att=0,
where C is the concentration; a is the diffusion coefficient [m2/s], and S 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 Δt be the time increment and NT≡[T/Δt] be the total step number, thus tn≡nΔt, where n is an integer. The trajectory of fluid particle is defined by
(5)X′(t)=u(X,t),in(tn-1,tn),X(tn)=x,
where u(X,t):Ω×(0,T)→ℝ3 is a smooth function denoting the velocity and X(t):(0,T)→ℝ3 is a solution of this ODE, representing the particle’s position at t.With definition in (5), a particle’s position at tn-1 can be approximated by X1(u(X,t),Δt):Ω→ℝ3 with the formula
(6)X(tn-1)=X1(u(X,tn-1),Δt)+O(Δt2),
where
(7)X1(u(X,tn-1),Δt)≡x-u(X,tn-1)Δt
denotes a first-order approximation; see [24]; therefore, the material derivative term in (1) can be written as follows:
(8)∂u∂t+(u·∇)u=un-un-1∘X1(un-1,Δt)Δt(x)+O(Δt).Here, the notation∘ designates the composition of functions. Similarly, (4) can also be replaced in an analogy way.


## 3.2. Stabilized FEM Scheme

Let the subscripth be the representative length of the triangulation, a triangulation ℑh of Ω consisting of tetrahedral elements, where
(9)Ω=⋃K∈ℑhK.The finite element spaces used in this work are as follows:(10)Xh≡{vh∈C0(Ω¯)3;vh|K∈P1(K)3,∀K∈ℑh},Mh≡{qh∈C0(Ω¯);qh|K∈P1(K),∀K∈ℑh},Vh(b→)≡{vh∈Xh;vh(P)=b→(P),∀P∈Γ1},Θh(b)≡{θh∈Mh;θh(P)=b(P),∀P∈Γ2},Vh≡Vh(0),Θh≡Θh(0),Qh=Mh.The weak form of the Navier-Stokes equations (1) after finite element discretization can thus be written in one-equation as follows:(11)find{(uhn,phn)}n=1NT∈Vh(u^)×Qh,suchthatforn=1,…,NT,(uhn-uhn-1∘X1(uhn-1,Δt)Δt,vh)+a(uhn,vh)+b(vh,phn)+b(uhn,qh)=(1ρfn,vh),∀(vh,qh)∈Vh×Qh,
where (·,·) means the inner product and the bilinear forms a(·,·) and b(·,·) are given by
(12)a(u,v)≡∫Ω2μρD(u):D(v)dx,(13)b(v,q)≡-∫Ω1ρq∇·vdx,
respectively.As is shown in (10) and (11), P1 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
(14)∑K∈ThτK(1ρ∇phn,-1ρ∇qh)K
is added to the left hand side of (11), where τK=min{Δt/2,ρhk2/24μ} is the stabilization parameter.The scheme is as follows:forn=1,…NT, ∀θh∈Θh and ∀(vh,qh)∈Vh×Qh,STEP 1:findChn∈Θh(C^)by(15)(Chn-Chn-1∘X1(uhn-1,Δt)Δt,θh)+(a∇Chn,∇θh)=(S,θh).STEP 2:find(uhn,phn)∈V(u^)×Qby(16)(uhn-uhn-1∘X1(uhn-1,Δt)Δt,vh)+a0(uhn,vh)+b(vh,phn)+b(uhn,qh)+∑K∈ℑhτK(∇phn,-∇qh)K=(1ρfn,vh).It should be noted that the element information termX(tn-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 uhn-1 is used for both equations; therefore, the searching results can be shared by them in one nonstationary step.


## 3.3. Boundary Conditions

The model is extended to three dimensions in this work and the width inx2 direction is set to 0.5 m, while the AGV’s width in x2 direction is 0.4 m. A diagram of the three-dimensional AGV model is given by Figure 2(17)u1=u2=u3=0onΓAGV,u1=-U1u2=0onΓup+Γdown+Γgrating,u1=-U1u2=0u3=-0.264onΓHEPA,u2=0onΓside.Figure 2
A 3D AGV.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 x2 direction is set to be zero (u2=0), as is shown in (17).In order to investigate the end wall effects produced by the three-dimensional model, the restriction inx2 direction is removed and the boundary conditions are given as follows:
(18)u1=u2=u3=0onΓAGV+Γside,u1=-U1u2=0onΓgrating,u1=-U1u2=0u3=-0.264onΓHEPA,u1=-U1onΓup+Γdown.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.


## 3.4. 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 Figure3.Figure 3
A domain decomposition system.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
(19)MBDD-1=QC+(I-QCS)QL(I-QCS),
where QC  and  QL are the coarse level preconditioner and the local level preconditioner, respectively, and S is the Schur complement of the linear system. In order to obtain QC, 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, fill-in occurs in the position of some zero elements due to the presence of any previous nonzero element in the same row. As a result, the factors of the coarse matrix become less sparse than the original coarse matrix and these make the solution expensive. A similar phenomena in structural analysis was reported by Ogino et al. [27].A strategy to neglect thefill-inat 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
(20)MIBDD-1=QC~+(I-QC~S)QL(I-SQC~).


## 4. Numerical Results

The criteria of the stationary stage is set by an element basedH1(Ω3)×L2(Ω)×H1(Ω) norm,(21)diff=∥(un,pn,Cn)-(un-1,pn-1,Cn-1)∥H1(Ω3)×L2(Ω)×H1(Ω)∥(un-1,pn-1,Cn-1)∥H1(Ω3)×L2(Ω)×H1(Ω)≤ε1,
where
(22)∥(u,p,C)∥H1(Ω3)×L2(Ω)×H1(Ω)≡1Re∥u∥H1(Ω)3+∥p∥L2(Ω)+∥C∥H1(Ω).
The convergence judgment of each CG loop is made by the Euclidean norm,
(23)∥g(i)∥2∥g(0)∥2≤ε2,
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 (Figure4(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 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.Comparison ofu1 and u3 at the plane of x2=0.5. ((a): 50×50×50, Re = 400; (b): 80×80×80, Re = 1,000).
(a)
(b)The efficiency of the incomplete balancing domain decomposition preconditioner is also tested by comparing it with other preconditioners and results are shown in Figure5.IBDD versus BDD.
(a)
(b)Figure5(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 QC is changed to QC~.   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.


### 4.2. 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.Figure 6
The meshing of the 3D AGV model.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 three-dimensional AGV model is demonstrated by Figure7.Figure 7
The domain decomposed 3D AGV model.A finer mesh of920×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.Velocity vectors around the AGV (a) current 3D results; (b) Kanayama et al.’s 2D results.
(a)
(b)Aiming at getting a stationary stage of the ADV at a constant speed,U1 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 three-dimensional 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: x2=0.2 (a front plane) and x2=0.3 (the central plane).Figure 9
Kanayama et al.’s experimental visualization.Velocity vectors around 3D AGV model.
(a)
(b)Compared with the pseudo 2D results in Figure8(a), it is not difficult to confirm the necessity of using three-dimensional 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 planex2=0.2 is smaller than that of the plane x2=0.3, as it is closer to Γside and more influence is received from the three-dimensional 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 Figures9 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 amount concentration on the upper face of AGV is illustrated in Figure11. By changing the velocity of external air entering the cleaning room, the distributions are compared and both the current results and Kanayama et al.’s agree quite well on the effect of the external air. Due to the improvement in accuracy, the current results look to be more smooth than Kanayama et al.’s.Figure 11
Contaminants distribution.The necessity of the vertical external airflow was further proved by numerical experiments with boundaries described in Section2.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.Isolines of contaminants concentration. (a) Without vertical external air; (b) with vertical external air.
(a)
(b)In this work, the concern about the proper setting of the speed of external vertical clean air under boundary conditions described in Section2.3 had been throughly studied. A series of concentration values of microcontaminants on surface of the cassette were obtained under various settings of u3, and the relationship between them is displayed in Figure 13.Figure 13
The effects ofu3.From Figure11, it is easy to see that with the increase of u3, the decreasing speed of contaminant concentration is slowing down; at u3=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 Figure14. 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.Figure 14
Isolines of vorticity.The model was divided into920×210×50 elements (mesh size h=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.67 GHz,Memory: 12 [GB].Due to the unconditional stability of the adapted Lagrange-Galerkin method,Δt was set to 0.1 s. The computation took about 20 hours on this small PC cluster.


## 4.1. Validation and Efficiency

The solver is tested by a 3D lid-driven cavity problem. The current results of Re = 400 (Figure4(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 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.Comparison ofu1 and u3 at the plane of x2=0.5. ((a): 50×50×50, Re = 400; (b): 80×80×80, Re = 1,000).
(a)
(b)The efficiency of the incomplete balancing domain decomposition preconditioner is also tested by comparing it with other preconditioners and results are shown in Figure5.IBDD versus BDD.
(a)
(b)Figure5(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 QC is changed to QC~.   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.


## 4.2. 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.Figure 6
The meshing of the 3D AGV model.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 three-dimensional AGV model is demonstrated by Figure7.Figure 7
The domain decomposed 3D AGV model.A finer mesh of920×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.Velocity vectors around the AGV (a) current 3D results; (b) Kanayama et al.’s 2D results.
(a)
(b)Aiming at getting a stationary stage of the ADV at a constant speed,U1 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 three-dimensional 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: x2=0.2 (a front plane) and x2=0.3 (the central plane).Figure 9
Kanayama et al.’s experimental visualization.Velocity vectors around 3D AGV model.
(a)
(b)Compared with the pseudo 2D results in Figure8(a), it is not difficult to confirm the necessity of using three-dimensional 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 planex2=0.2 is smaller than that of the plane x2=0.3, as it is closer to Γside and more influence is received from the three-dimensional 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 Figures9 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 amount concentration on the upper face of AGV is illustrated in Figure11. By changing the velocity of external air entering the cleaning room, the distributions are compared and both the current results and Kanayama et al.’s agree quite well on the effect of the external air. Due to the improvement in accuracy, the current results look to be more smooth than Kanayama et al.’s.Figure 11
Contaminants distribution.The necessity of the vertical external airflow was further proved by numerical experiments with boundaries described in Section2.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.Isolines of contaminants concentration. (a) Without vertical external air; (b) with vertical external air.
(a)
(b)In this work, the concern about the proper setting of the speed of external vertical clean air under boundary conditions described in Section2.3 had been throughly studied. A series of concentration values of microcontaminants on surface of the cassette were obtained under various settings of u3, and the relationship between them is displayed in Figure 13.Figure 13
The effects ofu3.From Figure11, it is easy to see that with the increase of u3, the decreasing speed of contaminant concentration is slowing down; at u3=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 Figure14. 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.Figure 14
Isolines of vorticity.The model was divided into920×210×50 elements (mesh size h=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.67 GHz,Memory: 12 [GB].Due to the unconditional stability of the adapted Lagrange-Galerkin method,Δt was set to 0.1 s. The computation took about 20 hours on this small PC cluster.


## 5.  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.



---

*Source: 570237-2013-02-26.xml*