변분법을 응용한 초크랄스키 공정에서 단결정 표면 온도 분포의 최적화.

Part I : 열 응력 문제

정 자훈, 강 인석

포항 공과 대학교 화학 공학과 및 기계/산업 공학부

Application of variational calculus to optimization of single crystal

surface temperature distribution in the Czochralski process.

Part I : Thermal stress problem

J.H. Jeong, I.S. Kang

Department of Chemical Engineering and Division of Mechanical/Industrial Engineering, Pohang University of Science and Technology

INTRODUCTION

In the Czochralski(CZ) process, the quality of single crystal is determined by crystalline defects. Since microscopic phenomena inside single crystal are strongly influenced by thermal conditions, many researches have aimed at optimal design of thermal surroundings. In particular, the increases of capabilities of numerical analysis have stimulated numerical computations to calculate the temperature distribution and to predict single crystal qualities. However, most numerical works for the design of processing conditions have been accomplished based on the simulations case by case [1, 2]. In this work, we want to obtain the optimal distribution of the crystal surface temperature in a systematic way.

In our optimization problem, we seek an optimal solution in the sense that the index of the crystal quality is maximized while the growth rate is also maximized. The objective function is formulated by considering the distribution of the von Mises stress because it is believed to provide the driving forces that induce crystalline dislocations. The radial uniformity of the temperature distribution in the crystal is also included in the objective function because it is directly related to the uniform properties of produced wafers. As the single crystal growth rate, the heat transfer rate at the crystal-melt interface toward the crystal phase is considered. In order to solve the optimization problem with the equality constraints described by partial differential equations, we apply the variational calculus and solve the Euler-Lagrange equations, which are given by the coupled partial differential equations, by using the iterative numerical scheme proposed in this work.

The optimal distributions of the crystal surface temperature obtained in this work may provide an insight for the optimal design of thermal surroundings such as the configurations of thermal shield and the heater/cooler positions. Futhermore, the methodology adopted in this work can be applied to other optimization problems for which state equations are given by partial differential equations.

VARIATIONAL FORMULATION

The optimization problem is given by

[eqn: {rm Min.} I ~=~ alpha_u int_{Omega} F_u {rm d} Omega + alpha_{rm S} int_{Omega} F_{rm S} {rm d}Omega - alpha_G ``G],	(1)

where [eqn: alpha_u], [eqn: alpha_{rm S}], [eqn: alpha_G] are positive weighting parameters and [eqn: F_u], [eqn: F_{rm S}], [eqn: G] are the performance measures that represent the radial uniformity of the temperature distribution, the magnitude of the von Mises stress, and the single crystal growth rate respectively. Thus we have

[eqn: F_u = left( {partial u}/{partial r} right)^2],  [eqn: F_{rm S} = left( tau_{vm} right)^2],  [eqn: G = int_{Gamma_1} ( {bold {rm n}} cdot nabla u ) {rm d} Gamma],	(2)

where [eqn: u] denotes by the dimensionless temperature in the crystal phase, [eqn: {bold {rm n}}] the outgoing unit normal vector. The von Mises stress, [eqn: tau_{v``m}] is defined by

[eqn: tau_{vm} = left( ``(1/2) ````({bold {rm S}}_d ````: {bold {rm S}}_d ``)} right)^{1/2}],  where[eqn: {bold {rm S}}_d equiv {bold {rm S}} - (1/3) ( {bold {rm S}}:{bold {rm E}}) {bold {rm E}}] ,	(3)

and [eqn: {rm {bold E}}] is the identity tensor.

As shown in Fig. 1, the crystal is assumed to by a straight cylinder. The temperature field is governed by the energy balance equation and the stress field is governed by Hooks' law and the force balance equations in static equilibrium [3]. Their dimensionless forms are given by

[eqn: {rm L}_1 ({bold {rm D}}, {rm {bold v}} ) = {bold {rm D}} - {1}over{2} left( (nabla {bold {rm v}}) + (nabla {bold {rm v}})^T right) = 0],	(4)

[eqn: {rm L}_2 ({bold {rm S}}, {rm {bold D}}, u`` ) = {bold {rm S}} - {{rm {bold D}}}over{(1+ nu)} - {nu ``({bold {rm D}}:{bold {rm E}}){bold {rm E}}}over{(1+nu)(1-2 nu)} + {beta`` (u-1){bold {rm E}}}over{(1-2 nu)}],	(5)

[eqn: {rm L}_3 ({bold {rm S}} ) = nabla cdot {bold {rm S}} = 0],	(6)

[eqn: {rm L}_4 ( u`` ) = nabla^2 u = 0],	(7)

where [eqn: L_1], [eqn: L_2], [eqn: L_3], [eqn: L_4] represent differential or algebraic operators,[eqn: {bold {rm v}}] the displacement vector, [eqn: {bold {rm D}}] the strain tensor, [eqn: nu``] the Poisson ratio, and [eqn: beta``] the dimensionless thermal expansion coefficient. When the data given in the literature are used, [eqn: nu = 0.25] and [eqn: beta=0.01] in the silicon single crystal growth [1]. Since the displacements are much less than unity, their effect on the system geometry is neglected. The gravitational force is also neglected in the calculation of the stress fields. We have the traction free conditions for all boundaries since the top and bottom plates and the side of the crystal phase are free.

In order to handle the equality constraints of (4)-(7), we introduce the augmented objective function of [eqn: I_a = I + I_{rm S} + I_{rm D} + I_{rm v} + I_u] , where

[eqn: I_{rm S} = < {bold {rm S}}^* , {rm L}_1 ({bold {rm D}}, {bold {rm v}} ) >],	[eqn: I_{rm D} = < {bold {rm D}}^* , {rm L}_2 ({bold {rm S}}, {rm {bold D}}, u`` ) >],

[eqn: I_{rm v} = < {bold {rm v}}^* , {rm L}_3 ({bold {rm S}} ``) >],	[eqn: I_{u} = < {u}^* , {rm L}_4 (u ) >].

Here [eqn: {bold {rm S}}^*], [eqn: {bold {rm D}}^*], [eqn: {bold {rm v}}^*], [eqn: {u}^*] are Lagrange multipliers, called adjoint variables.

As system of equations, Euler-Lagrange equations can be obtained from the variations of the augmented objective function. By taking the variations to [eqn: I_a] with respect to adjoint variables, we recover the state equations of (4)-(7). On the other hand, by taking the variations to [eqn: I_a] with respect to state variables, we have the adjoint equations, given by

[eqn: {bold {rm D}}^* - {1}over{2} left( (nabla {bold {rm v}}^* ) + (nabla {bold {rm v}}^* )^T ``right) + alpha_{rm S} left( {bold {rm S}} - {1}over{3} ( {bold {rm S}} : {bold {rm E}}) {bold {rm E}} right) = 0] ,	(8)

[eqn: {bold {rm S} }^* - {1}over{1+ nu} {bold {rm D}}^* - {nu}over{(1+ nu)(1-2 nu)}( {bold {rm D}}^* : {bold {rm E}}) {bold {rm E}} = 0] ,	(9)

[eqn: nabla cdot {bold {rm S}}^* = 0] ,	(10)

[eqn: nabla^2 {u}^* + {beta}over{(1-2 nu)} ({bold {rm D}}^* : {bold {rm E}}) + 2 alpha_u left( {partial^2 u}over{partial r^2} + {1}over{r} {partial u}over{partial r} right) = 0] .	(11)

In this work, we use the displacement formulation, where only [eqn: u],[eqn: {bold {rm v}}], [eqn: u^*],[eqn: {bold {rm v}}^*] appear in the mathematical boundary-value problem, whereas the remainder are determined afterwards by direct calculation.

RESULTS AND DISCUSSION

In the cylindrical coordinates with axisymmetric condition, we have six coupled partial differential equations for three state variables of [eqn: u], [eqn: v_z], [eqn: v_r] and three adjoint variables of [eqn: u^*], [eqn: v_z^*], [eqn: v_r^*]. In solving these equations, we use finite difference method and alternative iterative scheme. Fig. 2 shows the convergence of the numerical scheme when [eqn: alpha_u =1.0], [eqn: alpha_{rm S} = 20.0], [eqn: alpha_G = 0.1]. In the figure, the numerical errors for [eqn: u], [eqn: v_z], [eqn: v_r], [eqn: u^*], [eqn: v_z^*], [eqn: v_r^*] are denoted by [eqn: e_1], [eqn: e_2], [eqn: e_3], [eqn: e_4], [eqn: e_5], [eqn: e_6] respectively. From the figure, we can see that iterative calculations make all variables converge competitively each other. A prototype of the result for the optimal solutions is shown in Fig. 3, where the gray regions mean that the variables have negative values. It is interesting that the adjoint variables have similar spatial distributions with their corresponding state variables as if mirror images.

Fig. 4 shows the effects of [eqn: alpha_{rm S}] on the optimal surface temperature profiles when [eqn: alpha_u = 0.0], [eqn: alpha_G = 0.1]. It is noteworthy that the slope of the axial temperature gradient at the interface (z=0) becomes steeper as the weighting parameter [eqn: alpha_{rm S}] increases. Since the crystal growth rate is directly proportional to the conductive heat transfer at the interface, we can see that the crystal growth rate decreases as [eqn: alpha_{rm S}] increases. In other words, in order to reduce the level of thermal stress in the crystal phase, we have to maintain low crystal growth rate. On the other hand, Fig. 5 shows the effects of [eqn: alpha_u] when [eqn: alpha_{rm S} = 0.0], [eqn: alpha_G = 0.1]. It is noteworthy that effects of [eqn: alpha_u] are very similar to those of [eqn: alpha_{rm S}]. In fact, this qualitative similarity between them can be expected by deriving asymptotic solution in the limiting case of either [eqn: alpha_u rightarrow infinity] or [eqn: alpha_{rm S} rightarrow infinity]. Consequently, we can see that the increase of radial uniformity of the temperature distribution directly implies the decrease of the von Mises stress distribution. That is, the radial uniformity of the temperature distribution can be considered an indicating measure of the thermal stress distribution. This point agrees well with the efforts reported elsewhere, where system configurations and processing conditions are designed to lower the radial temperature gradients in the crystal phase by using auxiliary heater, cold sink and thermal shield around the crystal [1, 2].

CONCLUSION

The optimization of the crystal surface temperature distribution is performed for the single crystal growth in the Czochralski process by considering the crystal quality and the growth rate. A systematic scheme is proposed based on the calculus of variations and the method of Lagrange multiplier. Then the Euler-Lagrange equations are derived and solved by the iterative numerical scheme proposed in this work. Through this work, we obtain the optimal distribution of the crystal surface temperature to provide an insight for the optimal design of thermal surroundings.

ACKNOWLEDGEMENT

This work has been partly supported by the grant from the KOSEF through ARC at POSTECH. This work has also been supported by the BK21 program of the Ministry of Education of Korea.

REFERENCES

1. D.E. Bornside, T.A. Kinney, R.A. Brown, J. Crystal Growth, 108 (1991) 779-805.

2. T. Tsukada, K. Kakinoki, M. Hozawa, N. Imaishi, K. Shimamura, T. Fukuda, J. Crystal Growth, 180 (1997) 543-550.

3. B.A. Boley, H.J. Weiner, □Theory of thermal stress□, Krieger, Malabar, FL (1985).

Fig. . Schematic diagram

Fig. . Convergence of the numerical scheme.

Fig. 3. Prototype of the optimal distributions of the state variables and the adjoint variables.

Fig. 4. Effects of [eqn: alpha_{rm S}] on [eqn: u``|_r=1].

Fig. 5. Effects of [eqn: alpha_u] on [eqn: u``|_r=1].