A numerical investigation ofe!ects ofa moving operator on air#ow
patterns in a cleanroom
Suh-Jenq Yang
a, Wu-Shung Fu
b;∗aDepartment of Industrial Engineering and Management, Nan Kai College, Tsaotun, Nantou, 54210, Taiwan bDepartment of Mechanical Engineering, National Chiao Tung University, Hsinchu 30056, Taiwan, ROC
Received 25 May 2000; received in revised form 21 June 2000; accepted 31 July 2001
Abstract
The variations ofthe air#ow induced by a moving operator in a cleanroom installed with a curtain were studied numerically. This situation is cataloged to a class ofthe moving boundary problems. An arbitrary Lagrangian–Eulerian kinematic description method is utilized to describe the #ow 5eld and a penalty 5nite element formulation with moving meshes is adopted to solve this problem. The e!ects ofthe moving operator and curtain on the air#ow patterns under di!erent distances from the workbench to the curtain with the moving speed ofthe operator equal to 2:0 and Reynolds number Re = 500 are taken into account. The results show that recirculation zones are formed around the operator and workbench due to the movement of the operator. The recirculation zones are not favorable to the cleanroom because they may induce a local turbulent #ow and entrain and trap contaminants. These phenomena are remarkably di!erent from those of the moving operator assumed as stationary in the cleanroom. Based on the length of the curtain, it is useful for protecting the operator from the hazardous gases. c 2002 Elsevier Science Ltd. All rights reserved.
Keywords: Cleanroom; Curtain; Moving operator; ALE; Recirculation zones
1. Introduction
Cleanroom recently became an indispensable envi-ronment in the pharmaceutical, medical technology, and semiconductor industries for producing high quality and precision products. Since the external air entering the clean-room must be 5ltered by the HEPA or ULPA 5lter banks, the operators and equipment become the major sources of the contaminants in the cleanroom [1]. Most particles gen-erated by the above sources are continuously swept away by the air#ow from the ceiling of the cleanroom. Some residual particles, which may be suspended within recircu-lation zones or deposited on the products and equipment by gravitational settling, di!usion, collision, and electrostatic attraction, etc., are extremely di?cult to be removed by the air#ow. How to remove these residual particles e!ectively becomes an important issue in the semiconductor industry.
Concerning the motions ofthe air#ow and particles in the cleanroom, several studies investigated this subject. Er-mak and Buckholz [2] adopted a Monte Carlo method to simulate the e!ects ofthe air#ow on the particles, and the
∗Corresponding author. Tel.: +886-3-5712121; fax: 886-3-5720634.
E-mail address: wsfu@cc.nctu.edu.tw (W.-S. Fu).
results showed that the characteristics ofthe particle trans-port were dominated by the air#ow. Kuehn [3], Yamamoto et al. [4], Busnaina et al. [5], and Lemaire and Luscure [6] utilized the numerical modeling to study the air#ow and par-ticle transport in the cleanroom. Liu and Ahn [7] used the analogy between mass transfer and heat transfer to deter-mine particle deposition rates by di!usion. Settles and Via [8] adopted Schloeren observation to observe the #ow paths ofthe particles in the cleanroom. Furthermore, Marvell [9] summarized the factors of e!ect on the air#ow and contam-inant transport in a minienvironment system. Tannous [10] utilized the computational #uid dynamics method to inves-tigate the #ow 5eld ofa minienvironment in the cleanroom. However, for facilitating the analysis, most previous stud-ies regarded a moving operator as a stationary object, which resulted in the phenomena ofthe air#ow in the cleanroom being rather di!erent from the realistic situation. Besides, a curtain is usually employed in the cleanroom to divide the inlet ofthe air#ow into two sections and protect the operator from hazardous gases while processing. One sec-tion uses well-5ltered air at a higher #ow rate to sweep away the particles suspended near the working area and the other uses air 5ltered at a lower #ow rate for saving en-ergy. To the knowledge ofthe authors, the e!ects ofboth
0360-1323/02/$ - see front matter c 2002 Elsevier Science Ltd. All rights reserved.
Notation
h0 dimensional height ofthe operator [m]
H0 dimensionless height ofthe operator
h1 dimensional height ofthe cleanroom [m]
H1 dimensionless height ofthe cleanroom
h2 dimensional distance from the outlet of the
cleanroom to the operator [m]
H2 dimensionless distance from the outlet of the
cleanroom to the operator
h3 dimensional distance from the workbench to the
curtain [m]
H3 dimensionless distance from the workbench to
the curtain
h4 dimensional height ofthe workbench [m]
H4 dimensionless height ofthe workbench
p dimensional pressure [N=m2]
p∞ reference pressure [N=m2] P dimensionless pressure
Re Reynolds number
t dimensional time [s]
u; v dimensional velocities ofthe air#ow in x- and y-direction [m=s]
U; V dimensionless velocities ofthe air#ow in X - and Y -direction
ub dimensional moving velocity ofthe operator
[m=s]
Ub dimensionless moving velocity ofthe operator
ˆu dimensional mesh velocity in x-direction [m=s] ˆ
U dimensionless mesh velocity in X -direction v0 dimensional air#ow inlet velocity at section CD
[m=s]
V0 dimensionless air#ow inlet velocity at section
CD
v1 dimensional air#ow inlet velocity at section AB
[m=s]
V1 dimensionless air#ow inlet velocity at section
AB
w dimensional width ofthe cleanroom [m] W dimensionless width ofthe cleanroom w0 dimensional width ofthe operator [m]
W0 dimensionless width ofthe operator
w1 dimensional width ofthe air#ow inlet at section
AB [m]
W1 dimensionless width ofthe air#ow inlet at
sec-tion AB
w2 dimensional width ofthe curtain [m]
W2 dimensionless width ofthe curtain
w3 dimensional width ofthe air#ow inlet at section
CD [m]
W3 dimensionless width ofthe air#ow inlet at
sec-tion CD
w4 dimensional width ofthe workbench [m]
W4 dimensionless width ofthe workbench
w5 dimensional distance from the workbench to the
operator [m]
W5 dimensionless distance from the workbench to
the operator
x; y dimensional Cartesian coordinates [m] X; Y dimensionless Cartesian coordinates Greek symbols
penalty parameter
kinematic viscosity [m2=s]
density [kg=m3]
dimensionless time
dimensionless stream function Superscripts (e) element m iteration number T transpose matrix Other [ ] matrix { } column vector row vector
the moving operator and curtain on the air#ow are hardly investigated, and minimum literature is available on this subject.
Consequently, the aim ofthis paper is to investigate the variations ofthe air#ow patterns induced by the move-ment ofthe operator in the cleanroom numerically. In ad-dition, three di!erent lengths ofthe curtain are taken into consideration to examine the e!ects ofthe length ofthe curtain on the air#ow patterns. Due to the movement of the operator, this problem is classi5ed into a class ofthe
moving boundary problems, and is hardly analyzed by either the Lagrangian or Eulerian kinematic description method solely. An arbitrary Lagrangian–Eulerian (ALE) kinematic description method [11–13], which combines the character-istics ofthe Lagrangian and Eulerian methods, is an appro-priate method to describe this problem. In the ALE method, the computational meshes may move with the #uid (La-grangian), be held 5xed (Eulerian), or be moved in any other prescribed way. The details ofthe ALE method are delin-eated by Hirt et al. [11], Hughus et al. [12], and Ramaswamy
[13]. A Galerkin 5nite element method with moving meshes and a backward di!erence scheme, dealing with the time terms, are used to solve the governing equations. The results show that the length ofthe curtain and the movement ofthe operator signi5cantly in#uence the air#ow patterns in the cleanroom. These phenomena are quite di!erent from those regarded the moving operator as a stationary one.
2. Physical model
A two-dimensional vertical laminar #ow cleanroom as sketched in Fig. 1 is used. The width and height ofthe clean-room are w(= w1+ w2+ w3) and h1, respectively. A
rectan-gular block with height h0 and width w0is used to simulate
an operator in the cleanroom. The distance from the outlet ofthe cleanroom to the bottom surface ofthe operator is h2.
A workbench with height h4and width w4(= w1+w2) is set
on the left side of the cleanroom. A curtain with width w2is
mounted on the inlet ceiling and divides the inlet air#ow into two sections, AB and CD. The distance from the workbench to the curtain is h3. Two di!erent inlet air velocities are v1
and v0 #owing separately at the inlet sections AB and CD.
Initially (t = 0), the operator is stationary and the air#ow is #owing steadily, and the distance from the workbench to the operator is w5. As the time t ¿ 0, the operator moves to the
workbench with a constant velocity ub and remains beside
the workbench. Finally, the operator leaves the workbench
Fig. 1. Physical model.
and moves back to the original place. The interaction tween the inlet air#ow and moving operator a!ects the be-havior ofthe air#ow patterns and contaminant di!usion in the cleanroom. This problem becomes time-dependent and can be cataloged to a class ofmoving boundary problems. As a result, the ALE method is properly utilized to analyze this problem.
In order to facilitate the analysis, the following assump-tions are made.
(1) The #uid is air and the #ow 5eld is two-dimensional, incompressible and laminar.
(2) The #uid properties are constant and the e!ect ofthe gravity is neglected.
(3) The no-slip condition is held at the interface between the #uid and operator.
Based upon the characteristic scales of w0; v0 and v20, the
dimensionless variables are de5ned as follows: X =wx 0; Y = y w0; U = u v0; V = v v0; ˆ U = ˆu v0; Ub= ub v0; P =p − p∞ v2 0 ; =tv0 w0; Re = v0w0 ; (1)
where ˆu is the mesh velocity and ub is the moving velocity
ofthe operator.
According to the above assumptions and dimensionless variables, the dimensionless ALE governing equations are expressed as the following equations:
Continuity equation @U @X + @V @Y = 0; (2) Momentum equations @U @ + (U − ˆU) @U @X + V @U @Y = −@X@P +Re1 @2U @X2 + @2U @Y2 ; (3) @V @ + (U − ˆU) @V @X + V @V @Y = −@Y@P +Re1 @2V @X2 + @2V @Y2 : (4)
As the time ¿ 0, the boundary conditions are as follows: On the wall surfaces DE, FG, GH, and AH
U = V = 0: (5)
On the curtain surfaces BR, CS, and RS,
On the air#ow inlet section AB
U = 0; V = − 1:25: (7)
On the air#ow inlet section CD
U = 0; V = − 1:0: (8)
On the air#ow outlet section EF
@U=@Y = @V=@Y = 0: (9)
On the interface of the operator and the air#ow
U = Ub; V = 0: (10)
3. Numerical method
A Galerkin 5nite element method with moving meshes and a backward scheme to deal with the time terms are adopted to solve the governing equations (2)–(4). A penalty function [14] and the Newton–Raphson iteration algorithm are utilized to simplify the pressure and nonlinear terms in the momentum equations, respectively. The velocity terms are approximated by quadrilateral and nine-node quadratic isoparametric elements. The discretization processes ofthe governing equations are similar to those used by Fu et al. [15] and Huebner et al. [16]. Then, the momentum equations (3) and (4) can be expressed as follows:
ne
l
([A](e)+ [K](e)+ [L](e)){q}(e) +P= ne l {f}(e); (11) where ({q}(e)+P)T= U 1; U2; : : : ; U9; V1; V2; : : : ; V9m+1+P; (12)
[A](e)includes the (m)th iteration values of U and V at time
+ P,
[K](e)includes the shape function, ˆU and time di!erential
terms,
[L](e)includes the penalty function terms,
{f}(e)includes the known values of U and V at time
and (m)th iteration values of U and V at time + P. In Eq. (11), the terms with the penalty parameter, , are integrated by 2×2 Gaussian quadrature, and the other terms are integrated by 3 × 3 Gaussian quadrature. The value of penalty parameter used in this study is 106 and the frontal
method solver [17,18] is utilized to solve Eq. (11). Concerning the mesh velocity ˆU, it is linearly distributed and inversely proportional to the distance between the nodes ofthe computational meshes and the operator in this study. The mesh velocity near the operator is faster than that near the boundaries ofthe computational domain. In addi-tion, the boundary layer thickness on the operator surface is extremely thin and can be approximately estimated by Re−1=2[19]. To avoid the computational nodes in the vicinity ofthe operator slipping away from the boundary layer, the
mesh velocities adjacent to the operator are expediently as-signed equal to the velocity ofthe operator.
A briefoutline ofthe solution procedures are described as follows:
(1) Determine the optimal mesh distribution and number ofthe elements and nodes.
(2) Solve the values of U and V at the steady state and regard them as the initial values.
(3) Determine the time increment P and the mesh veloc-ities ˆU at every node.
(4) Update the coordinates ofthe nodes and examine the determinant ofthe Jacobian transformation matrix to ensure the one-to-one mapping to be satis5ed during the Gaussian quadrature numerical integration, otherwise, execute the mesh reconstruction.
(5) Solve Eq. (11), until the following criteria for conver-gence are satis5ed:
$m+1$m+1− $m +P ¡ 10−3; where $ = U; V (13) (6) Continue the next time step calculation until the
as-signed position ofthe operator is reached.
4. Results and discussion
The length ofthe curtain plays an important role for pro-tecting the wafer on the workbench, than three di!erent di-mensionless distances from the workbench to the curtain, H3= 6:0, 3.0, and 2.0, which correspond to cases 1, 2 and 3,
respectively, is listed in Table 1. At time = 0:0, the distance from the operator to the workbench is W5= 4:5. For
satisfy-ing the velocity boundary conditions at the inlet and outlet, the dimensionless lengths of H1(= h1=w0) and H2(= h2=w0)
are determined by numerical tests.
Generally, the walking speed ofthe operator is twice the inlet speed ofthe air#ow in the cleanroom. The phenomena ofthe moving speed ofthe operator being equal to 2:0 (the moving velocity ofthe operator, Ub, may be 2:0 or −2:0) and
Re = 500 are analyzed in detail. For obtaining an optimal computational mesh, three di!erent nonuniform distribution elements 3612, 4048, and 4804 (corresponding to 14816, 16572, and 19640 nodes, respectively) are tested for case 2 at the steady state. The results ofthe velocities U and
Table 1
The dimensionless geometric lengths ofthe cleanroom for three di!erent cases
W0 W1 W2 W3 W4 H0 H1 H2 H3 H4 case 1 1.0 2.3 0.2 12.5 2.5 6.0 30.0 18.0 6.0 6.0 case 2 1.0 2.3 0.2 12.5 2.5 6.0 34.0 18.0 3.0 6.0 case 3 1.0 2.3 0.2 12.5 2.5 6.0 30.0 18.0 2.0 6.0
Fig. 2. Comparison ofthe distributions ofthe velocities U and V along the line MN for case 2 at steady state under di!erent meshes. (a) X -U, and (b) X -V .
V distributions along the line MN as indicated in Fig. 1 are shown in Fig. 2. Based upon the results, the computational mesh with 4804 elements is adopted for case 2. Similarly, the computational meshes with 4644 and 4560 elements (corresponding to 18964 and 18664 nodes, respectively) are adopted for cases 1 and 3, respectively.
As for the selection of the time step P, three di!erent time steps 0:05, 0:01 and 0:005 at Ub= − 2:0 are tested for
case 2. The distributions ofthe velocities U and V along the line MN at time = 2:0 are shown in Fig. 3. The result ofthe velocities U and V for di!erent time steps is quite consistent. The time step P = 0:01 is chosen for all cases.
The dimensionless stream function is de5ned as
U =@@Y and V = − @@X: (14)
For illustrating the #ow 5eld clearly, the phenomena ofthe streamlines around the work region ofthe cleanroom are
Fig. 3. Comparison ofthe distribution ofthe velocities U and V along the line MN at time = 2:0 for case 2 at Ub= − 2:0 under di!erent time steps (a) X -U, and (b) X -V .
presented exclusively. However, it should be noted that the computational domain included a much larger region than what is displayed in the subsequent 5gures.
The transient developments ofthe streamline distributions for case 1 (H3= 6:0) are shown in Fig. 4. In this case, the
height ofthe top surface ofthe operator is lower than the bottom surface of the curtain. At time = 0:0, as shown in Fig. 4(a), the operator is stationary and the air#ow is #owing steadily. Recirculation zones are found near the top and lateral surfaces of the workbench. As time ¿ 0:0, the operator starts to move towards the workbench with a con-stant velocity Ub=−2:0, and the variations ofthe #ow 5eld
attain a transient state. As shown in Fig. 4(b), the space be-tween the workbench and operator is contracted gradually, which results in the inlet air#ow beginning to #ow over the rear region ofthe operator. Since the operator moves to-ward the workbench, the operator pushes the air#ow before the operator, and the direction ofthis air#ow is forced to be
Fig. 4. The transient developments ofthe streamline distributions around the work region ofthe cleanroom for case 1 (a) = 0:0, Ub= 0:0, (b) = 0:2, Ub= − 2:0, (c) = 1:0, Ub= − 2:0, (d) = 2:0, Ub= − 2:0, (e) = 2:2, Ub= 0:0, (f) = 4:5, Ub= 0:0, (g) = 7:0, Ub= 0:0, (h) = 7:2, Ub= 2:0, (i) = 8:0, Ub= 2:0, and (j) = 9:0, Ub= 2:0.
changed and to #ow toward the workbench remarkably. In the meantime, the air#ow in the vicinity ofthe rear region ofthe operator and the top and bottom surfaces ofthe op-erator simultaneously replenishes the vacant space induced by the movement ofthe operator. As a result, new recircu-lation zones are formed around the operator. As the time
increases, the space between the workbench and operator becomes narrower, then most inlet air#ow #ows through the rear region ofthe operator, as shown in Figs. 4(c)–(d). In this duration, the air#ow induced by the movement ofthe operator is forced to #ow over the top surface of the opera-tor, which causes the recirculation zones in the rear region ofthe operator to enlarge gradually. From a viewpoint of #uid mechanics, these recirculation zones are not favorable to the cleanroom because they may induce a local turbu-lent #ow near the operator and workbench and entrain and trap particles. Particles may escape from these recirculation zones due to the e!ects ofthe inertia force, gravitational set-tling, turbulent di!usion, Brownian di!usion, electrostatic force, or other forces and deposit on the workbench to pol-lute the products. These phenomena are very complex and can hardly be predicted by both the experimental and nu-merical methods.
For the duration oftime from 2.0 to 7:0, the operator stays beside the workbench (Ub= 0:0), as indicated in Figs.
4(e)–(g). Since the space between the workbench and op-erator is narrow, the air#ow passing through this space is slight and most ofthe air#ow from the inlet #ow through the rear region ofthe operator. As a result, the recirculation zones around the operator migrate to the downstream and shrink gradually by the air#ow from the inlet. But large re-circulation zones are observed around the top surface of the workbench. This #ow may cause the particles to deposit on the product.
As time ¿ 7:0, as shown in Figs. 4(h)–(j), the oper-ator begins to leave the workbench with a constant veloc-ity Ub= 2:0. Because ofthe operator moving toward the
right, the air#ow from the inlet is a!ected by the movement ofthe operator and also #ows toward the right, which causes the recirculation zones near the workbench and operator to be destroyed gradually. Furthermore, the space mentioned above becomes broad gradually, and the inlet air#ow easily passes through this space and replenishes the vacant space induced by the movement ofthe operator. Thus, new re-circulation zones appear near the top and bottom surfaces ofthe operator and extend to the lower region beside the workbench.
Fig. 5 shows the transient developments ofthe stream-line distributions for case 2 (H3= 3:0) in which the height
ofthe top surface ofthe operator is equal to the bottom sur-face ofthe curtain. At the beginning ofthe transient state, the operator moves toward the workbench with a constant velocity Ub= − 2:0 and these phenomena are similar to
case 1. As the time increases, the space between the work-bench and operator is contracted gradually, and some air-#ows from the inlet section AB circumvent the curtain and #ow over the rear region ofthe operator (Fig. 5(c)). As a result, the curtain protects the operator from the hazardous gases. Later, the air#ow #owing around the curtain interacts with the air#ow from the section CD, and large recirculation zones appear around the operator while small new recircu-lation zones are observed around the curtain, as shown in
Fig. 5. The transient developments ofthe streamline distributions around the work region ofthe cleanroom for case 2 (a) = 0:0, Ub= 0:0, (b) = 0:2, Ub= − 2:0, (c) = 1:0, Ub= − 2:0, (d) = 2:0, Ub= − 2:0, (e) = 2:2, Ub= 0:0, (f) = 4:5, Ub= 0:0, (g) = 7:0, Ub= 0:0, (h) = 7:2, Ub= 2:0, (i) = 8:0, Ub= 2:0, and (j) = 9:0, Ub= 2:0.
Fig. 5(d). Besides, due to the existence ofthe curtain, small recirculation zones around the top surface of the workbench cannot be destroyed by the air#ow, which is di!erent from that ofcase 1 as shown in Fig. 4(d).
During time from 2.0 to 7:0, as shown in Figs. 5(e)–(g), the operator remains beside the workbench (Ub= 0:0). Some
air#ows from the inlet section AB #ow through the space between the curtain and operator, and then #ow over the rear region ofthe operator. Thus, the recirculation zones around the curtain enlarge gradually and extend farther, as shown in Figs. 5(e)–(f). As the time increases, these recirculation zones are destroyed by the air#ow circumventing the curtain and the air#ow from the inlet section CD, as indicated in Fig. 5(g). Furthermore, due to the existence ofthe curtain, the recirculation zones near the top surface of the workbench are shrunk gradually, which is bene5cial to the contamination control.
As time ¿ 7:0, as shown in Figs. 5(h)–(j), the operator leaves the workbench with a constant velocity Ub= 2:0. In
this situation, the variations ofthe #ow 5elds are similar to case 1 as shown in Figs. 4(h)–(j).
Fig. 6 shows the transient developments ofthe streamline distributions for case 3 (H3= 2:0) in which the height ofthe
top surface of the operator is higher than the bottom surface ofthe curtain. The variations ofthe air#ow patterns for this case are more slightly drastic than that ofcase 2. As the operator moves toward the workbench, the space between the curtain and operator is contracted gradually. As Figs. 6(b)–(d) show, due to the obstruction ofthe curtain, the air#ow from the inlet section AB circumventing the curtain is di?cult and mass #ow rate ofthe air#ow through the space between the curtain and operator is smaller than that ofcase 2. Consequently, the recirculation zones around the top surface of the workbench is shrunk. Besides, the operator seems to be an obstruction in the way ofthe air#ow passing through the space between the curtain and workbench, then the variations ofthe recirculation zones behind the operator become more apparent than the former ones.
5. Conclusions
The e!ects ofa moving operator and curtain on the varia-tions ofthe air#ow patterns in a vertical laminar cleanroom are investigated numerically. The results can be summarized as follows:
1. The moving operator a!ects the variations ofthe air#ow patterns in the cleanroom very much. Recirculation zones are observed around the operator and workbench as the operator approaches the workbench. These phenomena are remarkably di!erent from those of the moving oper-ator regarded as a stationary object in the cleanroom of the previous studies.
2. The curtain can usually protect the operator from the hazardous gases as the distance between the workbench and curtain is small.
3. The curtain may con5ne the #owing ofthe air#ow from the inlet and force the recirculation zones around the top surface of the workbench shrunk gradually when the operator remained beside the workbench.
Fig. 6. The transient developments ofthe streamline distributions around the work region ofthe cleanroom for case 3 (a) = 0:0, Ub= 0:0, (b) = 0:2, Ub= − 2:0, (c) = 1:0, Ub= − 2:0, (d) = 2:0, Ub= − 2:0, (e) = 2:2, Ub= 0:0, (f) = 4:5, Ub= 0:0, (g) = 7:0, Ub= 0:0, (h) = 7:2, Ub= 2:0, (i) = 8:0, Ub= 2:0, and (j) = 9:0, Ub= 2:0.
Acknowledgements
The support ofthis work by the National Science Council ofTaiwan, R.O.C., under contract NSC90-2626-E-252-003 is gratefully acknowledged.
References
[1] Gilmore JS. Ultra5ne particles. In: Donovan RP, editor. Particle control for semiconductor manufacturing. New York: Marcel Dekker Inc., 1990. p. 79–103.
[2] Ermak DL, Buckholz H. Numerical integration ofthe Langerin equation: Monte Carlo simulation. Journal ofComputational Physics 1980;35:169–82.
[3] Kuehn TH. Computer simulation ofair#ow and particle transport in cleanrooms. Journal ofEnvironmental Science 1988;31(5):21–7. [4] Yamamoto T, Donovan RP, Ensor DS. Model study for
optimization ofcleanroom air#ow. Journal ofEnvironmental Science 1988;31(6):24–9.
[5] Busnaina AA, Abuzeid S, SharifMAR. Numerical modeling of #uid #ow and particle transport in clean room. Proceedings of the Ninth International Committee Contamination Control Societies Conference, Los Angeles, 1988. p. 600–7.
[6] Lemaire T, Luscuere P. Investigating computer modeling of cleanroom air#ow patterns. Microcontamination 1991;9(8):19–26. [7] Liu BYH, Anh K. Particle deposition on semiconductor wafers.
Aerosol Science and Technology 1987;6(3):215–24.
[8] Settles GS, Via GG. A portable Schlieren system for clean-room air-#ow analysis. Journal ofEnvironmental Science 1988;30(5):17–21. [9] Marvell G. Minienvironment air #ow dynamics. Solid State
Technology 1993;36(8):47–8.
[10] Tannous AG. Air #ow simulation in a minienvironment. Solid State Technology 1996;39(7):201–9.
[11] Hirt CW, Amsden AA, Cooks HK. An arbitrary Lagrangian–Eulerian computing method for all #ow speeds. Journal of Computational Physics 1974;14:227–53.
[12] Hughes TJR, Liu WK, Zimmermann TK. Lagrangian–Eulerian 5nite element formulation for incompressible viscous #ows. Computer Methods in Applied Mechanics and Engineering 1981;29:329–49. [13] Ramaswamy B. Numerical simulation ofunsteady viscous free
surface #ow. Journal of Computational Physics 1990;90:396–430. [14] Reddy JN, Gartling DK. The 5nite element method in heat transfer
and #uid dynamics. Ann. Arbor: CRC Press. Inc., 1994. p. 145–51. [15] Fu WS, Kau TM, Shieh WJ. Transient laminar natural convection in an enclosure from steady #ow state to stationary state. Numerical Heat Transfer A 1990;18:189–211.
[16] Hueber KH, Thornton EA, Byrom TG. The 5nite element method for engineers. 3rd ed. New York: Wiley, 1995. p. 406–15. [17] Irons BM. A frontal solution program for 5nite element analysis.
International Journal for Numerical Methods in Engineering 1970;2:5–32.
[18] Taylor C, Hughes TG. Finite element programming ofthe Navier– Stokes Equations. UK: Pineridge Press Ltd., 1981 [Chapter 6]. [19] Schlichting H. Boundary layer theory, 7th ed. New York:
