## Iterative Methods for Nonlinear Equations or Systems and Their Applications

View this Special IssueResearch Article | Open Access

Qing-He Yao, Qing-Yong Zhu, "Investigation of the Contamination Control in a Cleaning Room with a Moving AGV by 3D Large-Scale Simulation", *Journal of Applied Mathematics*, vol. 2013, Article ID 570237, 10 pages, 2013. https://doi.org/10.1155/2013/570237

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

**Academic Editor:**Alicia Cordero

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

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

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

#### 3. Formulation and Implementation

##### 3.1. Governing Equations

Let be the boundary of a three-dimensional polyhedral domain. is the Sobolev space of the first order and is the subspace of 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 such that for any , the following set of equations hold:
where
is the velocity [m/s]; is the pressure [N/m^{2}]; is the density (const.) [kg/m^{3}]; is the body force [N/m^{3}]; is the boundary velocity [m/s], and is the stress tensor [N/m^{2}] defined by
with the Kronecker delta and the viscosity [kg/ms].

Let ; the governing equation of microcontaminants diffusion is to find 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 is a smooth function denoting the velocity and is a solution of this ODE, representing the particle’s position at .

With definition in (5), a particle’s position at can be approximated by 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.

##### 3.2. Stabilized FEM Scheme

Let the subscript be the representative length of the triangulation, a triangulation 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 one-equation as follows: where means the inner product and the bilinear forms and are given by respectively.

As is shown in (10) and (11), 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 is the stabilization parameter.

The scheme is as follows: for , and STEP 1: STEP 2:

It should be noted that the element information term 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 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 in direction is set to 0.5 m, while the AGV’s width in 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 direction is set to be zero (), as is shown in (17).

In order to investigate the end wall effects produced by the three-dimensional model, the restriction in 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.

##### 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 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, *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 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

#### 4. Numerical Results

The criteria of the stationary stage is set by an element based norm, where The convergence judgment of each CG loop is made by the Euclidean norm, and in this work, and .

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

**(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 Figure 5.

**(a)**

**(b)**

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

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 Figure 7.

A finer mesh of 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.

**(a)**

**(b)**

Aiming at getting a stationary stage of the ADV at a constant speed, 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: (a front plane) and (the central plane).

**(a)**

**(b)**

Compared with the pseudo 2D results in Figure 8(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 . (2)Similarly, the eddy in the plane is smaller than that of the plane , as it is closer to 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 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 amount concentration on the upper face of AGV is illustrated in Figure 11. 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.

The necessity of the vertical external airflow was further proved by numerical experiments with boundaries described in Section . 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.

**(a)**

**(b)**

In this work, the concern about the proper setting of the speed of external vertical clean air under boundary conditions described in Section had been throughly studied. A series of concentration values of microcontaminants on surface of the cassette were obtained under various settings of , and the relationship between them is displayed in Figure 13.

From Figure 11, it is easy to see that with the increase of , the decreasing speed of contaminant concentration is slowing down; at , the concentration is below 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 elements (mesh size ) and analyzed parallelly by the scheme described in Section . 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, 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.

#### Acknowledgments

This work was supported by the National Science Foundation of China (NSFC), Grant 11202248, 91230114, and 11072272; the China Postdoctoral Science Foundation, Grant 2012M521646; and the Guangdong National Science Foundation, Grant S2012040007687.

#### References

- S. H. Tang, O. Motlagh, A. R. Ramli, N. Ismail, and D. N. Nia, “A novel GA-FCM strategy for motion learning and prediction: application in wireless tracking of intelligent subjects,”
*Arabian Journal for Science and Engineering*, vol. 37, no. 7, pp. 1929–1958, 2012. View at: Publisher Site | Google Scholar - S. E. Kesen and Ö. F. Baykoç, “Simulation of automated guided vehicle (AGV) systems based on just-in-time (JIT) philosophy in a job-shop environment,”
*Simulation Modelling Practice and Theory*, vol. 15, no. 3, pp. 272–284, 2007. View at: Publisher Site | Google Scholar - R. Ye, W. J. Hsu, and Z. H. Liu, “Web-based parallel simulation of AGVs using Java and JINI,” in
*Parallel Computing Technologies*, vol. 2127 of*Lecture Notes in Computer Science*, pp. 379–384, 2001. View at: Publisher Site | Google Scholar | Zentralblatt MATH - S. Rajotia, K. Shanker, and J. L. Batra, “Determination of optimal AGV fleet size for an FMS,”
*International Journal of Production Research*, vol. 36, no. 5, pp. 1177–1198, 1998. View at: Publisher Site | Google Scholar | Zentralblatt MATH - J. Lee, M. Tangjarukij, and Z. Zhu, “Load selection of automated guided vehicles in flexible manufacturing systems,”
*International Journal of Production Research*, vol. 34, no. 12, pp. 3383–3400, 1996. View at: Publisher Site | Google Scholar | Zentralblatt MATH - S. McClelland, “The cleanest robots in the world,”
*Industrial Robot*, vol. 13, no. 4, pp. 217–220, 1986. View at: Publisher Site | Google Scholar - L. Qiu, W. J. Hsu, S. Y. Huang, and H. Wang, “Scheduling and routing algorithms for AGVs: a survey,”
*International Journal of Production Research*, vol. 40, no. 3, pp. 745–760, 2002. View at: Publisher Site | Google Scholar | Zentralblatt MATH - S. Funabiki and K. Kodera, “A steering control of automated guided vehicles by the neural networks using a real-time tuning function,”
*IEEJ Transactions on Industry Applications*, vol. 118, no. 5, pp. 605–610, 1998. View at: Publisher Site | Google Scholar - C. W. Kim, J. M. A. Tanchoco, and P. H. Koo, “AGV dispatching based on workload balancing,”
*International Journal of Production Research*, vol. 37, no. 17, pp. 4053–4066, 1999. View at: Publisher Site | Google Scholar | Zentralblatt MATH - M. B. Zaremba, A. Obuchowicz, Z. A. Banaszak, and K. J. Jedrzejek, “A max-algebra approach to the robust distributed control of repetitive AGV systems,”
*International Journal of Production Research*, vol. 35, no. 10, pp. 2667–2687, 1997. View at: Publisher Site | Google Scholar - T. Nishi and Y. Tanaka, “Petri net decomposition approach for dispatching and conflict-free routing of bidirectional automated guided vehicle systems,”
*IEEE Transactions on Systems, Man and Cybernetics Part A*, vol. 42, no. 5, pp. 1230–1243, 2012. View at: Publisher Site | Google Scholar - C. H. Lin and L. L. Wang, “Intelligent collision avoidance by fuzzy logic control,”
*Robotics and Autonomous Systems*, vol. 20, no. 1, pp. 61–83, 1997. View at: Publisher Site | Google Scholar - C. Oboth, R. Batta, and M. Karwan, “Dynamic conflict-free routing of automated guided vehicles,”
*International Journal of Production Research*, vol. 37, no. 9, pp. 2003–2030, 1999. View at: Publisher Site | Google Scholar | Zentralblatt MATH - E. Sekine, Y. Hamamatsu, and T. Kongouji, “Analysis of traffic congestion in AGV system,”
*IEEJ Transactions on Industry Applications*, vol. 119, no. 4, pp. 515–522, 1999. View at: Google Scholar - S. J. Yang, “An arbitrary Lagrangian-Eulerian finite element method for interactions of airflow and a moving AGV in a cleanroom,”
*Finite Elements in Analysis and Design*, vol. 39, no. 5-6, pp. 521–533, 2003. View at: Publisher Site | Google Scholar - S. J. Yang, S. F. Chen, and W. S. Fu, “Numerical study of variations of airflow induced by a moving automatic guided vehicle in a cleanroom,”
*Journal of the Chinese Institute of Engineers*, vol. 25, no. 1, pp. 67–75, 2002. View at: Publisher Site | Google Scholar - A. G. Tannous, “Air flow simulation in a minienvironment,”
*Solid State Technology*, vol. 39, no. 7, pp. 201–209, 1996. View at: Publisher Site | Google Scholar - H. Kanayama, K. Toshigami, Y. Tashiro, M. Tabata, and S. Fujima, “Finite element analysis of air flow around an automatic guided vehicle,”
*Journal of Wind Engineering and Industrial Aerodynamics*, vol. 46-47, pp. 801–810, 1993. View at: Publisher Site | Google Scholar - T. Yamamoto and R. P. Donovan, “Model study for optimization of cleanroom airflows,”
*Journal of Environmental Sciences*, vol. 31, no. 6, pp. 24–29, 1988. View at: Google Scholar - T. H. Kuehn, “Computer simulation of airflow and particle transport in cleanrooms,”
*Journal of Environmental Sciences*, vol. 31, no. 5, pp. 21–27, 1988. View at: Google Scholar - J. Donea, S. Giuliani, and J. P. Halleux, “An arbitrary lagrangian-eulerian finite element method for transient dynamic fluid-structure interactions,”
*Computer Methods in Applied Mechanics and Engineering*, vol. 33, no. 1–3, pp. 689–723, 1982. View at: Publisher Site | Google Scholar | Zentralblatt MATH - H. Kanayama, H. Tsukikawa, and I. Ismail, “Simulation of hydrogen dispersion by the domain decomposition method,”
*Japan Journal of Industrial and Applied Mathematics*, vol. 28, pp. 43–53, 2011. View at: Publisher Site | Google Scholar | Zentralblatt MATH - Q. H. Yao, H. Kanayama, M. Ognio, and H. Notsu, “Incomplete balancing domain decomposition for large scale 3-D non-stationary incompressible flow problems,” in
*Proceedings of the 9th World Congress on Computational Mechanics and 4th Asian Pacific Congress on Computational Mechanics*, vol. 10, 2010. View at: Google Scholar - O. Pironneau,
*Finite Element Methods for Fluids*, John Wiley & Sons, Chichester, UK, 1989. View at: Zentralblatt MATH | MathSciNet - F. Brezzi and J. Douglas Jr., “Stabilized mixed methods for the Stokes problem,”
*Numerische Mathematik*, vol. 53, no. 1-2, pp. 225–235, 1988. View at: Publisher Site | Google Scholar | Zentralblatt MATH - J. Mandel, “Balancing domain decomposition,”
*Communications in Numerical Methods in Engineering*, vol. 9, no. 3, pp. 233–241, 1993. View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet - M. Ogino, R. Shioya, and H. Kanayama, “An inexact balancing preconditioner for large-scale structural analysis,”
*Journal of Computational Science and Technology*, vol. 2, no. 1, pp. 150–161, 2008. View at: Publisher Site | Google Scholar - F. Shoichi, T. Masahisa, and F. Yasuji, “Extension to three-dimensional problems of the upwind finite element scheme based on the choice of up- and downwind points,”
*Computer Methods in Applied Mechanics and Engineering*, vol. 112, no. 1–4, pp. 109–131, 1994. View at: Publisher Site | Google Scholar | Zentralblatt MATH - http://adventure.sys.t.u-tokyo.ac.jp/.

#### Copyright

Copyright © 2013 Qing-He Yao and Qing-Yong Zhu. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.