- Email: [email protected]
A new finite element approach to stress analysis in microfabrication technology
Microelectronics Reliability 41 (2001) 837±845
www.elsevier.com/locate/microrel
A new ®nite element approach to stress analysis in microfabrication technology q Slobodan Mijalkovic * Department of Microelectronics, Faculty of Information Technology and Systems, ECTM Laboratory-DIMES, Delft University of Technology, Mekelweg 4, 2628 CD Delft, The Netherlands Received 27 July 2000; received in revised form 1 October 2000
Abstract It is well known that the standard displacement-based ®nite element approach to viscoelastic stress analysis in microfabrication technology is susceptible to the volumetric locking and highly oscillatory pressure solutions. As an eective remedy for these problems, a new ®nite element approach based on separate weak statements for the displacement and pressure is proposed in this paper. Apart from the various existing techniques, the proposed approach is inherently free from the constraints imposed on the selection of discrete approximation spaces and permits a robust and stable discretisation based on piecewise linear Galerkin elements. The accuracy and stability of the proposed methodology has been tested in numerical experiments involving stress analysis problems from microfabrication deposition and native ®lm growing processes. Ó 2001 Elsevier Science Ltd. All rights reserved.
1. Introduction Microfabrication technology employs sequences of processes that utilise a variety of structural elements by embedding, butting and overlaying materials layers with dierent mechanical and electrical properties. Stress and strain are always present in such multilayer structures. Mechanical properties of the materials involved in the microfabrication are varying in the wide range from purely elastic solids to viscous ¯uids and therefore quite accurately modelled with the constitutive relationships of Maxwell viscoelasticity [1]. There is an increasing interest to simulate the stress and strain evolution in the microfabrication technology [2,3]. The principle aim is to predict and eliminate harmful eects that stress can impose on the produced microelectronics and micromechanical devices. Examples of stress related problems in defective devices could
q
An earlier version of this paper was printed at the First Small Systems Simulation Symposium (SSSS 2000), 4±5 September 2000, Nis, Yugoslavia. * Tel.: +31-15-278-7586; fax: +31-15-278-5691. E-mail address: [email protected] (S. Mijalkovic).
be the increased density of dislocations and point defects near silicon isolation regions, or the strain-induced bending of boron doped cantilevers in micromachining. Strong stress ®elds developed during microfabrication can severely aect the electrical characteristics and reliability of the produced devices. The accurate stress analysis could be useful in providing guidelines for the development of new technologies, that may take advantage of the stress presence. Namely, deliberately developed and carefully controlled stress and strain distributions can be used to enhance the performance of existing devices by tailoring the band structures, mobility, and material optical properties [4]. Finally, the stress analysis is essential for the accurate simulation of the microfabrication itself. In that sense, of special interest is the evaluation of the displacement ®elds during native ®lm growing processes [3]. Moreover, the stress ®elds modify diusivities and concentration of points defects aecting also the dopant redistribution during fabrication [5]. For many years, ®nite element method is the essential approach to the numerical stress analysis [6]. The attractive features of the method are its natural consistency with the variational formulation of the stress governing equations and the straightforward handling
0026-2714/01/$ - see front matter Ó 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 0 2 6 - 2 7 1 4 ( 0 1 ) 0 0 0 2 3 - 3
838
S. Mijalkovic / Microelectronics Reliability 41 (2001) 837±845
of geometrically complicated domains. However, it is recognised that the accuracy and stability of the standard displacement-based Galerkin ®nite element method could be seriously deteriorated in the approaching of the viscous ¯ow limit of viscoelastic materials. It is mainly contributed to the fact that the respective contributions of the dilatational and deviatoric stress components become in that case dierent by orders of magnitude. This phenomenon is also referred to as near incompressibility [6]. Particularly, the problem manifests itself as the nonrobustness of the displacement ®eld discrete approximation, broadly termed volumetric locking [7], and the unstable pressure discrete solution polluted by, so called, spurious pressure modes [8]. Various remedies have been proposed to overcome or reduce the eect of locking and pressure instability in the ®nite element implementation of the viscoelastic stress governing equations. Remaining in the framework of the displacement formulation the general approach is to apply less restrictive constraint on the incompressibility and thus reduce the ill-conditioned nature of the discrete equations. Perhaps the most widely accepted approach is to employ mixed methods based on dierent approximation spaces for the displacement and pressure [9]. However, the practical implementation of the mixed ®nite element methods is nontrivial and much more involved compared to standard piecewise linear elements. The main intention of this paper is to propose a method that eliminates the problem of the nearly incompressible materials already on the level of the weak formulation and allows a robust and stable implementation of piecewise linear elements for stress analysis problems in microfabrication. Some initial results concerning application of that new methodology to the problems of linear elasticity are presented in Ref. [10]. In this paper, the new ®nite element approach is extended to more general viscoelastic deformation stress models. The rest of the paper is organised as follows. Section 2 gives an overview of the stress governing equations with special emphasis on the near incompressibility problem. A new methodology for the weak statements formulation and ®nite element discretisation is presented in Section 3. It is applied in representative simulation examples in Section 4. Finally, a summary discussion is presented in Section 5.
2. Governing equations Consider a continuous material body occupying a bounded domain X with a piecewise smooth boundary C. The motion of the body is assumed to be governed by the momentum equation [11]. r rd rp
in X
1
and boundary conditions
pI rd n g on Cg ;
2
u h on Cu ;
3
where rd is the symmetric deviatoric stress tensor, p is the mean pressure, u is the displacement vector, g is the surface traction of the boundary segment Cg C, h is the displacement of the boundary segment Cu C (Cu \ Cg ;), and n is the outward unit normal vector on the boundary. Notice that acceleration terms and body forces are neglected in Eq. (1) which is quite justi®ed for the material deformation problems in microfabrication. The wide spectrum of the material deformation characteristics in microfabrication is commonly described by the constitutive relationships of Maxwell viscoelasticity. In the Maxwell's model, the elastic and viscous properties of the deviatoric stress component are combined by the dierential equation [1,11] 1 drd 1 dd rd ; 2G dt 2l dt
4
where d is the deviatoric strain tensor, G > 0 is the shear modulus and l > 0 is the material viscosity. Notice that Eq. (4) provides a continuous modelling of the material mechanical behaviour from the purely elastic deformation to the viscous ¯ow. Namely, for G>l it reduces to the Hooke's law for the elasticity while for l>G, we obtain the Newton law for the viscous ¯uids. For the dilatational stress component, associated with the pressure p, the viscous eects are neglected and it is described by [11] 1 dpI ds ; dt 3K dt
5
where s is the spherical strain tensor and K > 0 is the bulk modulus. In order to apply constitutive relationships of linear viscoelasticity to the large deformation problems in semiconductor processing, the total processing time is divided into many steps of smaller time periods. It is assumed that a small deformation strains are generates within a single time step Dt. This allows further simpli®cations of the Maxwell's model. Namely, for a short time period Dt, the rate of the strain tensor can be assumed to be constant in time. Consequently, the deviatoric stress tensor can be evaluated from Eq. (4) as [3] Dt rd 2Geff dd rd0 exp ; 6 s where Geff G
s 1 Dt
exp
Dt s
7
S. Mijalkovic / Microelectronics Reliability 41 (2001) 837±845
839
is the eective shear modulus, dd is the incremental strain tensor, rd0 is the total deviatoric stress tensor at the end of the previous time step and s
l G
8
is the Maxwellian relaxation time. It is also straightforward to evaluate the pressure from Eq. (5) as p
3Kds p0 ;
9
where ds is the incremental strain tensor and p0 is the pressure at the end of the previous time step. Since dd and ds are small deformation strain tensors, they can be quite accurately expressed by the ®rst order expressions [11]: dd 12 ru ruT ds 13 r uI;
1 r 3
uI;
10 11
where u is the incremental displacement ®eld. It is used to update the total strain ®eld at the end of each time step. Notice that with introduction of the eective shear modulus Geff , the large deformation problem of viscoelasticity is eectively reduced to the sequence of linear elasticity problems that should be solved at each time step. The principle diculty in the formulation of discrete models for the equations of linear elasticity appears in the case when the relative contribution of the deviatoric stress component is relatively small compared to the dilatation stress component. It is also referred to as the phenomenon of near incompressibility. The commonly used quantitative measure of incompressibility is the Poisson ratio m 0 6 m 6 0:5 [11]. The materials with m close to 0:5 are consider nearly incompressible. Although, the Poisson's ratio is safely away from the value 0.5 for most of the materials used in microfabrication, it is not the case when they are considered as material with eective rigidity Geff in incremental Maxwell viscoelasticity. In that case, it is possible to formulate an eective Poisson ratio meff
1 a ; 2a
12
where a
2Geff : 3K
Fig. 1. Eective Poisson ratio for silicon dioxide.
values of Dt=s, resulting from small values of s, which typically occurs at high processing temperatures, in combination with large time-step sizes, the eective Poisson's ratio asymptotically approaches the value 0.5 and silicon dioxide can be considered to be nearly incompressible. Besides the eective Poisson's ratio, the parameter a can be also used as a measure of eective material incompressibility. Notice from Eq. (13) that a 1 corresponds meff 0, while for a 0, which is equivalent to meff 0:5, we have an eectively incompressible material. It is found also more convenient to use a as a critical problem parameter in the description of the ®nite element implementation of the stress governing equations.
3. Finite element implementation In this section we will demonstrate how the problem of near incompressibility can be eliminated within the generalised Galerkin formulation of the stress governing equations. The formulation of the weak statement for the momentum equations is performed in a straightforward way. Multiplying the ith component of the momentum equation (1) by an arbitrary test function w, with w 0 on Cu , and integrating in X we obtain the residual statement Z Z r rd w rpw: 14 X
13
For Geff G, we have meff m. As an example, Fig. 1 shows the dependence of eective Poisson's ratio as a function of Dt=s for silicon dioxide. Notice that for large
X
Integrating by parts and taking into account the boundary condition (2) we obtain Z X
Z rd rw
X
Z prw
Cg
gw:
15
840
S. Mijalkovic / Microelectronics Reliability 41 (2001) 837±845
Z
Using Eqs. (6), (9) and (11), we ®nally have Z Z ru ruT rw 1 a p rw Geff X X Z Z Z Dt gw a p0 rw exp rd0 rw s Cg X X
X
16 as a weak Galerkin statement for the incremental displacement vector. The problems of nearly incompressible materials generally originate from the the pressure constitutive relationship (9) which can be also expressed as 3adp 2Geff r u 0;
17
where dp p p0 . Notice that when a ! 0, the pressure constitutive relationship degenerates into a divergencefree constraint r u 0. This fact makes the direct implementation of pressure from Eq. (17) in Eq. (16) inconvenient and require application of a special mixed weak statement formulated separately for Eq. (17). In order to avoid the diculties associated with the formulation of weak statements for Eq. (17), we can formulate an equivalent but robust weak statement starting again from the momentum equation. To this end, it is expressed as Geff r ru ruT 1 arp Dt a rp0 : exp 18 s Implementing the identity r ru ruT r r u 2r r u
19
into the divergence term of Eq. (18), multiplying Eq. (18) by rw, with w 0 on Cg , and integrating in X, we obtain Z Z 1 a rp rw 2Geff r r u rw ZX X Dt a rp0 rw Geff exp s X Z r r u rw: 20
rp rw
exp
Dt s
2a
Z
rp0 rw 1 2a X Z Geff r u n rw 1 2a Cu
22
as a weak statement for the pressure. Notice that Eq. (22) actually represents a weak statement for the Poisson type of equation Dt exp 2a 2 s r2 p r p0 in X; 23 1 2a and that the boundary integral term in Eq. (22) implements the appropriate natural pressure boundary condition. For a 0 and s Dt, Eq. (23) reduces to the well known pressure Poisson equation for the Stokes problem [12]. In that sense, Eq. (23) can be considered as generalisation of the pressure Poisson equation to more general problems of viscoelasticity. It should be emphasised that the weak statement (22) implicitly maintains the critical pressure constitutive relationship (17) in X. Namely, subtracting Eq. (22), after multiplication by 1 2a, and Eq. (20) we obtain a weak statement corresponding to the Laplace problem r2 3adp 2Geff r u 0 in X;
24
which means that Eq. (17) is maintained in X only if 3adp 2Geff r u 0 on C:
25
While Eq. (25) is implicitly applied with the natural boundary condition for pressure on Cu [12], it has to be enforced explicitly on Cg as an essential boundary condition for pressure. The diculties in practical implementation of the essential pressure boundary condition on Cg for a ! 0 can be easily circumvented by subtracting Eq. (25) and the normal component of the boundary condition (24) given by
X
Using integration by parts in the last integral term, it can be transformed into a boundary integral term as Z Z r r u rw r u n rw; 21 X
Cu
since r rw 0 in X for every test function w, and the vector r u n has only a tangential component on Cg where w 0. Employing now Eq. (17) in the second integral term of Eq. (20) we obtain
Fig. 2. Thin ®lm strips geometry and simulation domain.
S. Mijalkovic / Microelectronics Reliability 41 (2001) 837±845
1
ap 2Geff on un gn ap0
on Cg
26
841
where os us is the tangential derivative of the tangential incremental displacement vector on the boundary.
The weak statements (22) and (16) represent a generalised formulation of the problem and a basis for the ®nite element discretisation. It should be emphasised that the proposed methodology for the generalised formulation of the stress governing equations allow to employ equal order approximation spaces for the displacement and pressure including the simplest case of piecewise linear trial and test functions. In the simulation examples presented in Section 4, discretisation is performed on triangular grids with piecewise linear trial
Fig. 3. Pressure distribution: (a) a 0, (b) a 1.
Fig. 4. Lateral displacement: (a) a 0, (b) a 1.
where on un is the normal derivative of the normal incremental displacement on the boundary. For example, in 2D case we obtain as an essential boundary condition for the pressure 1 2ap 2Geff os us
gn 2ap0
on Cg
27
842
S. Mijalkovic / Microelectronics Reliability 41 (2001) 837±845
4.1. Film-edge induced stress analysis Modern integrated semiconductor fabrication technology relies on the use of thin ®lms for many purposes. For example, nitride ®lms are often patterned on the silicon substrates as isolation masks or as diusion barriers. Equidistant nitride strips patterned on the silicon substrate are also used in memory circuits. If the strips are narrow, high levels of stress are generated around the nitride mask regions. As a practical simulation example, we consider a material body ®lling the half space y P 0 and being covered with periodically repeating thin ®lm strips as shown in Fig. 2. The presence of thin ®lm strips with a built-in volumetric stress is physically associated with
Fig. 6. Pressure distribution near the ®lm edge for dierent mesh sizes.
Fig. 5. Vertical displacement: (a) a 0, (b) a 1.
and test functions for both displacement and pressure in weak statements (16) and (22). 4. Simulation examples The principal sources of the stress and strain ®elds in microfabrication are deposition of material layers with built-in intrinsic stress distribution and volumetric changes of natively growing material layers. Accordingly, we consider here simulation examples with the characteristic features of both stress sources.
Fig. 7. The cylindrical oxide domain.
S. Mijalkovic / Microelectronics Reliability 41 (2001) 837±845
concentrated line forces g acting tangentially to the body surface y 0. Because of the periodic structure of the ®lm strips, the re¯ecting boundary conditions are imposed along the boundary segments Cs . The challenging feature of the thin ®lm strip stress analysis problem is the singular behaviour of the displacement and pressure solutions near the ®lm edges. The typical discrete distributions of the lateral displacement component, vertical displacement components and pressure are shown in Figs. 3±5, respectively. Notice that the case of incompressible material a 0 is characterised with higher pressure, but lower lateral displacement values. Moreover, as can be seen in Fig. 4(a), the lateral displacement component has also negative values in the material bulk which compensate for the positive displacement of the material surface. Such negative lateral displacement values are not present in the compressible case shown in Fig. 4(b). The vertical displacement distribution shown in Fig. 5 is antisymmetric with maximum absolute values in the bulk, for a 0, and on the surface, for a 1. It has been veri®ed through a number of numerical experiments that the new ®nite element approach to stress analysis with linear basis functions for both displacement vector components and pressure, gives locking-free and stable discrete solutions in the whole range of the parameter a including the case of fully incompressible materials a 0. It should be emphasised that from the stability and accuracy point of view, for the proposed piecewise linear Galerkin approach the case a 1 or a 0 are fully equivalent to any other value of a in the range 0 6 a 6 1. In order to practically demonstrate simulation accuracy of the new ®nite element stress analysis formulation, Fig. 6 shows the comparison of the lateral distribution y 0 of the discrete pressure solutions which are obtained using dierent mesh sizes. The discrete solutions are also compared to an analytical solution which is possible to formulate in the case of semi-in®nite material samples exposed to concentrated
Fig. 8. The relative error of the displacement.
843
line forces [13]. The simulation is performed in fully incompressible case a 0 when standard displacement ®nite element formulation fails. Fig. 6 clearly demonstrates stable and uniform convergence of the discrete pressure solution near the near the thin ®lm edge with the concentrated line force where the analytical solution has a singular value. 4.2. Stress analysis in nonplanar thermal oxide The thermal oxidation processes produce oxide layers on the semiconductor surface to serve as a mask against
Fig. 9. Distribution of the x displacement component (a) standard displacement implementation, (b) a new implementation.
844
S. Mijalkovic / Microelectronics Reliability 41 (2001) 837±845
Fig. 10. Oxide domain and boundaries for local oxidation of silicon.
implant or diusion of dopant into silicon, to provide surface passivation, to isolate one device from another and to act as a component in MOS structures. As a ®rst practical example we consider here the thermal oxidation of the cylindrical (or ring) semiconductor structures. Fig. 7 shows an example of convex cylindrical domain used in the simulation. It is assumed that the chemical reaction that converts Si into SiO2 takes place at the oxide domain boundary Cu . Due to the volumetric expansion of the newly produced oxide, the boundary Cu is subject to a normal displacement h. The boundary Cg is assumed to be traction free while Cs denotes boundaries with re¯ecting boundary conditions. It should be emphasised that the cylindrical oxidesilicon structures are also often used for experimental calibration of the physical parameters for native growing thermal processes [2]. It also represents very useful
simulation example for simpli®ed analysis of more complex thermal oxidation topologies containing locally convex or concave geometrical shapes. As in the case of thin ®lm induced stress, the simple geometry provides also a possibility to construct an exact analytical solution for the displacement and pressure distributions in the cylindrically shaped oxide domains during thermal oxidation [2]. It is primarily used here for numerical investigation of the discretisation accuracy. To this end, Fig. 8 shows the distribution of the relative error of the discrete displacement along radial cylindrical coordinate. The comparison is made between the solutions obtained using standard displacement and new ®nite element simulations. The simulation is performed for the processing conditions corresponding to the incompressibility problem parameter a 0:01. Notice that the new ®nite element implementation of the stress governing equations produces several orders of magnitude smaller discretisation error. For even deeper levels of incompressibility, the nonrobustness of the displacement ®nite element implementation becomes even more emphasised. In order to practically demonstrate the nonrobustness of the displacement ®nite element implementation, Fig. 9 shows a the discrete solution for the x component of the displacement vector. The simulation is performed for the incompressibility parameter a 0:001. The second example is related to the local oxidation process (LOCOS) which is used for the semiconductor device isolation in integrated circuits. A typical LOCOS domain and boundaries are shown in Fig. 10. The boundaries Cu and Cg as well as normal displace-
Fig. 11. Distribution of pressure in the oxide during local oxidation of silicon (a) without boundary integral term in Eq. (22), (b) with boundary integral term in Eq. (22).
S. Mijalkovic / Microelectronics Reliability 41 (2001) 837±845
ment vector h have the same role as in the case of the cylindrical oxide structure. In a new ®nite element formulation of the stress governing equations, a special emphasise is put on the implementation of the correct Neumann boundary conditions for pressure [14]. It is actually implicitly incorporated into a weak statement (16). As an example Fig. 11 shows the distribution of pressure obtained without and with boundary integral term in Eq. (16) that accounts for Neumann pressure boundary conditions in the model problem. In both cases, a new ®nite element implementation of the stress governing equations has been employed and it is assumed that oxide is fully incompressible a 0. It should be emphasised that the usage of the inconsistent homogeneous Neumann pressure boundary condition in Eq. (22) results in the unphysical pressure ®eld shown in Fig. 11(a) and introduction of the numerical compressibility into simulation.
5. Conclusions A new implementation of the stress governing equations in semiconductor process simulation has been proposed in order to circumvent the problem of nearly incompressible materials. It is based on the formulation of two separate and irreducible weak statements for the displacement ®eld and pressure. In the strong formulation this approach could be also considered as a generalisation of the Poisson pressure equation to the wider class of stress analysis problems in semiconductor process simulation. The resulting Galerkin ®nite element discretisation enjoys stable and robust implementation of piecewise linear trial functions for both displacement and pressure. This signi®cantly simpli®es its practical implementation compared to the nonconforming, high order or mixed methods. A new implementation of the stress governing equations is tested in practical microfabrication simulation examples involving thin ®lm deposition and thermal oxidation.
845
References [1] Joppich W, Mijalkovic S. Semiconductor process modeling. In: Webster JG, editor. Wiley encyclopedia of electrical and electronics engineering, vol. 19. New York: Wiley; 1999. p 127±39. [2] Hu SM. Stress related problems in silicon technology. Int J Appl Phys 1991;15(6):R53±80. [3] Senez V, Collard D, Ferreira P, Baccus B. Two-dimensional simulation of local oxidation of silicon: calibrated viscoelastic ¯ow analysis. IEEE Trans Electron Dev 1996;43:720±31. [4] Jain SC. Germanium±silicon strained ®lms and heterostructures, Boston: Academic Press; 1994. [5] Park H, Jones KS, Slinkman JA, Law ME. Eect of hydrostatic pressure on dopant diusion in silicon. J Appl Phys 1995;78(6):3664±70. [6] Zienkiewicz OC, Taylor RL. The ®nite element methods: basic formulations and linear problems, vol. 1. London: McGraw-Hill; 1989. [7] Babuska I, Suri M. Locking eects in the ®nite element approximation of elasticity problems. Numer Math 1992;62:439±63. [8] Sani RL, Gresho PM, Lee RL, Griths DF. The cause and cure (?) of the spurious pressures generated by certain FEM solutions of the incompressible Navier± Stokes equations: Part 1. Int J Numer Meth Fluids 1981;1: 17-43. [9] Brezzi F, Fortin M. Mixed and hybrid ®nite element methods, New York: Springer; 1991. [10] Mijalkovic S. A piecewise linear galerkin approach to stress analysis of nearly incompressible materials. Commun Numer Meth Engng 2000;16(8):537±43. [11] Mase GE. Theory and problems of continuum mechanics, New York: McGraw-Hill; 1970. [12] Gresho PM. Incompressible ¯uid dynamics: some fundamental formulation issues. Annu Rev Fluid Mech 1991;21: 413±53. [13] Timoshenko SP, Goodier JN. Theory of elasticity, New York: McGraw-Hill; 1987. [14] Mijalkovic S, Joppich W. Derived boundary conditions for viscous thermal oxidation equations in pressure potential form, In: De Meyer K, Biesemans S, editors. Simulation of semiconductor processes and devices 1998, New York: Springer; 1998. p. 384±7.
www.elsevier.com/locate/microrel
A new ®nite element approach to stress analysis in microfabrication technology q Slobodan Mijalkovic * Department of Microelectronics, Faculty of Information Technology and Systems, ECTM Laboratory-DIMES, Delft University of Technology, Mekelweg 4, 2628 CD Delft, The Netherlands Received 27 July 2000; received in revised form 1 October 2000
Abstract It is well known that the standard displacement-based ®nite element approach to viscoelastic stress analysis in microfabrication technology is susceptible to the volumetric locking and highly oscillatory pressure solutions. As an eective remedy for these problems, a new ®nite element approach based on separate weak statements for the displacement and pressure is proposed in this paper. Apart from the various existing techniques, the proposed approach is inherently free from the constraints imposed on the selection of discrete approximation spaces and permits a robust and stable discretisation based on piecewise linear Galerkin elements. The accuracy and stability of the proposed methodology has been tested in numerical experiments involving stress analysis problems from microfabrication deposition and native ®lm growing processes. Ó 2001 Elsevier Science Ltd. All rights reserved.
1. Introduction Microfabrication technology employs sequences of processes that utilise a variety of structural elements by embedding, butting and overlaying materials layers with dierent mechanical and electrical properties. Stress and strain are always present in such multilayer structures. Mechanical properties of the materials involved in the microfabrication are varying in the wide range from purely elastic solids to viscous ¯uids and therefore quite accurately modelled with the constitutive relationships of Maxwell viscoelasticity [1]. There is an increasing interest to simulate the stress and strain evolution in the microfabrication technology [2,3]. The principle aim is to predict and eliminate harmful eects that stress can impose on the produced microelectronics and micromechanical devices. Examples of stress related problems in defective devices could
q
An earlier version of this paper was printed at the First Small Systems Simulation Symposium (SSSS 2000), 4±5 September 2000, Nis, Yugoslavia. * Tel.: +31-15-278-7586; fax: +31-15-278-5691. E-mail address: [email protected] (S. Mijalkovic).
be the increased density of dislocations and point defects near silicon isolation regions, or the strain-induced bending of boron doped cantilevers in micromachining. Strong stress ®elds developed during microfabrication can severely aect the electrical characteristics and reliability of the produced devices. The accurate stress analysis could be useful in providing guidelines for the development of new technologies, that may take advantage of the stress presence. Namely, deliberately developed and carefully controlled stress and strain distributions can be used to enhance the performance of existing devices by tailoring the band structures, mobility, and material optical properties [4]. Finally, the stress analysis is essential for the accurate simulation of the microfabrication itself. In that sense, of special interest is the evaluation of the displacement ®elds during native ®lm growing processes [3]. Moreover, the stress ®elds modify diusivities and concentration of points defects aecting also the dopant redistribution during fabrication [5]. For many years, ®nite element method is the essential approach to the numerical stress analysis [6]. The attractive features of the method are its natural consistency with the variational formulation of the stress governing equations and the straightforward handling
0026-2714/01/$ - see front matter Ó 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 0 2 6 - 2 7 1 4 ( 0 1 ) 0 0 0 2 3 - 3
838
S. Mijalkovic / Microelectronics Reliability 41 (2001) 837±845
of geometrically complicated domains. However, it is recognised that the accuracy and stability of the standard displacement-based Galerkin ®nite element method could be seriously deteriorated in the approaching of the viscous ¯ow limit of viscoelastic materials. It is mainly contributed to the fact that the respective contributions of the dilatational and deviatoric stress components become in that case dierent by orders of magnitude. This phenomenon is also referred to as near incompressibility [6]. Particularly, the problem manifests itself as the nonrobustness of the displacement ®eld discrete approximation, broadly termed volumetric locking [7], and the unstable pressure discrete solution polluted by, so called, spurious pressure modes [8]. Various remedies have been proposed to overcome or reduce the eect of locking and pressure instability in the ®nite element implementation of the viscoelastic stress governing equations. Remaining in the framework of the displacement formulation the general approach is to apply less restrictive constraint on the incompressibility and thus reduce the ill-conditioned nature of the discrete equations. Perhaps the most widely accepted approach is to employ mixed methods based on dierent approximation spaces for the displacement and pressure [9]. However, the practical implementation of the mixed ®nite element methods is nontrivial and much more involved compared to standard piecewise linear elements. The main intention of this paper is to propose a method that eliminates the problem of the nearly incompressible materials already on the level of the weak formulation and allows a robust and stable implementation of piecewise linear elements for stress analysis problems in microfabrication. Some initial results concerning application of that new methodology to the problems of linear elasticity are presented in Ref. [10]. In this paper, the new ®nite element approach is extended to more general viscoelastic deformation stress models. The rest of the paper is organised as follows. Section 2 gives an overview of the stress governing equations with special emphasis on the near incompressibility problem. A new methodology for the weak statements formulation and ®nite element discretisation is presented in Section 3. It is applied in representative simulation examples in Section 4. Finally, a summary discussion is presented in Section 5.
2. Governing equations Consider a continuous material body occupying a bounded domain X with a piecewise smooth boundary C. The motion of the body is assumed to be governed by the momentum equation [11]. r rd rp
in X
1
and boundary conditions
pI rd n g on Cg ;
2
u h on Cu ;
3
where rd is the symmetric deviatoric stress tensor, p is the mean pressure, u is the displacement vector, g is the surface traction of the boundary segment Cg C, h is the displacement of the boundary segment Cu C (Cu \ Cg ;), and n is the outward unit normal vector on the boundary. Notice that acceleration terms and body forces are neglected in Eq. (1) which is quite justi®ed for the material deformation problems in microfabrication. The wide spectrum of the material deformation characteristics in microfabrication is commonly described by the constitutive relationships of Maxwell viscoelasticity. In the Maxwell's model, the elastic and viscous properties of the deviatoric stress component are combined by the dierential equation [1,11] 1 drd 1 dd rd ; 2G dt 2l dt
4
where d is the deviatoric strain tensor, G > 0 is the shear modulus and l > 0 is the material viscosity. Notice that Eq. (4) provides a continuous modelling of the material mechanical behaviour from the purely elastic deformation to the viscous ¯ow. Namely, for G>l it reduces to the Hooke's law for the elasticity while for l>G, we obtain the Newton law for the viscous ¯uids. For the dilatational stress component, associated with the pressure p, the viscous eects are neglected and it is described by [11] 1 dpI ds ; dt 3K dt
5
where s is the spherical strain tensor and K > 0 is the bulk modulus. In order to apply constitutive relationships of linear viscoelasticity to the large deformation problems in semiconductor processing, the total processing time is divided into many steps of smaller time periods. It is assumed that a small deformation strains are generates within a single time step Dt. This allows further simpli®cations of the Maxwell's model. Namely, for a short time period Dt, the rate of the strain tensor can be assumed to be constant in time. Consequently, the deviatoric stress tensor can be evaluated from Eq. (4) as [3] Dt rd 2Geff dd rd0 exp ; 6 s where Geff G
s 1 Dt
exp
Dt s
7
S. Mijalkovic / Microelectronics Reliability 41 (2001) 837±845
839
is the eective shear modulus, dd is the incremental strain tensor, rd0 is the total deviatoric stress tensor at the end of the previous time step and s
l G
8
is the Maxwellian relaxation time. It is also straightforward to evaluate the pressure from Eq. (5) as p
3Kds p0 ;
9
where ds is the incremental strain tensor and p0 is the pressure at the end of the previous time step. Since dd and ds are small deformation strain tensors, they can be quite accurately expressed by the ®rst order expressions [11]: dd 12 ru ruT ds 13 r uI;
1 r 3
uI;
10 11
where u is the incremental displacement ®eld. It is used to update the total strain ®eld at the end of each time step. Notice that with introduction of the eective shear modulus Geff , the large deformation problem of viscoelasticity is eectively reduced to the sequence of linear elasticity problems that should be solved at each time step. The principle diculty in the formulation of discrete models for the equations of linear elasticity appears in the case when the relative contribution of the deviatoric stress component is relatively small compared to the dilatation stress component. It is also referred to as the phenomenon of near incompressibility. The commonly used quantitative measure of incompressibility is the Poisson ratio m 0 6 m 6 0:5 [11]. The materials with m close to 0:5 are consider nearly incompressible. Although, the Poisson's ratio is safely away from the value 0.5 for most of the materials used in microfabrication, it is not the case when they are considered as material with eective rigidity Geff in incremental Maxwell viscoelasticity. In that case, it is possible to formulate an eective Poisson ratio meff
1 a ; 2a
12
where a
2Geff : 3K
Fig. 1. Eective Poisson ratio for silicon dioxide.
values of Dt=s, resulting from small values of s, which typically occurs at high processing temperatures, in combination with large time-step sizes, the eective Poisson's ratio asymptotically approaches the value 0.5 and silicon dioxide can be considered to be nearly incompressible. Besides the eective Poisson's ratio, the parameter a can be also used as a measure of eective material incompressibility. Notice from Eq. (13) that a 1 corresponds meff 0, while for a 0, which is equivalent to meff 0:5, we have an eectively incompressible material. It is found also more convenient to use a as a critical problem parameter in the description of the ®nite element implementation of the stress governing equations.
3. Finite element implementation In this section we will demonstrate how the problem of near incompressibility can be eliminated within the generalised Galerkin formulation of the stress governing equations. The formulation of the weak statement for the momentum equations is performed in a straightforward way. Multiplying the ith component of the momentum equation (1) by an arbitrary test function w, with w 0 on Cu , and integrating in X we obtain the residual statement Z Z r rd w rpw: 14 X
13
For Geff G, we have meff m. As an example, Fig. 1 shows the dependence of eective Poisson's ratio as a function of Dt=s for silicon dioxide. Notice that for large
X
Integrating by parts and taking into account the boundary condition (2) we obtain Z X
Z rd rw
X
Z prw
Cg
gw:
15
840
S. Mijalkovic / Microelectronics Reliability 41 (2001) 837±845
Z
Using Eqs. (6), (9) and (11), we ®nally have Z Z ru ruT rw 1 a p rw Geff X X Z Z Z Dt gw a p0 rw exp rd0 rw s Cg X X
X
16 as a weak Galerkin statement for the incremental displacement vector. The problems of nearly incompressible materials generally originate from the the pressure constitutive relationship (9) which can be also expressed as 3adp 2Geff r u 0;
17
where dp p p0 . Notice that when a ! 0, the pressure constitutive relationship degenerates into a divergencefree constraint r u 0. This fact makes the direct implementation of pressure from Eq. (17) in Eq. (16) inconvenient and require application of a special mixed weak statement formulated separately for Eq. (17). In order to avoid the diculties associated with the formulation of weak statements for Eq. (17), we can formulate an equivalent but robust weak statement starting again from the momentum equation. To this end, it is expressed as Geff r ru ruT 1 arp Dt a rp0 : exp 18 s Implementing the identity r ru ruT r r u 2r r u
19
into the divergence term of Eq. (18), multiplying Eq. (18) by rw, with w 0 on Cg , and integrating in X, we obtain Z Z 1 a rp rw 2Geff r r u rw ZX X Dt a rp0 rw Geff exp s X Z r r u rw: 20
rp rw
exp
Dt s
2a
Z
rp0 rw 1 2a X Z Geff r u n rw 1 2a Cu
22
as a weak statement for the pressure. Notice that Eq. (22) actually represents a weak statement for the Poisson type of equation Dt exp 2a 2 s r2 p r p0 in X; 23 1 2a and that the boundary integral term in Eq. (22) implements the appropriate natural pressure boundary condition. For a 0 and s Dt, Eq. (23) reduces to the well known pressure Poisson equation for the Stokes problem [12]. In that sense, Eq. (23) can be considered as generalisation of the pressure Poisson equation to more general problems of viscoelasticity. It should be emphasised that the weak statement (22) implicitly maintains the critical pressure constitutive relationship (17) in X. Namely, subtracting Eq. (22), after multiplication by 1 2a, and Eq. (20) we obtain a weak statement corresponding to the Laplace problem r2 3adp 2Geff r u 0 in X;
24
which means that Eq. (17) is maintained in X only if 3adp 2Geff r u 0 on C:
25
While Eq. (25) is implicitly applied with the natural boundary condition for pressure on Cu [12], it has to be enforced explicitly on Cg as an essential boundary condition for pressure. The diculties in practical implementation of the essential pressure boundary condition on Cg for a ! 0 can be easily circumvented by subtracting Eq. (25) and the normal component of the boundary condition (24) given by
X
Using integration by parts in the last integral term, it can be transformed into a boundary integral term as Z Z r r u rw r u n rw; 21 X
Cu
since r rw 0 in X for every test function w, and the vector r u n has only a tangential component on Cg where w 0. Employing now Eq. (17) in the second integral term of Eq. (20) we obtain
Fig. 2. Thin ®lm strips geometry and simulation domain.
S. Mijalkovic / Microelectronics Reliability 41 (2001) 837±845
1
ap 2Geff on un gn ap0
on Cg
26
841
where os us is the tangential derivative of the tangential incremental displacement vector on the boundary.
The weak statements (22) and (16) represent a generalised formulation of the problem and a basis for the ®nite element discretisation. It should be emphasised that the proposed methodology for the generalised formulation of the stress governing equations allow to employ equal order approximation spaces for the displacement and pressure including the simplest case of piecewise linear trial and test functions. In the simulation examples presented in Section 4, discretisation is performed on triangular grids with piecewise linear trial
Fig. 3. Pressure distribution: (a) a 0, (b) a 1.
Fig. 4. Lateral displacement: (a) a 0, (b) a 1.
where on un is the normal derivative of the normal incremental displacement on the boundary. For example, in 2D case we obtain as an essential boundary condition for the pressure 1 2ap 2Geff os us
gn 2ap0
on Cg
27
842
S. Mijalkovic / Microelectronics Reliability 41 (2001) 837±845
4.1. Film-edge induced stress analysis Modern integrated semiconductor fabrication technology relies on the use of thin ®lms for many purposes. For example, nitride ®lms are often patterned on the silicon substrates as isolation masks or as diusion barriers. Equidistant nitride strips patterned on the silicon substrate are also used in memory circuits. If the strips are narrow, high levels of stress are generated around the nitride mask regions. As a practical simulation example, we consider a material body ®lling the half space y P 0 and being covered with periodically repeating thin ®lm strips as shown in Fig. 2. The presence of thin ®lm strips with a built-in volumetric stress is physically associated with
Fig. 6. Pressure distribution near the ®lm edge for dierent mesh sizes.
Fig. 5. Vertical displacement: (a) a 0, (b) a 1.
and test functions for both displacement and pressure in weak statements (16) and (22). 4. Simulation examples The principal sources of the stress and strain ®elds in microfabrication are deposition of material layers with built-in intrinsic stress distribution and volumetric changes of natively growing material layers. Accordingly, we consider here simulation examples with the characteristic features of both stress sources.
Fig. 7. The cylindrical oxide domain.
S. Mijalkovic / Microelectronics Reliability 41 (2001) 837±845
concentrated line forces g acting tangentially to the body surface y 0. Because of the periodic structure of the ®lm strips, the re¯ecting boundary conditions are imposed along the boundary segments Cs . The challenging feature of the thin ®lm strip stress analysis problem is the singular behaviour of the displacement and pressure solutions near the ®lm edges. The typical discrete distributions of the lateral displacement component, vertical displacement components and pressure are shown in Figs. 3±5, respectively. Notice that the case of incompressible material a 0 is characterised with higher pressure, but lower lateral displacement values. Moreover, as can be seen in Fig. 4(a), the lateral displacement component has also negative values in the material bulk which compensate for the positive displacement of the material surface. Such negative lateral displacement values are not present in the compressible case shown in Fig. 4(b). The vertical displacement distribution shown in Fig. 5 is antisymmetric with maximum absolute values in the bulk, for a 0, and on the surface, for a 1. It has been veri®ed through a number of numerical experiments that the new ®nite element approach to stress analysis with linear basis functions for both displacement vector components and pressure, gives locking-free and stable discrete solutions in the whole range of the parameter a including the case of fully incompressible materials a 0. It should be emphasised that from the stability and accuracy point of view, for the proposed piecewise linear Galerkin approach the case a 1 or a 0 are fully equivalent to any other value of a in the range 0 6 a 6 1. In order to practically demonstrate simulation accuracy of the new ®nite element stress analysis formulation, Fig. 6 shows the comparison of the lateral distribution y 0 of the discrete pressure solutions which are obtained using dierent mesh sizes. The discrete solutions are also compared to an analytical solution which is possible to formulate in the case of semi-in®nite material samples exposed to concentrated
Fig. 8. The relative error of the displacement.
843
line forces [13]. The simulation is performed in fully incompressible case a 0 when standard displacement ®nite element formulation fails. Fig. 6 clearly demonstrates stable and uniform convergence of the discrete pressure solution near the near the thin ®lm edge with the concentrated line force where the analytical solution has a singular value. 4.2. Stress analysis in nonplanar thermal oxide The thermal oxidation processes produce oxide layers on the semiconductor surface to serve as a mask against
Fig. 9. Distribution of the x displacement component (a) standard displacement implementation, (b) a new implementation.
844
S. Mijalkovic / Microelectronics Reliability 41 (2001) 837±845
Fig. 10. Oxide domain and boundaries for local oxidation of silicon.
implant or diusion of dopant into silicon, to provide surface passivation, to isolate one device from another and to act as a component in MOS structures. As a ®rst practical example we consider here the thermal oxidation of the cylindrical (or ring) semiconductor structures. Fig. 7 shows an example of convex cylindrical domain used in the simulation. It is assumed that the chemical reaction that converts Si into SiO2 takes place at the oxide domain boundary Cu . Due to the volumetric expansion of the newly produced oxide, the boundary Cu is subject to a normal displacement h. The boundary Cg is assumed to be traction free while Cs denotes boundaries with re¯ecting boundary conditions. It should be emphasised that the cylindrical oxidesilicon structures are also often used for experimental calibration of the physical parameters for native growing thermal processes [2]. It also represents very useful
simulation example for simpli®ed analysis of more complex thermal oxidation topologies containing locally convex or concave geometrical shapes. As in the case of thin ®lm induced stress, the simple geometry provides also a possibility to construct an exact analytical solution for the displacement and pressure distributions in the cylindrically shaped oxide domains during thermal oxidation [2]. It is primarily used here for numerical investigation of the discretisation accuracy. To this end, Fig. 8 shows the distribution of the relative error of the discrete displacement along radial cylindrical coordinate. The comparison is made between the solutions obtained using standard displacement and new ®nite element simulations. The simulation is performed for the processing conditions corresponding to the incompressibility problem parameter a 0:01. Notice that the new ®nite element implementation of the stress governing equations produces several orders of magnitude smaller discretisation error. For even deeper levels of incompressibility, the nonrobustness of the displacement ®nite element implementation becomes even more emphasised. In order to practically demonstrate the nonrobustness of the displacement ®nite element implementation, Fig. 9 shows a the discrete solution for the x component of the displacement vector. The simulation is performed for the incompressibility parameter a 0:001. The second example is related to the local oxidation process (LOCOS) which is used for the semiconductor device isolation in integrated circuits. A typical LOCOS domain and boundaries are shown in Fig. 10. The boundaries Cu and Cg as well as normal displace-
Fig. 11. Distribution of pressure in the oxide during local oxidation of silicon (a) without boundary integral term in Eq. (22), (b) with boundary integral term in Eq. (22).
S. Mijalkovic / Microelectronics Reliability 41 (2001) 837±845
ment vector h have the same role as in the case of the cylindrical oxide structure. In a new ®nite element formulation of the stress governing equations, a special emphasise is put on the implementation of the correct Neumann boundary conditions for pressure [14]. It is actually implicitly incorporated into a weak statement (16). As an example Fig. 11 shows the distribution of pressure obtained without and with boundary integral term in Eq. (16) that accounts for Neumann pressure boundary conditions in the model problem. In both cases, a new ®nite element implementation of the stress governing equations has been employed and it is assumed that oxide is fully incompressible a 0. It should be emphasised that the usage of the inconsistent homogeneous Neumann pressure boundary condition in Eq. (22) results in the unphysical pressure ®eld shown in Fig. 11(a) and introduction of the numerical compressibility into simulation.
5. Conclusions A new implementation of the stress governing equations in semiconductor process simulation has been proposed in order to circumvent the problem of nearly incompressible materials. It is based on the formulation of two separate and irreducible weak statements for the displacement ®eld and pressure. In the strong formulation this approach could be also considered as a generalisation of the Poisson pressure equation to the wider class of stress analysis problems in semiconductor process simulation. The resulting Galerkin ®nite element discretisation enjoys stable and robust implementation of piecewise linear trial functions for both displacement and pressure. This signi®cantly simpli®es its practical implementation compared to the nonconforming, high order or mixed methods. A new implementation of the stress governing equations is tested in practical microfabrication simulation examples involving thin ®lm deposition and thermal oxidation.
845
References [1] Joppich W, Mijalkovic S. Semiconductor process modeling. In: Webster JG, editor. Wiley encyclopedia of electrical and electronics engineering, vol. 19. New York: Wiley; 1999. p 127±39. [2] Hu SM. Stress related problems in silicon technology. Int J Appl Phys 1991;15(6):R53±80. [3] Senez V, Collard D, Ferreira P, Baccus B. Two-dimensional simulation of local oxidation of silicon: calibrated viscoelastic ¯ow analysis. IEEE Trans Electron Dev 1996;43:720±31. [4] Jain SC. Germanium±silicon strained ®lms and heterostructures, Boston: Academic Press; 1994. [5] Park H, Jones KS, Slinkman JA, Law ME. Eect of hydrostatic pressure on dopant diusion in silicon. J Appl Phys 1995;78(6):3664±70. [6] Zienkiewicz OC, Taylor RL. The ®nite element methods: basic formulations and linear problems, vol. 1. London: McGraw-Hill; 1989. [7] Babuska I, Suri M. Locking eects in the ®nite element approximation of elasticity problems. Numer Math 1992;62:439±63. [8] Sani RL, Gresho PM, Lee RL, Griths DF. The cause and cure (?) of the spurious pressures generated by certain FEM solutions of the incompressible Navier± Stokes equations: Part 1. Int J Numer Meth Fluids 1981;1: 17-43. [9] Brezzi F, Fortin M. Mixed and hybrid ®nite element methods, New York: Springer; 1991. [10] Mijalkovic S. A piecewise linear galerkin approach to stress analysis of nearly incompressible materials. Commun Numer Meth Engng 2000;16(8):537±43. [11] Mase GE. Theory and problems of continuum mechanics, New York: McGraw-Hill; 1970. [12] Gresho PM. Incompressible ¯uid dynamics: some fundamental formulation issues. Annu Rev Fluid Mech 1991;21: 413±53. [13] Timoshenko SP, Goodier JN. Theory of elasticity, New York: McGraw-Hill; 1987. [14] Mijalkovic S, Joppich W. Derived boundary conditions for viscous thermal oxidation equations in pressure potential form, In: De Meyer K, Biesemans S, editors. Simulation of semiconductor processes and devices 1998, New York: Springer; 1998. p. 384±7.