WO2002069391A1 - Method of forming a thin wax layer on a mounting block for mounting a semiconductor wafer - Google Patents

Abstract

A method of the invention is directed to forming thin layers of wax on surfaces of mounting blocks (10) for mounting semiconductor wafers on the mounting blocks so that each wax layer (18) has a desired wax layer shape and thickness selected to facilitate shaping of each wafer during subsequent polishing of the wafer. The method includes selecting values of at least some variables having an effect on the shape and thickness of the wax layer .The way layer is modelled numerically multiple times while changing at least one of the values to determine values for the desired wax layer shape and thickness. Wax layers are formed on the surfaces of the mounting blocks by using the determined values of the variables to achieve the desired wax layer shape and thickness for producing a consistent, desired shape of the wafer after polishing.

Description

METHOD OF FORMING A THIN WAX LAYER ON A MOUNTING BLOCK FOR MOUNTING A SEMICONDUCTOR WAFER

Background of the Invention

The present invention relates generally to depositing fluid layers on rotating disks, and more particularly to depositing wax layers on mounting blocks for mounting semiconductor wafers for single-side polishing.

A conventional process for polishing a single side of a semiconductor wafer includes "wax mounting" the wafer on a disk-shaped mounting block prior to polishing the wafer. A layer of wax is formed on a planar surface of the mounting block, and a back surface of the wafer is pressed against the wax to fix the wafer against the block.

Thereafter, the front surface of the wafer is polished by mounting the block in a single-side polishing apparatus and then forcing the block downward so that the front surface of the wafer is forced against a polishing pad mounted on the apparatus. The shape and thickness of the wax layer affects the shape of the polished wafer because the wax layer is interposed between the planar surface of the block and the wafer. Thus, any non- uniformities in the wax layer may affect the shape of the polished wafer.

Several parameters affect the shape and thickness of the wax layer. Among the parameters known to affect the shape and thickness are the speed at which the mounting block is spun, the length of time it is spun, and the temperature of the block. Conventionally, an operator or engineer might change one of the parameters to try to effect a change in the wax layer shape and thickness. However, the exact effect that the parameters have upon the shape and thickness is not known. Therefore, significant time is wasted in trial and error experimentation with the parameters to try to find the desired wax layer shape and thickness. Moreover, it is extremely difficult to discover the parameters which will yield exactly the desired shape and thickness simply by trial and error. One of the difficulties in discovering the parameters is that commonly available measurement tools in the industry cannot accurately determine the shape and thickness of the wax layer on the mounting block. Conventionally, one relies on the shape of the wafer after polishing to indirectly determine whether the shape and thickness of the wax layer was adequate. Such indirect determination is inaccurate since the shape of the wafer is also dependent on many other factors, such as the shape of the polishing pad and the magnitude of pressure applied during polishing.

With respect to the method of modelling the wax, several studies of fluid flows and spin coating are of general interest, though none of the studies have focused on wax layer modelling. Such studies are summarized as follows.

A. Purely Hydrodynamic Studies

Prandtl (1904) carried out the pioneering studies of fluid flows in which there is mass addition or mass removal at a boundary surface. The first theoretical description of the spin coating process was an early work by Emslie, Bonner and Peck (1958) where a Newtonian fluid was assumed to be rotating onto an infinite plane. In this study the final flow pattern was obtained assuming that a uniform thick film rigidly rotates with the disk, wherein the centrifugal forces per unit volume were balanced by the viscous shear forces caused by the gradient of radial velocities across the film. In Emslie et al., other influences, such as pressure distribution, gravitational forces, and Coriolis forces, were neglected. Emslie balanced the viscous and centrifugal forces per unit volume by using Stokes' equation and cylindrical coordinates with an origin in the center of rotation. Adding the continuity equation, the thickness as a function of elapsed time can be obtained where the film thickness H does not depend on the radial coordinate:

Where H₀, ω, v, t are respectively the initial film thickness, the disk angular velocity, the fluid kinematic viscosity and the time elapsed. Heat and mass transfers (e.g., evaporation of a solvent) into or out of the system are neglected.

Higgins (1986) analyzed the flow of a Newtonian liquid placed on an impulsively rotated disk by assuming a uniform thickness and by using the method of matched asymptotic expansions. Regh and Higgins (1988) analyzed the same flow using a finite difference method where the detailed behavior of the transient film flow and the inertia effect are expressed as a function of the Reynolds number.

B. Spin coating of Non-Newtonian fluids Acrivos, Shah, and Petersen (1960) extended the analysis of Emslie et al. (1958) to the flow of a non-Newtonian liquid that was characterized by the power law, showing that a uniformly thick liquid film cannot be obtained because of a nonlinear radial distribution of the viscous shear forces. Janecke and Shuldt (1984) used the equations obtained by Emslie et al. for the spin coating of non-Newtonian fluids. Janecke and Shuldt studied the coating flow of non-Newtonian fluid using both the power law and Carreau constitutive relationship observing that a thickness gradient is generated in the film. The Carreau liquid, which has been found to be suitable for polymer solutions, forms a uniformly thick film during centrifugation due to its Newtonian-like behavior at low shear rates while the Bingham and power-law liquids do not form uniformly thick films. Further results show that the film thickness and its uniformity are chiefly determined by the rheological properties of the fluids and not by the initial profile prior to spinning.

C. Evaporative models

For an accurate description of the entire spin coating process, the evaporation of the solvent and the increase in viscosity resulting from the increase of solute concentration must be taken into account. This was not considered in the purely hydrodynamic studies mentioned above.

Von Karman (1921) studied the problem of flow, mass and heat transfer on a rotating disk. He related the experimentally measured mass-transfer rates from a rotating disk in an infinite environment under laminar and turbulent conditions to the corresponding heat transfer process by means of an analogy method. Kreith, Taylor and Chong (1959) collected experimental data on the effect of an adiabatic surface placed at various distances parallel to the disk. Sparrow and Gregg (1959) extended the analysis to fluids of any Prandtl number.

Meyerhofer (1978) and Janecke (1983) proposed models incorporating the solvent evaporation and viscosity change using film-averaged equations. Flack et al. (1984) numerically simulated the process taking into account the non-Newtonian character of the solution and the radial distribution of film thickness. (The radial component of velocity was assumed to be inversely proportional to the radius, in contrast to all previous studies that found or assume a proportional relationship). Meyerhofer, Janecke and Flack also considered the changes in fluid viscosity and solvent diffusivity with changing polymer concentration. Results obtained from the model show that polymer film thickness is controlled by convective radial flow of the resist solution and solvent evaporation. The former process governs film thickness evolution during the early stages of the process, while the latter becomes significant in later stages.

Sukanek (1985) extended Meyerhofer's model (1978) and Bornside and Scriven (1987) performed an extensive review of the literature and described a mathematical model that included the solvent evaporation and the non-Newtonian behavior. However, many simplifications, such as using a constant evaporation rate, a constant solute concentration at the surface, and film-averaged equations, were also employed in these models. An accurate analysis of the atmospheric gas flow on the film cannot be adequately performed using such simplifications.

Ohara et al. (1989) intended to clarify in details the film thinning process by also considering the surrounding gas-phase to further improve the existing spin coating modeling processes. Regh and Higgins (1992) presented one of the most complete modeling works of the literature, reducing the equation system to a one-dimensionai problem. D. Experimental studies of Spin Coating

Many experimental studies have been presented in the literature reporting the resulting final film thickness, whereas Daughton and Givens (1979) reported a transient film thickness. Their studies were restricted to fluids and processes commonly associated with the semiconductor industry. High-speed photography was used to understand the dynamics of the spin coating process. Fluids with viscosity of 28, 43, 1 00, 2040 cP were used to determine the effects of fluid volume, initial film speed, final spinning speed, spinning acceleration and final spinning time on the resulting film thickness and thickness uniformity. These authors also presented an empirical expression, which describes the film thickness as a function of the measured parameters:

H + exp(-(αωf + β ))J (2.2)

Where H is the film thickness, μ is the liquid viscosity, ω is the angular velocity of the disk and ko, a, β, b must be empirically determined. Some values for fluids of different viscosity were reported as follows

μ (cP) β a β k₀

28 0.51 0.25 0J5 0.091

43 0.51 0.25 0.63 0.102

1100 0.81 0.110

2040 0.82 0.06 -0.71 0.110

Summary of the Invention

Among the several objects and features of the present invention may be noted the provision of a method which improves the shape of a polished wafer, the provision of such a method capable of more accurately predicting the shape and thickness of a wax layer deposited on a mounting block; and the provision of such a method which improves the yield of acceptable wafers. Briefly, a method of the invention is directed to forming thin layers of wax on surfaces of mounting blocks for mounting semiconductor wafers on the mounting blocks so that each wax layer has a desired wax layer shape and thickness selected to facilitate shaping of each wafer during subsequent polishing of the wafer. The method comprises selecting values of at least some variables having an effect on the shape and thickness of the wax layer. The wax layer is modelled numerically multiple times while changing at least one of the values to determine values for the desired wax layer shape and thickness. Wax layers are formed on the surfaces of the mounting blocks by using the determined values of the variables to achieve the desired wax layer shape and thickness for producing a consistent, desired shape of the wafer after polishing.

Another aspect of the invention is a method of predicting the shape and thickness of a wax layer formed on a mounting block for mounting a semiconductor wafer on the block comprising mounting the wafer on the mounting block and coating a surface of the wafer with a lubricating fluid. Values of at least some variables having an effect on the shape and thickness of the wax layer are selected, and wax is deposited on a surface of the wafer and the wafer is spun using the selected values to form a predictive wax layer on the wafer. The silicon wafer is demounted from the mounting block and the shape and thickness of the predictive wax layer is measured to predict the shape and thickness of the wax layer formed on the mounting block at the selected values. Other objects and features of the present invention will be in part apparent and in part pointed out hereinafter.

Brief Description of the Drawings

Fig. 1 A is a perspective view of a mounting block and wax deposited thereon; Fig. 1 B is an elevation of the block after wax is spread over the block;

Fig. 1 C is a schematic view of the block and a semiconductor wafer prior to mounting the wafer on the block;

Fig. 1 D is a schematic view of the wafer mounted on the block; Fig. 2A is graph showing an axial speed profile predicted by a model of the invention;

Fig. 2B is a graph showing a radial speed profile predicted by the model; Fig. 3A is a graph of a radial temperature profile predicted by the model; Fig. 3B is a graph of an axial temperature profile predicted by the model;

Figs. 4-7 are graphs of the thickness and shape of the wax layer predicted by the model using varying input data; and

Fig. 8 is a flowchart illustrating the model Fig. 9 is a flowchart illustrating an indirect measurement method. Corresponding reference characters indicate corresponding parts throughout the several views of the drawings.

Detailed Description of the Preferred Embodiment

Referring now to the drawings and in particular to Figs. 1 A-1 D, a disk- shaped mounting block, generally referred to as 10, has a back surface 11 and an extremely flat mounting surface 12. A semiconductor wafer W has a first mounting side 14 which faces the mounting block 10 when mounted thereon (Fig. 1D), and a second side 16 which will be polished after the wafer is mounted on the mounting block. The block 10 is made of ceramic material, though other materials are contemplated.

In a method of the invention, a thin wax layer 18 is formed on the mounting surface 12 of the mounting block 10 for mounting the semiconductor wafer W on the block such that the layer has a desired shape and thickness which is selected to facilitate shaping of the wafer during polishing of the wafer. The method comprises the steps of selecting values of at least some variables having an effect on the shape and thickness of the wax layer 18 and then modelling the wax layer numerically multiple times while changing at least one of the values to determine values of the variables for the desired wax layer shape. A wax layer 18 is then formed on the surface of the mounting block 10 by using the determined variables to achieve the desired wax layer shape and thickness for producing a consistent desired shape of the wafer W after polishing. The step of forming the wax layer 18 generally includes spin coating or wax deposition of the wax layer followed by heating of the block 10. More specifically, the block is mounted on a suitable device for spinning the block (not shown). The block 10 is cleaned by applying a cleaning solution, such as a peroxide solution, to the center of the block and spinning the block for a period of time. Referring to Fig. 1A, the spinning device spins the block 10 at an initial angular velocity (IAV) while liquid wax is applied, as by a nozzle (not shown), onto the surface of the block generally at the center of the block. Preferably, some of the cleaning solution remains on the block 10 when the liquid wax is deposited such that the block is wet when the wax is initially deposited. The temperature of the wax is controlled when the wax is applied to the block 10, and after a certain time, the block 10 and wax will reach an equilibrium temperature referred to as the block temperature which may also be controlled during the process. The block temperature may be measured using a pyrometer or by non-contact means such as an infrared video camera. The block 10 is spun for a period of time (referred to as spinning time) until the wax is spread over the surface 12 by centrifugal force (Fig. 1B shows the wax layer 18 after spinning). The angular velocity is increased during the spinning time at a fixed rate for an acceleration time period until the block 10 reaches a final angular velocity (FAV). The FAV is maintained for the remainder of the spinning time. In addition to spreading the wax, spinning causes evaporation of some of the solvent in the green wax. The block 10 is heated to evaporate solvents in the wax and thereby harden the wax. Preferably, the block 10 is heated by steam flow applied to the back surface 11 of the block. Meanwhile, the wafer W is cleaned to substantially eliminate particles on its surface. A vacuum is drawn around the block 10 and the wafer W to inhibit bubbles forming in the wax, and the block is pushed onto a mounting side of the wafer such that the hardened wax fixes the wafer on the mounting block. The block 10 and wafer W are then transferred to a polishing machine for polishing the second side of the wafer. Typically, the block 10 will be secured in a vacuum clamp of the polishing machine, and the block and wafer W will be forced against a polishing pad so that the second side of the wafer is polished. It is expected that any irregularity in the shape and thickness of the wax layer 18, i.e., any deviation from planarity, will be transferred to the wafer W during polishing and that such irregularity will appear in mirror image in the polished second side of the wafer. For example, if the wax layer 18 is concave, the second side is expected to be convex after polishing. However, other variables, such as polishing pressure and the polishing pad, may affect the shape of the second side of the wafer W. For example, if more pressure is applied by the polishing machine to the edges of the wafer W than the sides of the wafer during polishing, the second side of the wafer will be convex. Similarly, if the polishing pad is not planar, such non-uniformity will affect the planarity of the wafer W.

As stated above, values of at least some of the variables having an effect on the shape and thickness of the wax layer 18 are selected. The variables and the ranges of values selectable are listed in Table 1. The shape and thickness of the wax layer 18 is thereafter numerically modelled using the model described below. Generally, the model integrates the mass conservation equation in a two dimensional domain using the orthogonal collocation method in order to predict the shape and thickness of the wax layer 18.

TABLE 1

VARIABLE RANGE

Initial Angular Velocity (IAV) 0 / 1 ,000 rpm

Final Angular Velocity (FAV) 50 / 1 ,000 rpm Spinning Time 2 / 100 seconds

Acceleration Time 2 / 100 seconds

Ambient Air Temperature 18 / 40 ° C

Relative Humidity 0 / 100%

Wax Temperature up to 100 ° C Block Temperature up to 100 ° C Block Temperature Gradient up to 20 ° C

(from Block center to edge) Block Diameter up to 2 meters

Wax Layer Model

The model equations are the conservation equations of mass, energy and momentum. A cylindrical coordinate system is chosen because of the axial symmetry of the mounting block 10. The radial and axial velocity (v_r, v_z) are taken to be independent by the tangential (or azimuthal) coordinate. The nomenclature is provided in Table 2.

TABLE 2 (Nomenclature)

P wax density

9 gravity acceleration

V vectorial wax velocity t time

P pressure μ wax viscosity c_p specific heat of the wax

T wax temperature q heat conductive flux (heat transfer rate) k_r heat conductivity of the wax

C, concentration of the i-th component

Jι mass diffusive flux

0, diffusivity of the i-th component in the liquid mixture radial wax velocity

Vz axial wax velocity ₃ azimuthal (tangential) wax velocity r radial coordinate z axial coordinate s azimuthal (tangential) coordinate ω angular velocity k_dis_k thermal conductivity of the disk δ_disk half of the disk's thickness i_herm heat transfer coefficient λ_ev, latent heat of evaporation of the i-th component

T_b ^s temperature in the gas phase bulk k_c,ι mass transfer coefficient

Co, '^s mass fraction of the component i-th in the gas phase at the gas-liquid interface ω,^*'' mass fraction of the component i-th in the liquid phase at the gas-liquid interface ω,^ά,s mass fraction of the component i-th in the gas phase bulk

ω,^w mass fraction of the component i-th in the liquid phase bulk

P°(T) vapor pressure at temperature T MM_Tj average molecular weight of the liquid mixture

MM_T,g average molecular weight of the gas mixture

MM,- molecular weight of the component i-th

H thickness v_z ^s axial velocity at the free surface v_r ^s radial velocity at the free surface

Tas_k block temperature

T_E temperature at the edge of the block

Tc temperature in the center of the block

R__isk block radius τ constant temperature time

Re Reynolds number

Pr Prandtl number The equations of the model are as follows:

dv

*- = ω ²r + (1.1) dt p dz²

^"_ 1 a f δT a²τ ^{p C}»ιr⁼ k_τ (1.3) r 9r 9^r , ⁺ 8z^> -

The unknown variables are: the radial and axial velocity v_r, v_z, the temperature T of the wax, the mass concentration of the component "i" C, and the thickness H. The other variables (described below) are physical parameters and are known. The variables r and z are respectively the radial and the axial coordinate. The variable ω,^{* 9}, is the mass fraction of component / in the gas phase at the liquid-vapour interface and is calculated with the following thermodynamic equation:

where ω, is directly related to the mass concentration according to the equation:

The mass fraction in liquid phase at the interface is directly related to the mass concentration at the free surface (in liquid phase).

It is to be noted that the effects of convective heat transfer are assumed negligible and are ignored in Equation 1.3 and the effects of convective mass transfer are likewise assumed negligible and are ignored in Equation 1.4. Such assumptions are found to be accurate within a short time after introduction of the wax on the block 10. Since the equations are partial differential equations, they must be completed with the following boundary conditions: no slip conditions at the block surface: v_f = 0 at z = 0 (1.8)

v_z = 0 at z = 0 (1.9)

traction condition (equality of the stresses) at the free surface:

^dV-^r = Q at z = H (1.10) dz

null velocity in the center of the block (singular point): v_r = 0 at r = 0 (1.1 1 )

symmetry condition at the center of the block

— = 0 at r = 0 (1.12) ^dr thermal axial symmetry condition at the center of the block:

— = 0 at r = 0 (1.13) dr '

energy flow condition at the block surface: dT ₌ _ k_disk 3(T - T_disk ) _{a{ z = 0 (1 M)} dz k_r δ_disk

energy flow condition at the gas-liquid interface:

concentration axial-symmetry condition at the center of the block : dC,

= 0 at r = 0 (1.16) dr

null mass flow at the block surface:

^X- = 0 at z = 0 (1.17) dz

mass flow condition at the free surface

Because of the dependence of the unknown variables v_r, v₂, T and on time, the equations need an initial condition, which is defined as follows:

T=T°(r,z) at t=0

v_r=v_r°(r,z) at t=0

v_z=v_z°(r,z)=0 at t*0

Note that at time t=0 when wax is first introduced on the block 10, the boundary condition of Equation 1.8 is not completely accurate due to the moving boundary problem when the wax-block interface is growing. Since the moving boundary problem ceases very quickly, i.e., once the interface stops growing, there is little effect on the overall accuracy of the model.

Physical Parameters Estimation

Viscosity of the wax is dependent on temperature and water content and is suitably interpolated through the following exponential relationship.

//7rμ = 3.0268 + 6.122861 — - 1.6312954 ) + 38.007 [ ^ -0.0386424) (1.19)

where the water content w is expressed in mass fraction and the temperature is in °K.

Estimation of the activity coefficients γ used to determine the thermodynamic properties of the liquid vapor equilibrium is described with respect to the following regression equation:

\n(_{Y l}) = xl[A + 2(B -A)-_Xl] (1.20)

where A = 0.299753 and B = -13.224378.

The following physical and chemical properties are estimated using computational correlation. Gas diffusivity of water in air and of isopropyl alcohol in air are calculated using the relationship described by Chapman and Enskog [Prausnitz et al. (1988)]. The two diffusivities are:

Dwater, air = 2.13798 X 10^"5 m²/S Dc3H70H, air= 1.02641 x 10^'5 m²/s The viscosity and the thermal conductivity of the air are calculated using linear interpolation between 280 K and 300 K [Perry et al. (1984)]: μair- (0.175 + 5x10^" (Tamb-280)) 10^"4 Pa s k_air = (0.0247 + 8x10- (T_amb-280)) 10^"3 kW (m K)^"1 The thermal conductivity of the liquid is assumed to be equal to the thermal conductivity of pure water.

Heat and mass transfer coefficients estimation von Karman (1921) first analyzed the laminar flow near a rotating disk situated in a large body of quiescent fluid. Kreith et al. (1959) and Sparrow and Gregg (1959) studied the problems of mass and heat transfers from a rotating disk. Sparrow and Gregg found the following limiting relations for very low and very high Prandtl numbers (Pr):

= 0.62048 - PX* Pr→ ∞ (1.21 b)

where k_a,-_r is the thermal conductivity of the air and v_a/r is the kinematic viscosity of the air.

Note that the coefficient h is proportional to ω^α5 and that for laminar flow the transfer coefficients do not depend on the radial coordinate. The two asymptotic expressions have been combined to obtain the mass and heat transfer coefficients for Pr near to 1.

Adopting the Colburn analogy, the mass transfer problem is solved substituting the

Schmidt coefficient (Sc) for Pr.

Approximate calculations of rates of heat and mass transfer from a disk rotating in turbulent flow are available from Kreith (1959); the transition between laminar and turbulent flow occurs at about Re_c = 2.8 X 10⁵ where: Re = pω R disk

At a Reynolds greater than Re_c the flow is completely turbulent and the heat and mass transfer coefficients increase with increasing radius. In this example, the Reynolds number is calculated to be lower than about 5 X10⁴. Therefore the above laminar flow approximation can be used.

Integration of the equations with the orthogonal collocations method

The orthogonal collocations method is preferred to integrate the partial differential equations (PDE). The method uses the assumption that the solutions of the equations can be described through an expansion in polynomial form, in particular, the Lagrange polynomial is preferably used.

"Collocation points" (x_1f x₂, ..., x_n) are chosen in the domain of the mounting block. The Lagrange polynomial of degree n with a pole at xu is defined as:

An important characteristic of the Lagrange polynomial is that it is equal to 1 if it is calculated at x« and is equal to 0 if it is calculated in every xι with i≠k. The orthogonal collocations method is based on the assumption that the solutions of the equations written above are of the kind:

where a are unknown coefficients of the Lagrange series expansion to be found with the collocation method. The next step to transform the PDE into a solvable form is to assume that the supposed solution constituted by Equation 2.1 is exactly the real solution of the PDE in all the collocation points. In this manner and recalling the properties of Lagrange polynomials it is possible to reduce the set of partial differential equations to a set of Ordinary Differential Equations (ODE) and the boundary conditions to a set of Algebraic Differential Equations. In fact, all the partial derivatives with respect to the axial coordinate z and to the radial coordinate r are reduced to a sum of unknown variables.

For example, applying the polynomial approximation Equation 2.1 to the energy balance (1.3), the temperature dependence on rand z is:

where T_y are the unknown values of T at the collocation points. The partial derivative with respect to the time at the point fa, z_k) is: T(r_h,z_k) _ d

∑∑T, _;(z_A ) L,(r_Λ) dZ h (2.2). dt dt K.ι=ι /=1 dt

The partial derivative with respect to the radial coordinate is: i Ξl - T I (_T ^^L>fe) - τ SL,(r_h) (2.3) dr ι=1 =1 dr dr

Note that — '- M is computable for each polynomial L,.

Similarly, the second partial derivative with respect to the axial coordinate is:

The PDE (1.3) is reduced to a set of mxn ordinary differential equations (ODE) including mxn unknown T_tJ. This system assumes the form:

The mxn equation (2.5) is valid only in the internal points of the domain, while at the boundary points the boundary conditions have to be rewritten with the same substitutions Applying this technique to all the PDE's, a system of Algebraic Differential Equations is obtained. This system can be solved with an integrator such as SDASPK Routine, Petzold L.R., Brown P.N., Hindmarsh A.C., Ulrich C.W. (1990) available from http://www-users.cs.umn.edu/~petzold/.

Model of heater

In a simplified model of the heater, it is assumed that the shape of the layer 18 does not change during heating. However, the thickness decreases due to the evaporation of the solvent. The equations in the simplified model of the heating process are global mass balances at each collocation point on radial coordinates. The temperature of the block 10 during the heating process is assumed to be constant with respect to radial and axial coordinates because of the high thermal conductivity of the block 10 and the small thickness of the wax layer 18. The temperature of the block 10 increases during the heating process following a given temperature profile or ramp. The equations of the simplified model are as follows:

• Mass balance : dH dC.

where C, is the partial density of component i, H is the thickness and J,- is the mass flux due to the evaporation;

• Additivity of the volumes:

where d is the partial density of the specie i-th in the mixture while p,- is the density of the pure specie i. Modelling using computer processor and FORTRAN code

Preferably, the modeling step is performed using a computer processor, such as a personal computer. The model is programmed in a FORTRAN program (see appendix for listing) and graphically displayed using Excel software available from Microsoft Corporation, Redmond, Washington. Referring to Fig. 8, in step 101, an operator inputs selected values of the variables (operating conditions) plus the chemical and physical parameters of the wax (further described below) into the processor as input data. Collocation points are entered at step 102 and first derivatives of the Lagrange polynomial are evaluated at step 103 with respect to the collocation points. Using the first derivatives, the first derivative of the v, T, C field is evaluated with respect to time at the initial operating conditions and parameters in step 104. The program sets t equal to zero at step 105. The ordinary differential equation system generated in steps 103 and 104 are evaluated by the SDASPK routine at step 106. Still within the routine, the program determines a new time step Δt used to calculate new time t' at steps 107 and 108. Note that the SDASPK routine is shown in the appendix. At step 109, the program determines whether t is equal to t_f. If it is not, the SDASPK routine is repeated. When t is equal to t_f, the program outputs data at step 110. Preferably, the model outputs radial thickness and shape as a function of time, radial velocity vs. time and radial temperature vs. time. If the wax layer thickness and shape, as modelled, is not the desired wax layer thickness and shape, at least one of the input values is changed and the program is run repeatedly until a final model shows the desired thickness and shape. Thereafter, the wax layer 18 is formed on the block 10 as described above using the determined values of the variables according to the final model. Preferably, all the steps of the method are performed automatically. Examples

The program was run using the input data in Table 3.

TABLE 3

The program calculated the axial speed profile at points along the block radius at t=0.8s (Fig. 2A) and the radial speed profile at several axial positions (Fig. 2B). Radial speed shows a quadratic dependence with coordinate z once r is fixed. The temperature radial and axial profiles at several times are shown in Figs. 3A, 3B. Figure 3A shows that the axial temperature gradient is very small. Figure 3B shows that the initial wax temperature was greater than the block temperature and the air temperature. Figure 4 is a graph of the thickness of the wax layer 18 versus the radius of the block 10. Figure 4 shows that the shape of the wax layer 18 will be generally convex.

The program was run again using the same data except that final angular velocity (FAV) was changed from 320 to 600 rpm. Figure 5 is a graph of thickness versus radius at the increased FAV showing that the thickness is decreased but the shape is not significantly changed.

Another run was performed with the same data as the first run except that the block temperature was increased to 30°C. Increasing the temperature increases the evaporation rate so that the wax layer 18 is thinner, as shown in Fig. 6, and also causes the wax layer to be more convex.

The program was run a fourth time with the same data as the first run except that the thermal gradient in the block was inverted, i.e., the temperature at the center was 22.4°C and the temperature at the edge of the block was 26.4°C. As shown in Fig. 7, inversion of the gradient causes the shape of the wax layer 18 to change from convex to concave.

Sensitivity Analysis

The sensitivity of the wax layer thickness and roll-off (a measure of shape) was studied by defining sensitivity coefficients as follows:

dlnH

³thkJ lnx_j

_ d\n(roll off)

' roll-off, j ln^ where H is the wax layer thickness and x/ indicate the considered independent variable. Sensitivity coefficients (Table 4 were determined using the example data given above. TABLE 4

The coefficients provide a general idea of the effect that changes in certain variables should have on shape and thickness. The sensitivity coefficients are used to aid in choosing appropriate values of variables to more quickly arrive at the desired wax layer shape and thickness. As shown, wax layer thickness is more significantly influenced by FAV, and is somewhat influenced by spinning time and block temperature. Note that block temperature refers to the temperature of the block after it reaches equilibrium after dissipation of the thermal gradient. Wax roll off, a measure of shape, is significantly influenced by block temperature and thermal gradient and somewhat by acceleration time. Note that the initial profile of the wax layer 18 (the profile of the wax immediately after deposition on the block) does not influence the thickness and only slightly affects the roll off.

In a preferred embodiment of the invention, multiple program runs are performed using a different value of FAV for each run to determine the FAV for the desirable wax layer thickness and shape. Preferably, the program allows the user to vary the value of one variable in each run, and to quickly perform program runs so that the user may view up to 30 program runs per hour. Note that the other variables may or may not be held constant in such runs. Similarly, multiple program runs may be performed using a different spinning time for each run to determine the spinning time for the desirable wax layer. Also, multiple program runs may be performed using different block temperature and/or thermal gradient to determine the values of block temperature and thermal gradient for the desirable layer thickness and shape. Typically, one variable (spinning time, block temperature or thermal gradient) is varied while holding the other variables constant in successive runs. However, it is contemplated that the user may vary two or more variables simultaneously to determine the variables for the desirable layer. After determining the variables for the desirable wax layer, a wax layer 18 is formed on a mounting block 10 as described above such that the variables are controlled at the determined values during the process. Accordingly, the desired wax layer 18 is formed on the actual block 10 and the wafer W is mounted thereon. Thereafter, the wafer W is polished in a suitable single-side polishing machine. Preferably, the shape of the wafer W is measured and data from such measurement, such as the roll-off value for the wafer is communicated to the computer having the program thereon. Successive runs of the program are thereafter performed to optimize the wax layer shape based in part on the measurement of the actual wafer W. For example, if the measurement data shows that the wafer W is convex, the values of some of the variables are changed in the program to suitably change the predicted wax layer shape to correct the problem. Then, the wax layer 18 is formed on the block 10, a second wafer W is mounted and polished as described above, and the shape of the second wafer W is measured. The measurement data is again communicated to the computer, and successive runs of the program may be performed if the measurement data is not as desired.

Measurement procedure

An indirect measurement method is used to determine the thickness and shape of the wax layer 18 formed on the block 10 and to thereby confirm the model. A silicon wafer W is mounted on the mounting block 10 in the manner described above, i.e., forming a wax layer 18 on the block 10, heating the layer and pressing the wafer W onto the wax. Another wax layer 18 is formed on the wafer W, which is still mounted on the block 10, using the same selected variables as were used to form the wax layer on the mounting block. Thereafter, wafer W is demounted from the block 10 and the back side of the wafer (the side which faced the mounting block) is cleaned with isopropyl alcohol. The wafer W is transferred to a measuring machine capable of measuring the thickness of the wafer including the wax layer 18 thereon. A suitable machine is a Model ADE 3500 E3 measuring system, made by ADE Corporation of Westwood, Massachusetts 02090. The machine includes two capacitive probes spaced a known distance apart. Referring to Fig. 9, the wafer W is measured at step 201 by interposing the wafer between the probes, and each probe measures the distance between itself and the wafer surface or the wax layer 18 thereon at numerous points along the wafer. By comparing the data from the probes, point data is generated in polar coordinates at step 202, and the point data is transferred to a computer processor where it is transformed to cartesian coordinates (x, y, thk) at step 203. The wafer W is cleaned to remove the wax layer 18, and in step 204, the wafer is measured again and point data in polar coordinates is generated at step 205. The point data is transferred to a computer processor where it is transformed to cartesian coordinates (x, y, thk) at step 206. The processor is programmed with a suitable program for comparing the point data, such as a program in MATLAB code. At steps 207-208, the processor calculates the point-to-point difference between the two sets of data and generates a three-dimensional matrix in x, y and thk coordinates for the wax layer 18. At step 209, triangular interpolation is performed on the data using the MATLAB "griddata" function to uniformly re-distribute the three-dimensional matrix. At step 210, the a three- dimensional matrix is output wherein the thk values represent the thickness of wax at cartesian coordinates on the wafer W. Optionally, the processor outputs a chart of the matrix to graphically show the wax layer thickness at many points in a "topographical map" of the wafer surface at step 211. Thus, the program outputs the wax layer shape and thickness.

It is to be noted that the wax layer 18 used during wafer processing is formed on the mounting block 10, while the thickness measurement maps a wax layer formed on a wafer W. The underlying assumption is that the wax layer shape and thickness will be substantially identical regardless of whether the layer 18 is formed on a wafer W or a mounting block 10. The assumption is believed justified due to the fact that before the wax flows on the rotating disk (block or wafer), a thin film of water is present on the disk. The water "lubricates" the motion of the wax and the interaction between the wax and the material below are negligible.

Appendix A print-out of the program code for the model and SDASPK solver. AFPEKDΪX

C

C SPINCOM.F

C C it contains all the common statements for the program drygel.for

C

C all the vectors and matrices are automatically dimensioned by

C means of the following dimensioning parameters:

C

C NTPC number of collocation points

C

C ASSIGNED THE ABOVE PARAMETERS, ALL THE VECTORS AND MATRICES c ARE AUTOMATICALLY DIMENSIONED **********************************

C c

PARAMETER(NTPC=10) PARAMETER(NRAMP=100) C

C CHARACTER*70 TITLE

CHARACTER*70 DATE CHARACTER*9 PROBLEM CHARACTER* 1 IVIDEO CHARACTER* 10 RGUESS CHARACTER*9 SYSCOR

C

COMMON/ K2/IDST,NRT c

COMMON/CST/TITLE,DATE,PROBLEM,RVIDEO,RGUESS,SYSCOR

C

COMMON/CCA/NTX,NJX,N0X,NLX,NTY,NJY,N0Y,NLY COMMON/CCB/ROOTX(NTPC),AL 1 X(NTPC,NTPC),AL2X(NTPC,NTPC) C0MM0N/CCB/R00TY(NTPC),AL1Y(NTPC,NTPC),AL2Y(NTPC,NTPC) COMMON/CCD/ALX,BEX,ALY,BEY,WQX(NTPC,NTPC),WQY(NTPC,NTPC) COMMON/CCE/NTXML,NTXM2,NTYML ,NTYM2 COMMON/PR£C/ABSTOL,RELTOL c

COMMON NXN NNX,NΓNTER,TFΓN,NPRCON COMMON DATIPHY RHO_IOMEGAL,XMU0,XKC,XK,HTERM,CPLIQ,TAMB

COMMON/PHY2/ALPHA,RHOGAS,PMSEC,PMARIA,PMH20,YVAPI,BIOT,FMH20B COMMON PHY3 DIFH20,DELHEV,HTERML,UMΓREL C0MM0N PHY4 DIFIS0,XKCIS0,XKIS0,DELHEV2,PMIS0

C

COMMON/DIM/RMAX,RMAX 1 ,RCENTR,HMAX,HMAX1 NRRIFNZRIF,TAU,CRIF,TRIF COMMONNARA^fNTPC,NTPC)NZ(NTPC,NTPC),NRl(NTPC,NTPC) COMMON/VAR2/H(NTPC),Hl (NTPC)

COMMON/VAR3/CONC(NTPC,NTPC),CONC 1 (NTPC,NTPC),CISO(NTPC,NTPC) COMMON/VAR4/T(NTPC,NTPC),Tl (NTPC,NTPC),CIS01 (NTPCNTPC) COMMON VAR5/HFERMI,AFERMI,BFERMI,PPP02,EQCNTR,RCNTR

COMMON/ANGPAR/OMEGAI,OMEGAF,TACC,SPSAT,XKSAT,TSPIN COMMON/SAT/TSATI,TSATF,TMSC,TΓMETERM,T0(NTPC) C0MMON/USRJTLC0NV,FLWATER,FLIS0,CSI(NTPC),REPL,BIOTT,BIOTISO C 10 c = c -

C C

PROGRAM SPIN13bis 5

C BY F.RUBINI & G.L.VALENTE

C MODELLING OF SPIN COATING AND HEATING OPERATION

C (1999)

C 0 IMPLICIT REAL(A-H,0-Z)

EXTERNAL RES-JAC.PSOL

INCLUDE 'spincom.f

C

C sdasspk dimensioning according with full jacobian matrix and b C Newton method

C

PARAMETER(NNEQ=5*NTPC**2+NTPC)

PARAMETER(LIW=40+NNEQ)

PARAMETER(LRW=50+9*NNEQ+NNEQ**2) L COMMON/MIA/NPRTNT

C

DIMENSION RWORK(LRW),IWORK(LI )

DIMENSION INFO(20)

DIMENSION Y(NNEQ),YPRIME(NNEQ)₁DELTA(NNEQ) ^o DIMENSION RTOL(NNEQ),ATOL(NNEQ)

DIMENSION BASE( JNEQ),PARDIF(NNEQ,NNEQ) C FACCIO LTNPUT DEI DATI

OPEN(UNIT=14,FILE='VAROP.DAT')

OPEN(UNIT=1 ,FILE='DATI.DAT') 1 OPEN(UNIT=12,FILE='datiphy.dat') RHOGAS=1.16 !kg/m3

PMSEC=156.61 IKG/KMOL

PMARIA=29. ! KG7KMOL

PMH20=18. ! KG KMOL

PMISO=60. ! KG/KMOL

CALL RINPUT

CLOSE (12) CLOSE (14) CLOSE (15) CALL SINITO

NEQ=NNX IPAR=NEQ

C

C

C OUTPUT FILES

OPEN(UNIT=l 0,FILE=*RESULT.RIS') OPEN(UNIT=l l,FILE='SPESSORE.RIS') OPEN(UNIT=13,FrLE='TEMPER.RIS') OPEN(UNTT=15,FILE=^,CONCEN.RIS') OPEN(UNIT=40,FILE=MARANG.RIS')

C

C sdaspk inputs & initializations

C

OPEN(UNIT=20,FILE='dassl.dat^,) OPEN(UNIT=21,FILE='dassl2.daf)

CALL START(LNFO,IRES,RTOL,ATOL,RWOR ,LRW₎NNEQ,NEQ) CLOSE (20) CLOSE (21) C C CALLS COLLOC FOR THE ORTHOGONAL COLLOCATION METHOD CALL COLLOC C

C WRITES THE COLLOCATION POINTS VALUES ON OUTPUT FILES WRITE (11,*) (ROOTX(L)*RMAX*1000, L=1,NTX) WRITE (13,*) (ROOTX(L)*RMAX*1000, L=1,NTX)

WRITE (15,*) (ROOTX(L)*RMAX*1000, L=1,NTX) OPEN(UNIT=22,FILE='incmp.DAT')

C

C FIRST SET YPRIME(I)=0., THEN SET YPRIME c ACCORDING TO THE E.D. (NAV/ST) INSIDE "INIZIA"

C

DO 1=1, NNX YPRIME(I)=0. ENDDO CALL INLZIA

CALL DCBA(NEQ,Y,YPPJME)

CL0SE(22)

TX=0. DTIME=TFΓN/NΓNTER c

C TIME INTEGRATION: LOOP SDASPK

C c NPRINT=1

DO 1 K=1,NTNTER TXOUT=DTΓME*K 100 CALL SDASPK(RES,NEQ,TX,Y,YPRME,TXOUT,INFO,RTOL,ATOL,ROID, + RWORK,LRW,RWORK,LIW,RPAR,LPAR,JAC,PSOL) IFILE=10

IF(IDID.LT.0) THEN

CALL CRISIS(NEQ,IDrD,INFO,Y) ^• GOTO 100 ENDLF C printing of the partial results

IF (K.EQ.(NPRTNT*NPRCON)) THEN WRITE(*,*) TXOUT*TAU CALL OUTPAR(TX) ENDΓF ^• 1 CONTINUE CALL OKAY

CLOSE(IO) CLOSE(ll) CLOSE(13) CLOSE(15) CLOSE(40) END C

C SUBROUTINE RES(TX,NYPRIME,DELTA,IRES,RPAR,ΪPAR) C

C the subroutine defines the residuals of the system of

C algebraic-differential equations DELTA(i)=G(T,Y/YPRIME)

C NEQ=IPAR

C

IMPLICIT REAL(A-H,0-Z)

INCLUDE 'spincom.f

DIMENSION Y(IPAR),YPRIME(IPAR),DELTA(IPAR) c NEQ=IPAR C

CALL ABCD(NEQ,Y,YPRIME) C C NDELMAX IS THE CURRENT DELTA VALUE

NDELMAX=0

C BOUNDARY CONDITIONS ON Z=0 VR(0)=0, VZ(0)=0

DO 10 I=l,NTX

IEQ=I DELTA(IEQ)=VR(I,1)

FRANCO=DELTA(I)

IEQ=I+NTX

DELTA(IEQ)=VZ(I,1) 10 CONTINUE NDELMAX=2*NTX

DO 11 J=2,NTYM1

IEQ=(J-1)+NDELMAX

DELTA(IEQ)=VR(1,J)

D1VZR=0. DO N=l,NTX

D1VZR=D1VZR+AL1X(1,N)*VZ(N,J)

ENDDO

IEQ=(J- l)+NDELMAX+(NTY-2)

DELTA(IEQ)=D1VZR 11 CONTINUE

NDELMAX=NDELMAX+2*(NTY-2)

C THE BOUNDARY CONDITION (B.C.) IN R=0. IS D1HR=0.

D1HR=0. DO N=l,NTX D1HR=D1HR+AL1X(1,N)*H(N)

ENDDO

IEQ=NDELMAX+1 DELTA(IEQ)=D1HR NDELM AX=NDELM AX+ 1 C ON R=0 AND Z=H THE B.C. FOR VZ IS Dl VZR=0

D1VZR=0.

DO N=l.NTX

D 1 VZR=D 1 VZR+ AL 1 X( 1 ,N)* VZ(N,NTY)

ENDDO IEQ=NDELMAX÷1

DELTA(IEQ)=D1VZR

NDELMAX=NDELMAX+ 1 C B.C. ON Z=H IS Dl VRZ(H)=0

DO 12 1=1, NTX D1VRZ=0. C I IS THE POINT (ON r), M IS THE POLYNOMIAL D0 14 M=1,NTY

D1VRZ=D1VRZ+AL1Y(NTY,M)*VR(I,M) 14 CONTINUE

IEQ=I+NDELMAX DELTA(IEQ)=D1VRZ 12 CONTINUE

NDELMAX=NDELMAX+NTX

C B.C. D1CZ=0. Z=0.

IF (EQCNTR.GE.2) THEN IONLY IF EQCNTR>=2 USE EVAPORATION

DO 1=1, NTX D1CZ=0.

D1CISOZ=0. DO M=l,NTY

D1CZ=D1CZ+AL1Y(1,M)*C0NC(I,M) D1CIS0Z=D1CIS0Z+AL1Y(1,M)*CIS0(I,M) ENDDO

ΓEQ=NDELMAX+I DELTA(IEQ)=D1CZ IEQ=NDELMAX+NTX+I DELTA(IEQ)=D 1CISOZ ENDDO

NDELMAX=NDELMAX+2*NTX

C B.C. D1CR=0. IN R=0.

DO J=2,NTYM1

D1CR=0. DICTSORO.

DO N=l,NTX

D1CR=D1CR+AL1X(1,N)*C0NC(N,J

DlCISOR=DlCISOR+ALlX(l,N)*CISO(N,J)

ENDDO IEQ=NDELMAX+J-1

DELTA(IEQ)=D1CR

IEQ=NDELMAX+NTYM2+J-1

DELTA(IEQ)=DlCISOR

ENDDO NDELMAX=NDELMAX÷-2*NTYM2

C B.C. D1CZ=-BI0T*RH0GAS/CRIF*(FMH20G-FMH20B) IN Z=l

DO 1=1 , NTX D1CZ=0. D1CISOZ=0. DO M=1,NTY

D 1 CZ=D 1CZ+AL1 Y(NTY,M)*CONC(I,M) 5 D1CIS0Z=D1CIS0Z+AL1Y(NTY,M)*CIS0(I,M)

ENDDO C

FMH20L=CONC(LNTY)*CRTF/RHO FMISOL-CISO(I,NTY)*CRTF/RHO

1 XLIQI=FMH20L/PMH20/(FMH20L/PMH20+FMISOL/PMISO+(l -FMH20L-FMISOL)

1 /PMSEC)

XISOI=FMISOL/PMISO/(FMH20L/PMH20+FMISOL/PMISO+(l-FMH20L-FMISOL) 1 /PMSEC)

TI=T(I,NTY)*TRTF 15 YVAPI=P0_H2O(TI)/l .*XLIQI*GAMMA(XLIQI)

YISOI=P0_ISO(TI)/l .*XISOI

FMH20G=YVAPI*PMH20/(YVAPI*PMH20+YISOI*PMISO+(l-YVAPI-YISOI) 1 *PMARIA)

FMISOG=YISOI*PMISO/(YNAPI*PMH20+YISOI*PMISO+(l-YVAPI-YISOI) 20 1 *PMARIA)

C

ΓEQ=ΝDELMAX+I

BIOT=XKC/DLFH20*HMAX*H(I)

DELTA(IEQ)=-D 1 CZ-BI0T*RH0GAS/CRIF*(FMH20G-FMH20B)

25 C

IEQ=NDELMAX+NTX+I

BIOTISO=XKCISO/DIFISO*HMAX*H(I)

DELTA(IEQ)=-D1CISOZ-BIOTISO*RHOGAS/CRTF*(FMISOG-0.)

ENDDO 30 NDELMAX=NDELMAX+2*NTX

ENDIF !END OF THE EQCNTR CHECK

C B.C. D1TZ=-HTERM*H(I)/XK*(T(1,J)-TSAT) IN Z=0.

IF (EQCNTR.EQ.3) THEN 1=1 35 DO 1=1, NTX

IF (TX.LT.TIMETERM) THEN TSAT=(TSATI-T0(I))*TX/TIMETERM+T0(I) ELSE

TSAT=TSATI « G ENDIF

DTS AT=3/SPS AT*(T(1, 1 )-TS AT) Q=XKSAT*DTSAT

D1TZ=0.

D0 M=1,NTY

D1TZ=D1TZ+AL1Y(1,M)*T(LM) ENDDO

IEQ=NDELMAX+I

DELTA(IEQ)=-D 1TZ-Q*HMAX*H(I)/XK

ENDDO

NDELMAX=NDELMAX+NTX c B.C. DITR=O. ΓNR=O.

DO J=2,NTYM1

D1TR=0.

DO N=l,NTX

D 1 TR=D 1 TR+AL 1 X( 1 ,N)*T(N, J) ENDDO

IEQ=NDELMAX+J-1

DELTA(IEQ)=D1TR ENDDO NDELMAX=NDELMAX+NTYM2 C B.C. D1TZ=-BIOTT*(T(I,NTY)-TAMB)+ADΓM1*D1CZ+ADIM2*D1CISOZ

IN Z=l

DO 1=1, NTX D1TZ=0. D1CZ=0. D1CISOZ=0.

DO M=l,NTY

D1TZ=D1TZ+AL1Y(NTY,M)*T(I,M) D 1 CZ=D 1 CZ+AL l Y(NTY,M)*CONC(I,M) D1CIS0Z=D1CIS0Z+AL1Y(NTY,M)*CIS0(I,M) ENDDO

BIOTT=HTERM*H(I)*HMAX XK ADIMI=DΓFH20*CRΓF*DELHEV/(XK*TRΓF) ADIM2=DRFIS 0* CRIF*DELHEV2/(XK*TRJF) IEQ=NDELMAX+I DELTA(IEQ)=-D 1TZ-BI0TT*(T(I,NTY)-TAMB)+ADIM 1 *D 1 CZ+ADIM2*D 1 CISOZ

ENDDO

ND ELMAX=ND ELMAX+NTX ENDIF ! END OF THE EQCNTR CHECK

C C EQUATIONS (VALID FOR POINTS INSIDE THE GRID)

C c

C

C CONTINUITY EQUATIONS

DO 16 I=2,NTX R=ROOTX(I)

D0 18 J=2,NTY C FOR EACH i j DEFINE dVR dr AND dVZ/dz

DIVRR=0.

DO N=1,NTX D1VRR=D1VRR+AL1X(I,N)*VR(N,J)

ENDDO

D1VZZ=0.

D0 M=1,NTY

D1VZZ=D1VZZ+AL1Y(J,M)*VZ(LM) ENDDO

IEQ=NDELMAX+(J-2)*NTXM1+I-1

C REMOVE DIMENSIONS OF THE EQUATIONS

VALRUB=VZRΓF*RMAX/(HMAX*VRRLF)

DELTA(IEQ)=(D1VRR+VR(I,J) R)+VALRUB/H(I)*D1VZZ 18 CONTINUE 16 CONTINUE

NDELMAX=NDELMAX+NTYM1 *NTXM 1

C

C — NAVIER-STOKES EQUATION IN R DIRECTION (ONLY INTERNAL POINTS)

DO 20 I=2,NTX

R=ROOTX(I)

DO 22 J=2,NTYM1 C FOR EACH i,j DEFINE d2VR/dz2 D2VRZ=0.

D1VRZ=0.

D1VZZ=0.

DO M=1,NTY

D2VRZ=D2VRZ+AL2Y(J,M)*VR(I,M) D1VRZ=D1VRZ+AL1Y(J,M)*VR(I₎M)

D 1 VZZ=D 1 VZZ+AL 1 Y(J,M)* VZ(LM)

ENDDO

D1VRR=0.

DO N=l,NTX D1VRR=D1VRR+AL1X(I,N)*VR(N,J)

ENDDO

TEMPE=T(I,J)*TRΓF CONCEN=CONC(I,J)*CRIF RE=RHO*HMAX*H(I)*VRRIF/XMU(TEMPE,CONCEN) IEQ=NDELMAX+(J-2)*(NTX-1)+(I-1)

DELTA(IEQ)=-VR1(I,J)+(0MEGA(TX)**2*R)*(TMAXATIRIF)**2+1/ 1 (RE*(HMAX/RMAX)*H(I) )*D2VRZ 3 22 CONTINUE 20 CONTINUE

C

C H EQUATION C

NDELMAX=NDELMAX+NTXM1 *NTYM2 DO 25 I=2,NTX

D1HR=0.

DO N=l,NTX 15 D1HR=D1HR+AL1X(L1ST)*H(N)

ENDDO

FMH20L=CONC(I,NTY)*CRIF/RHO

FMISOL=CISO(I,NTY)*CRTF/RHO

XLIQI=FMH20L/PMH20/(FMH20L/PMH20+FMISOL/PMISO+(l-FMH20L-FMISOL) 20 1 PMSEC)

XISOI=FMISOL/PMIS0/(FMH2OL/PMH2O+FMISOL/PMIS0+(l-FMH2OL-FMISOL) 1 PMSEC)

TI=T(I,NTY)*TRLF

YVAPI=P0_H2O(TI)/l.*XLIQI*GAMMA(XLIQI) ²⁵ YISOI=P0_ISO(TI)/l.*XISOI

FMH20G=YVAPI*PMH20/(YVAPI*PMH20+YISOI*PMISO+(l-YVAPI-YISOI) 1 *PMARIA)

FMISOG=YISOI*PMISO/(YVAPI*PMH20+YISOI*PMISO+(l-YVAPI-YISOI) 1 *PMARIA) 30 C

IEQ=NDELMAX+(I-1) C GIANNI NUMBER

GIANNI=RMAX*VZRIF/(VRRIF*HMAX) IF (EQCNTR.EQ.1) THEN 35 DELTA(IEQ)=H1(I)-GIANNI*VZ(I,NTY)+D1HR*VR(I,NTY)

ELSE

DELTA(IEQ)=Hl(I)-GIANNI^:i:VZ(I,NTY) -DlHR*VR(I,NTY)+TAU/HMAX*RHOGAS* 1 (XKC/1000.*(FMH2OG-FMH2OB)+XKCISO/789.*(FMISOG-0.)) _{ή G} ENDIF

25 CONTINUE

NDELMAX=NDELMAX÷NTXM1 c

C EQUATION FOR C : FICK C

IF (EQCNTR.GE.2) THEN DO 27 I=2,NTX

R=ROOTX(I) DO 29 J=2,NTYM1 D2CZ=0. D2CISOZ=0. 0 DO M=l,NTY

D2CZ=D2CZ+AL2Y(J,M)*CONC(I,M) D2CISOZ=D2CISOZ+AL2Y(J,M)*CISO(I,M) ENDDO D1CR=0. 5 D1CISOR=0.

D2CR=0. D2CISOR=0. DO N=l,NTX

D 1 CR=D 1 CR+AL 1X(LN)*C0NC(N, J) 0 D1CIS0R=D1CIS0R+AL1X(I,N)*CIS0(N,J)

D2CR=D2CR+AL2X(I,N)*CONC(N,J)

D2CISOR=D2CISOR+AL2X(I,N)*CISO(N,J) ENDDO

IEQ=NDELMAX+(J-2)*(NTX-1)+(I-1) 5 DELTA(ffiQ)=-CONCl(I,J)+TAU*DIFH20*(l RMAX**2*(DlCR/R+D2CR)+l/

1 (H(I)**2*HMAX**2)*D2CZ)

TEQ=NDELMAX+NTXM1 *NTYM2+(J-2)*(NTX- 1 )+(I- 1 )

DELTA(IEQ)=-CIS01(I,J)+TAU*DLFISO*(1 RMAX**2*(D1CISOR/R+D2CISOR) 1 +l/(H(I)**2*HMAX**2)*D2CrSOZ) _C 29 CONTINUE

27 CONTINUE

NDELMAX=NDELMAX+2*NTXM 1 *(NTYM2)

ENDIF !END OF THE CONTROL OF EQCNTR

5 C

C EQUATION FOR T : FOURIER EQUATION C

IF (EQCNTR.EQ.3) THEN DO 37 I=2,NTX _fl R=ROOTX(I)

DO 39 J=2,NTYM1 D2TZ=0. DO M=1,NTY

D2TZ=D2TZ+AL2Y(J,M)*T(I.M) ENDDO D1TR=0. D2TR=0. DO N=l,NTX

D 1 TR=D 1 TR+AL 1 X(I,N) *T(N, J) D2TR=D2TR+AL2X(I,N)*T(N,J)

ENDDO

IEQ=NDELMAX+(J-2)*(NTX-1)+(I-1) ALPHA=XK/(RHO*CPLIQ)

DELTA(IEQ)=-T1(I,J)+TAU*ALPHA*(1/TMAX**2*(D1TP R+D2TR)+ 1 1/(H(I)**2*HMAX**2)*D2TZ)

39 CONTINUE

37 CONTINUE

NDELMAX=NDELMAX+NTXM 1 *NT YM2

ENDIF !END OF THE CONTROL OF EQCNTR

FLCONV=0.

CMED=0.

TMEDIA=0.

DO M=l,NTY

Z=ROOTY(M)

FLCONV=FLCONV+WQY(NTX,M)*VR(NTX,M)*VRRIF*HMAX*H(NTX) CMED=CMED+WQY(NTX,M)*CONC(NTX,M)*CRTF TMEDIA=TMEDIA+WQY(NTX,M)*T(NTX,M)*TRΓF ENDDO REPL=RHO*FLCONV/XMU(TMEDIA,CMED) RETURN

END

SUBROUTINE DCBA(NEQ,Y,YPRTME) c

C CHANGE FROM THE CURRENT VARIABLES (V,H,-.) TO THE Y FOR SDASPK C this routine does the opposite job of the routine ABCD

IMPLICIT REAL(A-H,0-Z) DIMENSION Y(NEQ),YPRIME(NEQ)

INCLUDE 'spincom.f c

DO 1 JX=1,NTX D0 2 JY=1,NTY JS=(JX-1)*NTY+JY

Y(JS)=VR(JX,JY) 2 CONTINUE

1 CONTINUE

NN1=NTX*NTY D0 3 JX=1,NTX

JS=NN1+JX Y(JS)=H(JX) CONTINUE NN2=NN1+NTX D0 4 JX=1,NTX D0 5 JY=1,NTY JS=NN2+(JX-1)*NTY+JY Y(JS)=VZ(JX,JY)

10 CONTINUE

4 CONTINUE

NN3=NN2+NTX*NTY

IF (EQCNTR.GE.2) THEN

DO 10 JX=l,NTX

DO 11 JY=1,NTY

15

JS=NN3+( JX- 1 ) *NTY÷JY

Y(JS)=CONC(JX,JY)

11 CONTINUE

10 CONTINUE

NN4=NN3+NTX*NTY

20

DO 20 JX=l,NTX

DO 21 TY=1,NTY

JS=NN4+(JX-1)*NTY+JY

Y(JS)=CISO(JX,JY)

21 CONTINUE ^έb 20 CONTINUE

NN5=NN4+NTX*NTY

ENDIF

IF (EQCNTR.EQ.3) THEN

DO 50 JX=1,NTX

50 DO 51 JY=1,NTY

JS=NN5+(JX-1)*NTY+JY

Y(JS)=T(JX,JY)

51 CONTINUE

50 CONTINUE

^3b ENDIF

C TIME DERIVATIVES

D0 6 IX=1,NTX

D0 7 JY=1,NTY

JS=(JX-1)*NTY+JY

40 YPRIME(JS)=VR1(JX,JY

7 CONTINUE

6 CONTINUE

D0 8 JX=l,NTX

JS=NN1+JX YPRIME(JS)=H1(JX) 8 CONTINUE

IF (EQCNTR.GE.2) THEN DO 40 JX=l,NTX D0 41 JY=1,NTY

JS=NN3+(JX-1)*NTY+JY YPRLME(JS)=CONC 1 (JX, JY) 41 CONTINUE

40 CONTINUE D0 45 JX=1_;NTX

D0 46 JY=1,NTY IS=NN4+(JX-1)*NTY+JY YPRLME(JS)=CIS01(JX,JY) 46 CONTINUE 45 CONTINUE

ENDIF

IF (EQCNTR.EQ.3) THEN DO 60 JX=1,NTX DO 61 JY=1,NTY JS=NN5+(JX-1)*NTY+JY YPRTME(JS)=T1(JX,JY) 61 CONTINUE 60 CONTINUE ENDIF RETURN END C c

SUBROUTINE ABCD(NEQ,Y,YPRIME) C CHANGE FROM Y (SDASPK) TO CURRENT VARIABLES C c

DIMENSION Y(NEQ),YPRIME(NEQ)

INCLUDE 'SPΓNCOM.F

DO 1 JX=1,NTX DO 2 JY=l,NTY JS=(JX-1)*NTY+JY VR(JX,JY)=Y(JS) 2 CONTINUE

1 CONTINUE

NN1=NTX*NTY

DO 3 JX=1,NTX

JS=NN1÷JX H(JX)=Y(JS) CONTINUE NN2=NN1+NTX D0 4 JX=1, TX D0 5 JY=1,NTY

JS=NN2+(LX- 1 )*NTY+JY VZ(IX,JY)=Y(JS) CONTINUE CONTINUE NN3=NN2+NTX*NTY

IF (EQCNTR.GE.2) THEN DO 10 JX=l,NTX DO 11 JY=1,NTY JS=NN3+(JX-1)*NTY+TY CONC(JX,TY)=Y(IS) CONTINUE CONTINUE

NN4=NN3+NTX*NTY D0 15 JX=1,NTX DO 16 JY=1,NTY

JS=NN4+(JX-1)*NTY+JY CISO(JX,JY)=Y(JS) CONTINUE CONTINUE NN5=NN4+NTX*NTY

ENDIF

IF (EQCNTR.EQ.3) THEN DO 30 JX=1,NTX DO 31 JY=1,NTY JS=NN5+(JX-1)*NTY+JY

T(JX,JY)=Y(JS) CONTINUE CONTINUE ENDIF TIME DERIVATIVES

D0 6 JX=1,NTX DO 7 JY=l,NTY JS=(JX-1)*NTY+JY VR1(JX,JY)=YPRTME(JS) CONTINUE CONTINUE DO 8 JX=1,NTX

JS=NN1+JX H1(JX)=YPRIME(JS) 8 CONTINUE

IF (EQCNTR.GE.2) THEN DO 20 JX=1,NTX DO 25 JY=1,NTY JS=NN3+(JX-1)*NTY+JY

CONC 1 (JX, JY)=YPRIME(JS) 25 CONTINUE

20 CONTINUE

DO 45 JX=1,NTX DO 46 JY=1,NTY

JS=NN4+(JX-1)*NTY+JY CIS01(JX,JY)=YPRIME(JS) 46 CONTINUE

45 CONTINUE ENDIF

TF (EQCNTR.EQ.3) THEN DO 70 JX=l,NTX D0 75 JY=1,NTY JS=NN5+(JX-1)*NTY+JY T1(JX,JY)=YPRIME(JS)

75 CONTINUE

70 CONTINUE ENDLF RETURN END

SUBROUTINE INTZIA C

C REMOVING DIMENSIONS TO THE VARIABLES C CURRENT VARIABLES INITIALIZATION

C

INCLUDE 'SPINCOM.F

NEQ=NNX IFILE=22 DO 1 JY=1,NTY

READ(IFILE,*)(VR(JX,JY),JX=1,NTX) 1 CONTINUE DOK=l,NTX DO L=l,NTY VR(K,L)=VR(K,L)/VRRIF pippo=vr(k,l) ENDDO ENDDO DO 2 JY=1,NTY

READ(IFILE,*)(VZ(JX,JY),JX=1 ,NTX) CONTINUE

DOK=l,NTX DOL=l,NTY

VZ(K,L)=VZ(K,L)VZRIF ENDDO ENDDO

DO200JY=l,NTY READ(IFILE,*)(CONC(JX,JY), JX=1 ,NTX)

FRANCO=CONC( 1 , JY) CONTINUE

DO250JY=l,NTY

READ(IFILE,*)(CISO(JX,JY),TX=l,NTX) CONTINUE DOK=l,NTX DOL=l,NTY

CONC(K,L)=CONC(K,L)/CRTF CISO(K,L)=CISO(K,L)/CRIF ENDDO

ENDDO DO 300 JY=1,NTY

READ(ΓFΓLE,*)(T(JX,JY),JX=I,NTX) CONTINUE DOK- NTX

DOL=l,NTY T(K,L)=T(K,L)/TRIF ENDDO ENDDO TO INITIALIZATION FOR THE B.C. ON THE BLOCK

DOK=l,NTX

T0(K)=T(K,1)

ENDDO

IF(PIPP02.EQ.1)THEN READ(IFILE,*) (H(JX),JX=1,NTX)

DO =l,NTX

H(K)=H(K)/HMAX

ENDDO

ENDIF IF (PIPPO2.EQ.0) THEN

DO M=1,NTX

R=ROOTX(M)

H(M)=FERMI(R)/HMAX

WRITE(*,*) H(M)*HMAX ENDDO ENDIF C

C D2VRZ CALCULATION DO 8 I=2,NTX

R=ROOTX(I)

DO 9 J=2,NTYM1

D2VRZ=0.

D1VRZ=0. D1VZZ=0.

DO M=l,NTY

D2VRZ=D2VRZ+AL2Y(J,M)*VR(I,M)

D 1 VRZ-D 1 VRZ+AL1 Y(J,M)* VR(I,M)

D1VZZ=D1VZZ+AL1Y(J,M)*VZ(I,M) ENDDO

D1VRR=0.

DO N=l,NTX

D1VRR=D1VRR+AL1X(I,N)*VR(N,J)

ENDDO TEMPE=T(LJ)*TRTF

CONCEN=CONC(I,J)*CRTF

RE=RHO*HMAX*VRRIF/XMU(TEMPE,CONCEN)*(r3MAXΛRMAX)

VR1(I ) 0MEGA(TX)**2*R)*(RMAX/VRRIF)**2+1/(RE*H(I)**2)*D2VRZ

TEST=VR1(I,J) 9 CONTINUE

D1HR=0.

DO N=l,NTX

D1HR=D1HR+AL1X(I,N)*H(N)

ENDDO FMH20L=CONC(LNTY)*CRLF/RHO

FMISOL=CISO(I,NTY)*CRIF/RHO

XLIQI=FMH20L PMH2O/(FMH2OL PMH2O+FMISOL PMISO+(l-FMH20L-FMISOL) 1 /PMSEC)

XISOI=FMISOL/PMISO/(FMH20L/PMH20+FMISOL/PMISO+(1-FMH20L-FMISOL) 1 /PMSEC)

TI=T(I,NTY)*TRIF

YVAPI=P0_H2O(TI)/l.*XLIQI*GAMMA(XLIQI) YISOI=P0_ISO(TI)/l .*XISOI

FMH20G=YVAPI*PMH20/(YVAPI*PMH20+YISOI*PMISO+(l-YVAPI-YISOI) 1 *PMARIA)

FMISOG=YISOI*PMISO/(YVAPI*PMH20+YISOI*PMISO+(1-YVAPI-YISOI) 1 *PxMARIA)

GIA ΠSΓI=RMAX*VZRΓF/(VRRIF^"I'HMAX) IF (EQCNTR.EQ.1) THEN

H1(I)=GIANNI*VZ(I,NTY)-D1HR*VR(I,NTY)

ELSE

H1(I)=GIANNI*VZ(I,NTY)-D1HR*VR(I,NTY)-TAU/HMAX*RH0GAS*( 1 XKC/1000.*(FMH2OG-FMH2OB)+XKCTSO/J89.*(FMISOG-0.))

ENDIF

TEST2=H1(I) 8 CONTINUE

C D0 27 I=2,NTX

DO 29 J=2,NTYM1

D2CZ=0.

D2CISOZ=0.

DO M=l,NTY D2CZ=D2CZ+AL2Y(J,M)*CONC(I,M)

D2CISOZ=D2CISOZ+AL2Y(J,M)*CISO(I,M)

ENDDO

D1CR=0.

D2CR=0.

D1CISOR=0.

D2CISOR ⁾.

DO N=1,NTX

DlCR=DlCR+ALlX(I,N)*CONC(N,J)

D2CR=D2CR+AL2X(I,N)*CONC(N,J) D1CIS0R=D1CIS0R+AL1X(I,N)*CIS0(N,J)

D2CISOR=D2CISOR+AL2X(I,N)*CISO(N,J)

ENDDO

IEQ=NDELMAX+(J-2)*(NTX-1)+(I-1)

C0NC1(I,J)=TAU*DIFH20*(1/RMAX**2*(D1CR R+D2CR)+1/ 1 (H(I)**2*HMAX**2)*D2CZ)

CIS01(I,J)=TAU*DIFISO*(l RMAX**2*(DlCISORyR+D2CISOR)+l/ 1 (H(I)**2*HMAX**2)*D2CISOZ) 29 CONTINUE 27 CONTINUE

C

DO 37 I=2,NTX D0 39 J=2,NTYM1 D2TZ=0. DO M=l,NTY

D2TZ=D2TZ+AL2Y(J,M)*T(I,M) ENDDO D1TR=0. D2TR=0. DO N=l,NTX

D 1 TR=D 1 TR+ AL 1 X(I,N)*T(N, J) D2TR=D2TR+AL2X(I,N)*T(N,I) ENDDO

IEQ=NDELMAX+(J-2)*(NTX-1 )+(I- 1 ) T1(I,J)=TAU*ALPHA*(1/RMAX**2*(D1TR/R+D2TR)+ 1 1/(H(I)**2*HMAX**2)*D2TZ)

39 CONTINUE

37 CONTINUE

RETURN END

C

SUBROUTINE RTNPUT

C

C DATA IMPORT FROM FILE

IMPLICIT REAL(A-H,0-Z) INCLUDE 'spincom.f

READ(12,*)RH0 READ(12,*)HMAX1

HMAX=HMAX1 *2.54/l 00 ! IN M

READ( 12,*) VRRIF READ(12,*)VZRIF C TRANSPORT PROPERTIES

READ(12,*)DIFH20

READ(12,*)DΓFISO

READ(12,*)CRΓF

READ(12,*)XKSAT READ(12,*)SPSAT

READ(12,*)XK C

C TERMODYNAMICAL PROPERTIES

READ(12,*)CPLIQ READ(12,*)DELHEV

R£AD(12,*)DELHEV2 READ(12,*)TPJF READ( 12,*)N0X,N1 X,ALX,BEX READ( 12,*)N0Y,N 1 Y,ALY,BEY READ(15,*) RιVlAXl,I TIPROF,HFERMLAFERxMI,BFERiVΪI,TSATITAiVrB,UMIREL,

1 NTNTER,NPRCON,NJX,NJY,TIMETERM,PLPP02 TSATI=TSATI+273.15 TAMB=TAMB+273.15 RMAX=RMAX1*2.54/100 TAU=RMAX/VRRIF

READ(14,*)TACC,OMEGA1,OMEGAF,TSPIN,EQCNTR TACC=TACC/TAU OMEGAI=OMEGAl*2*3.1415926/60 I IN RAD/S

OMEGAF=OMEGAF*2*3.1415926/60 ! IN RAD/S

TSPLN=TSPIN/TAU TIMETERM=TIMETERM/TAU FMH2OB=UMΓREL/100.*P0_H2O(TAMB)/1.* 18./29. c

C EXCHANGE COEFFICIENT CALCULATION

PROJ1 !ho calcolato per bene Pr tra 280-300 K

XMUAΓR=(0.175+5.E-4*(TAMB-280))* 1.E-4 !Pa*s

XKAIR=(0.0247+8.E-5*(TAMB-280))* l.E-3 !KW/(M*K) RE=(OMEGAF*4*RMAX**2*RHOGAS)/XMUATR

XNU1=0.88447*RE**.5*PR XNU2=0.62048*RE**.5*PR**0.33333 C C ASYMPTOTES COMBINATION XNU=XNUl*XNU2/(XNUl+XNU2)+1.8*((.3*2*RMAX*RHOGAS/XMUATR)*PR)**.3333

HTERM=XKALR*XNU/(2*RMAX) DWALR=0.01013*TAMB**lJ5*1.15957E-7 !M^Λ2/S DISALR=0.01013*TAMB**lJ5*4J48938E-8

SCW=XMUAΓR/RHOGAS/DWAΓR SCISO=XMUAIR/RHOGAS/DISAIR

SHW=1/(1/(0.88447*RE**.5*SCW)+1/(0.62048*RE**.5*SCW**0.33333)) 1 +1.8*((.3*2*RMAX*RHOGAS/XMUALR)*SCW)**.3333

SHISO=1/(1/(0.88447*RE**.5*SCISO)+1/(0.62048*RE**.5*SCISO**0.333)) 1 +1.8*((.3*2*RMAX*RHOGAS/XMUAIR)*SCISO)**.3333 XKC=SHW*DWAIR/(2*RMAX)

XKCISO=SHISO*DISAIR/(2*RMAX)

TAMB=TAMB/TRIF

TSATI=TSATI/TRIF TFΓN=TSPIN RETURN

END C

SUBROUTINE SINITO C C this element evaluates variables initializations

C NT ARE THE COLLOCATION POINTS (SEE DASSL.DAT)

C

IMPLICIT REAL(A-H,0-Z) INCLUDE 'spincom.f DTΓME=TFΓN/NINTER c c

NTX=NTX+N0X+N1X NTY=NJY+N0Y+N1Y

NTXM1=NTX-1

NTYM1=NTY-1

NTXM2=NTX-2

NTYM2=NTY-2 IF (EQCNTR.EQ.1)THEN

NPARZ=1

ELSELF (EQCNTR.EQ.2)THEN

NPARZ=3

ELSEIF (EQCNTR.EQ.3)THEN NPARZ=4

ENDIF

NNX=( 1 +NP ARZ)*NTX*NTY+NTX

TX=0.

RETURN

END C

SUBROUTINE OUTPAR(TX) C WRITING ON THE FILE OF THE CALCULATED POINT TX C the subroutine writes the calculation outputs at TX on file C 010 ('resO.res'). This file is the main results file. It C contains the full set of informations. C

IMPLICIT REAL(A-H,0-Z)

INCLUDE 'spincom.f

COMMON/MIA/NPRΓNT

CHARACTER* 10 STRTGG,STRTSS,STRTBB

NPRΓNT=NPRINT+I

C

IFILE=10 IFILE2=11 IFILE3=13 IFILE4=20 IFILE5=15 IFILE6=40 C

C PRINT ALL RESULTS ON RESULT.RIS C WRΠΈ(IFΓLE,*) (TX*TAU) NTEMP=NTX*NTY+NTY DO JY=l,NTY

WRITE(IFILE,*) (VR(JX,JY) * VRRIF ,JX= 1 ,NTX) ENDDO DO JY=l,NTY

WRITE(IFILE,*) (VZ(JX,JY)*VZRΓF,JX=1,NTX)

ENDDO

DO JY=l,NTY

WRITE(IFLLE,*) (CONC(JX,JY)*CRIF, JX=1 ,NTX) ENDDO

DO JY=l,NTY

WRITE(IFILE,*) (CISO(JX,JY)*CRLF,JX=l ,NTX) ENDDO DO JY=1,NTY WRITE(ΓFILE,*) (T(JX,JY)*TRΓF,JX=I,NTX)

ENDDO

WRITE(IFILE,*) (H(JX)*HMAX,JX=1,NTX) c

C PRINT THICKNESS ON SPESSORE.RIS c

WRITE(IFILE2,*) TX*TAU,(H(JX)*HMAX* 1.E6, JX=1 ,NTX)

C

C PRINT TEMPERATURE ON TEMPER.RIS C WRITE(IFILE3,*) TX*TAU,(T(JX,NTY)*TRIF-273.15,JX=1,NTX) C

C PRINT CONCENTRATIONS ON CONCEN.RIS C

WRITE(IFILE5,*) TX*TAU,(CONC(JX,NTY)*CRLF,JX=l ,NTX) C r

C PRINT HEATER VARIABLES

C

OPEN(UNIT=20,FILE='FORHEAT.RIS',STATUS='REPLACE')

DO 200 JY=l,NTY WRITE(IFILE4,*)(CONC(JX,JY),JX=l,NTX)

WRITE(IFILE4,*)(CONC 1 (JX,JY),JX=1 ,NTX) 200 CONTINUE

DO 250 JY=l,NTY

WRITE(IFILE4,*)(CISO(JX,JΥ) X=l,NTX) WRITE(IFILE4,*)(CISO 1 (JX,JY),JX=1 ,NTX)

250 CONTINUE

DO 300 JY=l,NTY

WRITE(IFILE4,*)(T(JX,JY),JX=1,NTX)

WRITE(IFILE4, *)(T 1 ( JX, JY), JX= 1 ,NTX) 300 CONTINUE WRITE(IFILE4,*) (H(JX),JX=1,NTX) WRITE(LFΓLE4,*) (H1(JX),JX=1,NTX) CL0SE(LFILE4) c C PRINT BIOT AND MARANGONI ON ADIM.RIS C

C TMED=(T(I,NTY)+T(I,l))/2

C DST=0.31843*l.E-3

C ALFA=XK/(RHO*CPLIQ) C MARANG=(-DST)*(T(I,NTY)-T(I,1))*H(I)*HMAX (ALFA*XMU(TMED))

WRITE(ΓFILE6,*) TX*TAU,REPL_;BIOT,BIOTISO,BIOTT

RETURN

END SUBROUTINE OKAY0

C THIS ROUTINE WRITES ON THE SCREEN THE SUCCESS OF THE

INTEGRATION

C

CHARACTER* 1 CICCIO INCLUDE 'SPINCOM.F

WRITE(*,*)

WRITE(*,*)

WRITE(*,*) 'VISUALSPIN 1.3BIS'

WRITE(*,*) WRITE(*,*) 'BY G.L.VALENTE & F. RUBINI'

WRITE(*,*)

WRITE(*,*) 'POLITECNICO Dl MILANO'

WRITE(*,*) 'DPT. CHIMICA-FISICA APPLICATA'

WRITE(*,*) 'MODELING OF SPIN COATING AND HEATING OPERATION TN 1 POLISHING PROCESS'

WRITE(*,*)

WRITE(*,*) 'The run has finished successfully! ! ' write(*,*) 'press 0 to continue'

READ(*,*) CICCIO RETURN

END C

SUBROUTINE START(INFO,IRES,RTOL,ATOL,RWORK,LRW_!NNEQ,NEQ) C C SDASPK initializations, according with instructions

C data are contained in file 021 ('dassl.dat') c DASSL.DAT CONTAINS SDASPK SETTINGS!! !! IMPLICIT REAL(A-H,0-Z) DIMENSION TNFO(20),RWORK(LRW) DIMENSION RTOL(NNEQ),ATOL(NNEQ) CHARACTER*9 STEP0,STEPM CHARACTER*8 YSIGN CHARACTER* 14 CONSYP CHARACTER* 12 JACNUM INCLUDE 'spincom.f

C first call

INFO(1)=0

C VECTOR error tolerances RTOL,ATOL

LNFO(2)=l

READ(20,*) RTOLV,ATOLV

15 NN2=NTX+NTX*NTY

DO 1 J=1,(NTX*NTY)

RTOL(J)=RTOLV

ATOL(J)=ATOLV

RTOL(J+NN2)=RTOLV 20 ATOL(J+NN2)=ATOLV

1 CONTINUE READ(20,*) RTOLH,ATOLH NN1=NTX*NTY

D0 2 J=1,NTX _z5 RTOL(J+NNl)=RTOLH

ATOL(J+NNl )=ATOLH

2 CONTINUE READ(20,*) RTOLC TOLC NN3=NN2+NTX*NTY

₃₀ D0 3 I=1,(2*NTX*NTY)

RTOL(J+NN3)=RTOLC ATOL(J+NN3)=ATOLC

3 CONTINUE READ(20,*) RTOLT,ATOLT

„ -. NN4=NN3+2*NTX*NTY

^ DO 4 J=l,(NTX*NTY)

RTOL(I+NN4)=RTOLT

ATOL(I+NN4)=ATOLT

4 CONTINUE

40 C writing density index

READ(21 ,*)LDST,IVrDEO

C NO intermediate integration stops

Q the integration proceed to TOUT

INFO(3)=0 45 C integration can proceed beyond TSTOP ΓNFO(4)=O

C evaluation of the jacobian

READ(21,*) JACNUM

IF(JACNUM.EQ.'NUMERIC JAC) INFO(5)=0 IF(IACNUM.EQ.'ANALYTIC JAC) INFO(5)=l

C Jacobian structure = full matrix

LNFO(6)=0 C decision for step max

READ(21,*) STEPM,HM O IF(STEPM.EQ.'AUTO STMX') INFO(7)=0

ΓF(STEPM.EQ.ΈXT STMX') THEN ΓNFO(7)=I

RW0RK(2)=HM

ENDΓF 5 C decision for 1st step

READ(21,*) STEP0,H0 ΓF(STEPO.EQ.ΆUTO ISTP*) ΓNFO(8)=O IF(STEP0.EQ.'EXT ISTP') THEN ΓNFO(8)=I 0 RWORK(3)=H0

ENDIF C set the maximum integration order to default MAXORD = 5

LNFO(9)=0

Q restraint on non negativity of Y(I) 5 READ(21,*) YSIGN

IF(YSIGN.EQ.'Y.GE.0.0') TNFO(10)=1

IF(YSIGN.EQ.'Y.LT.O.O') INFO(10)=0 C G(TX,Y,YPRIME) = 0 (consistent) at TX=0

READ(21,*) CONSYP 0 IF(CONSYP.EQ.'YPRIME YES CON') INFO(l 1)=0

IF(CONSYP.EQ.'YPRIME NOT CON') TNFO(ll)=l C newton iteration with a direct method

INFO(12)=0

INFO(13)=0

35 C

ΓNFO(15)=O c

IRES=0 c

40 RETURN

END c

SUBROUTINE COLLOC

C

-;5 C the subroutine coordinates the call necessaries to the C orthogonal collocation method. Such a method is suitable for

C the solution of PDE.

C The subroutine evaluates the zeros of the Jacobi polynomials

C besides the first and the second derivatives of the C Lagrange polynomials. Integration interval [0,1].

C

C NB: convention adopted for Lagrange polynomia:

C

C ALl(point,polynomia) C AL2(point,polynomia)

C c

INCLUDE 'spincom.f C DIMENSION VETl (NTPC),VET2(NTPC),ROOT(NTPC)

DIMENSION DIF1 (NTPC),DΓF2(NTPC),DIF3(NTPC) DIMENSION VRF(NTPC) C

NRMAX=NTPC ND=NTPC

ZERO=1.E-20 LD=3

C jacobi's zeros & derivatives along x direction

CALL JCOBI( ro,NJX,N0X,NlX_J rTX,ALX,BEX,DIFl,DIF2,DIF3,ROOT) DO 1 1=1, NTX

ROOTX(I)=ROOT(I)

CALL DFOPR(ND,NTX,I,DIF 1 ,DIF2,DIF3,ROOT, VETl , VET2) DO 2 K=1,NTX AL1X(I,K)=VET1(K) AL2X(I,K)=VET2(K)

2 CONTINUE

CALL RAD AU(ND-NTY,N0 Y,N 1 Y,ID, ALY,BEY,ROOT,DIF 1 NET 1 )

DO 30 K=1,ΝTX

WQX(I,K)=VET1(K) 30 CONTINUE

1 CONTINUE

C c

C test on derivatives calculations C DO 30 I=l,NTX

C VRFI=.0

C DO 40 K=l,NTX

C VRFI=VRFI+AL1X(I,K)

C40 CONTINUE C VRF(I)=VRFI

C30 CONTINUE

C c C jacobi's zeros & derivatives along y direction

CALL JCOBI( uD,NJY,N0Y,NlY,NTY,ALY,BEY,DrFl,DIF2,DIF3_>ROOT) DO 3 1=1, NTY ROOTY(I)=ROOT(I)

CALL DF0PR(ND,NTY,I,DIF1,DIF2,DIF3,R00TNET1NET2) D0 4 K=1,ΝTY

AL1Y(I,K)=VET1(K) AL2Y(I,K)=VET2(K)

4 CONTINUE c Lobatto quadrature (0,1) CALL RADAU(ND,NTY,N0Y_>N1 Y,LD,ALY,BEY,R00T,DLF1NET1)

DO 5 K=l,ΝTY WQY(I,K)=VET1(K)

5 CONTINUE

3 CONTINUE C printing of the calculations

LFILE=10

IF(IDST.GE.3) THEN WRITE(ΓFΓLE,50)

WRITE(IFΓLE,51) WRITE(TFLLE,217)' — x direction ~

',NJX,ALX,BEX

WRITE(LFILE,218)(ROOTX(L),L=l,NTX) WRITE(IFILE,219) DO 7 1=1, NTX WRITE(IFILE,220)(AL1X(I,K),K=1,NTX)

7 CONTINUE

WRITE(IFILE,221) DO 8 1=1, NTX WRITE(IFILE,220)(AL2X(I,K),K=1,NTX) 8 CONTINUE

C WRITE(IFILE₎52)(VRF(I),I=1,NTX)

ENDIF

IF(IDST.GE,3) THEN

WRITE(IFTLE,217)' — y direction - ',NJY,ALY,BEY WRITE(IFILE,218)(ROOTY(L),L=l,NTY)

WRITE(IFILE,219) DO 9 1=1, NTY WRITE(IFILE,220)(AL1Y(I,K),K=1,NTY) 9 CONTINUE WRITE(IFILE,221) DO 10 1=1, NTY WRITE(TFILE,220)(AL2Y(I,K),K=1 ,NTY) 10 CONTINUE C WRITE(IFILE,52)(VRF(I),I=1,NTX)

ENDIF C WRITE(IFILE,50)

C

50 F0RMAT(1XJ9('*')) 51 F0RMAT(/1X,'*** ORTHOGONAL COLLOCATION METHOD PARAMETERS: ' + '(SUBROUTINE COLLOC) *****^•)

217 FORMAT(/lX,'** PARAMETERS OF THE JACOBI POLYNOMIAL **', + 5X,A21,

+ /1X,'DEGREE =',12, IX,' W = (( 1 -X)**ALPHA)*X**BETA', + /1X,'ALPHA = ',F5.2,5X,'BETA = ',F5.2)

218 FORMAT(/lX,'ZEROS; ',10(1X,1PE10.3))

219 FORMAT(/lX,'LAGRANGE POLYNOMIAL FIRST ORDER DERIVATIVES',/)

220 F0RMAT(1X,9(1PE11.4,2X))

221 F0RMAT(/1X, 'LAGRANGE POLYNOMIAL SECOND ORDER DERTVATINES',/) 52 F0RMAT(1X,'VRF»»»»',8(E11.4,2X))

C

RETURN END

C C c

SUBROITTINE JCOBI(ND,N,NO,N1,NT,AL,BE,DIF1,DIF2,DIF3,ROOT) DIMENSION DIFl(NT3),DIF2(NO),DIF3(ND)_IROOT(ND)

C AB=AL+BE

AD=BE-AL AP=BE*AL

DIF 1 ( 1 )=(AD/(AB+2)+l)/2 DTF2(1)=0. IF(N.LT.2) GO T0 15

DO 10 I=2,N Z1=I-1 Z=AB÷2*Z1

DIFl(I)=(AB*AD/(Z*(Z+2))+l)/2 IF(I.NE.2) GO TO 11

DIF2(I)=(AB+AP+Z 1 )/(Z*Z* (Z+l)) GO TO 10 11 Z=Z*Z

Y=Z1*(AB+Z1) Y=Y*(AP+Y)

DΓF2(I)=Y/Z/(Z-I) CONTINUE X=0. L=l

DO 20 1=1, N XD=0.

XN=1. XD1=0. XN1=0.

DO 30 J=1,N XP=(X-DIF 1 (J)) *XN-DJJ2(J)*XD XP 1 =(X-DTF 1 (J))*XN 1 -DΓF2(J)*XD 1 +XN XD=XN XD1=XN1

XN=XP XN1=XP1 CONTINUE ZC=1. Z=XN/XN1

ΓF(I.EQ.I) G0 T0 21

DO 22 J=2,l ZC=ZC-Z/(X-ROOT(J-l)) Z=Z/ZC X=X-Z

LF (L.GT.500) GO TO 1 L=L+1

IF (ABS(Z/X).GT.l.E-5) GO TO 25 ROOT(I)=X X=X+.0001 CONTINUE CONTINUE

IFΓNO.EQ.O) GO TO 35

DO 31 1=1, N J=N+1-I ROOT( J+ 1 )=ROOT( J) ROOT(1)=0. IF(Nl.EQ.l)ROOT(NT)=l.

DO 40 1=1, NT X=ROOT(I)

DIF1(I)=1.

DIF2(I)=0.

DIF3(I)=0.

DO 40 J=1,NT IF(J.EQ.I) GO TO 40 Y=X-ROOT(J)

DIF3(I)=Y*DIF3(I)+3*DTF2(I)

40 CONTINUE RETURN END

C C

SUBROUTrNE DFOPR(ND,NT,I,DIFl,DIF2,DIF3,ROOTNETlNET2) c

DIMENSION DIFl(ΝD)₅DIE2(ΝD),Drf3(ΝD)_;ROOT(ΝD)NETl(ΝD)NET2f ro) DO 20 J=l,ΝT

LF(J.NE.I) GO T0 21 VETI(I)=DΓF2(I) DΓFI(I)/2 VET2(I)=DIF3(I)/DΓF1(I)/3

GO TO 20

21 CONTINUE

Y=ROOT(I)-ROOT(J) VETI(J)=DΓPI(I)/DΓFI(J)/Y VET2(J)=VET1(J)*(DΓF2(I)/DIF1(I)-2/Υ)

20 CONTINUE RETURN END

C c c

SUBROUTINE RADAU(ND,NT,N0,Nl,ro,AL,BE,ROOT,DIFl,VET)

IMPLICIT REAL(A-H,0-Z) DIMENSION DrFl(ΝD),ROOT(ΝD)NET(ΝD)

C ID = 3 LOBATTO QUADRATURE (BOTH 0 AND 1 ENDPOINTS) C

S=0.

DO 41 1=1, NT X=ROOT(I)

IF(ID-2) 10,20,30 10 AX=X

IF(N0.EQ.0) AX=1./AX GOTO 40 20 AX=1.-X

IF(N1.EQ.0) AX=1./AX GOTO 40 30 AX=l. 40 VET(I)=AX/DIF1(I)**2. 41 CONTINUE

IF(ID.NE.2) VET(NT)=VET(NT)/(1.+AL)

IF(ID.GT.l) VET(1)=VET(1)/(1.+BE)

DO 50 1=1, NT S=S+VET(I) 50 CONTINUE

DO 60 1=1, NT VET(I)=VET(I)/S 60 CONTINUE

RETURN

END

C THE FOLLOWING SET SUBROTINES ARE RELATED TO SDASPK SOLVER.

C THEY MAINLY ANALYZE THE OUTPUTS, THE ERROR MESSAGES AND TRY C TO RECOVER THE INTEGRATION LF MINOR FAILURES HAPPEN.

C C

PSOL(NEQ,T,Y,YPRIME,SAVR,WK_!CJ,WGHT,WP,IWP,B,EPLLN, + ΓER,RPAR,IPAR)

IMPLICIT REAL(A-H,0-Z) DIMENSION Y(NEQ),YPRIME(NEQ) RETURN

END

C c

SUBROUTINE SDHELP(IDID,IFIL) IMPLICIT REAL(A-H,0-Z)

WRITE(IFIL,30) IF(IDID.GT.O) WRITE(IFIL,49) IF(IDID.EQ.l) WRITE(IFIL,51) LF(IDID.EQ.2) WRITE(IFIL,52) IF(IDID.EQ.3) WRITE(IFIL,53)

IF(IDID.LT.O) WRITE(IFIL,50) IF(IDID.EQ.-l) WRITE(IFIL,54) IF(IDID.EQ.-2) WRITE(TFIL,55) IF(IDID.EQ.-3) WRITE(IFIL,56) IF(IDID.EQ.-5) WRITE(IFIL,57)

IF(IDID.EQ.-6) WRITE(IFIL,58) IF(IDID.EQ.-7) WRITE(IFIL,59) IF(IDID.EQ.-8) WRITE(IFIL,60) LF(IDID.EQ.-9) WRITE(IFIL,61)

IF(IDID.EQ.-IO) WRITE(TFIL,62)

IF(IDID.EQ.-11) WRITE(IFIL,63)

IF(IDID.EQ.-12) WRITE(ΓFIL,64) ΓF(IDID.EQ.-13) WRITE(IFIL,65)

IF(IDID.EQ.-H) WRITE(IFIL,66)

IF(IDLD.EQ.-33) WRITE(IFIL,67) C

30 F0RMAT(1XJ8('*')) 49 F0RMAT(1X,'TASK COMPLETED')

50 FORMAT(lX,'TASK INTERRUPTED*)

51 FORMAT(lX,'succesfully intermediate step, T<TOUT EDLD = 1 ')

52 FORMAT(lX,'integration succesfully completed T=TOUT LDID = 2')

53 FORMAT(lX,'integration succesfully completed T>TOUT IDID = 3') 54 FORMAT(lX, 'too large amount of work (500 steps) LDID = -1')

55 FORMAT(lX,'tolerances too stringent LDID = -2')

56 FORMAT(lX,'local error test not satisfied ATOL=0 , IDID = -3')

57 FORMAT(lX,'repeated failures in JAC IDLD = -5')

58 FORMAT(lX,'SDASPK, error test failures, last step IDID = -6') 59 FORMAT(lX,'non linear system solv. could not convege LDLD=-7')

60 FORMAT(lX, 'matrix of partial derivatives singular LDLD = -8')

61 FORMAT(lX,'non linear system solver not converges LDLD = -9')

62 FORMAT(lX,'non linear solver fails because LRES=-1 LDLD =-10')

63 FORMAT(lX,TRES=-2 encountred, intergration halted LDLD =-11') 64 FORMAT(lX,'SDASPK fails to compute YPRIME at T=0 LDLD =-12')

65 FORMAT(lX,'unrecoverable error in PSOL routine LDID =-13')

66 FORMAT(lX,*KRYLOV linear system solver fails IDID =-14')

67 FORMAT(lX,'the code has encontred fatal errors IDID =-33') C RETURN

END

C

C

SUB ROUTINE SDTEST(RWORK,LRW,IWORK,LIW) IMPLICIT REAL(A-H,0-Z)

DIMENSION RWORK(LRW),IWORK(LIW)

C step size used in the last step

HOLD=RWORK(7) order of the method for the next step KORD=IWORK(7)

Q order of the method used in the last step

KΓLD=ΓWORK(8) f number of steps taken so far

NST=IWORK(l l) C number of calls of RES so far NRE=IWORK(12) C number of calls of JAC so far

NJE=IWORK(13)

C number of error test failures so far NETF=IWORK(14)

C number of non linear conv. failures so far

NCFN=IWORK(15) Q number of conv. failures in linear iter.ns

NCFL=IWORK(16) C length oflWORK actually required

LENIW=IWORK(17) C length of RWORK actually required

LENRW=IWORK(18)

C total number of non linear iterations NNI=IWORK(19)

RETURN END

C c

SUBROUTINE CRISIS(NEQ,LDID,LNFO,Y) C

C this routine try to solve a not fatal crisis of the solver C using the SDASPK instructions

IMPLICIT REAL(A-H,0-Z) DIMENSION INFO(20)N(NEQ)

LF(IDID.EQ.-l) THEN INFO(l)=l

WRTTE(*,*)'about 500 steps were taken, go fo 1000'

RETURN

ENDIF

IF0DID.EQ.-2) THEN LNFO(l)=l

WRITE(*,*)'tolerances _{were re}ieased, continue'

RETURN ENDIF IF(IDID.EQ.-3) WRITE(*,*)'SET ALL ATOL.NE.0.' IF(IDID.EQ.-5) WRITE(*,*)'JAC fails with the KRYLOV method'

IF(IDID.EQ.-6) WRITE(*,*)'error in last step (SINGULARITY)' IF(IDID.EQ.-7) WRITE(*,*)'ill-conditioned jacobian, restart' IF(IDID.EQ.-8) WRITE(*,*)PD is singular, redundant equations¹ IF(IDID.EQ.-9) WRITE(*,*)'iH-posed problem or discontinuity' IF(IDID.EQ.-IO) WRITE(*,*)TRES=-1' ΓF(ΓDID.EQ.-H) WRITE(*,*)'IRES=-2' LF(ΓDΓD.EQ.-I2) THEN

WRITE(*,*)'not able to evaluate YPRIME, G(Y,YPRLME,T).NE.O' C LFILE=22 C REWIND 22

C CALL DCBA(NEQ,Y)

C CALL SAVEIT(NEQ,LFILE,Y)

WRITE(22,*)^,G(Y, YPRIME, Y).GT.TOL - BEST SOLUTION: RUN AGAIN' C CLOSE (22) ENDIF

LF(LDID.EQ.-13) WRITE(*,*)'fatal error in PSOL' LF(IDLO.EQ.-14) WRITE(*,*)'KRYLOV fails' LF(LDID.EQ.-33) WRITE(*,*)'integration cannot continue' C STOP c

C RETURN END c SUBROUTINE JAC(TX,Y,YPRLME,PD,CJ,RPAR,NEQ)

C

C PD(ij) = d G(i)/d Y(j) + CJ* d G(i)/d YPPJME(j)

C

C Gi(T,Y,YPRIME)=0. C

IMPLICIT REAL(A-H,0-Z)

INCLUDE 'spincom.f

DIMENSION Y(NEQ),YPRIME(NEQ),PD(NEQ,NEQ) C

C B.C. ON Z=0

DO 101=1, NTX

10 CONTINUE

DO 12 1=1, NTX IEQ=I+NTX

IV AR=NTX*NTY+NTX+NT Y* (I- 1 )+ 1 PD(IEQ,IVAR)=1. 12 CONTINUE

NDELMAX=2*NTX

C B.C. ON R=0

DO 14 J=2,NTYM1 IEQ=(J-1)+NDELMAX IVAR=J PD(IEQ,IVAR)=1. IEQ=(J- l)+NDELMAX+(NTY-2) DO N=l,NTX

IN R=NTX*NTY+NTX -(N-1)*NTY+J PD(IEQ,IVAR)=AL1X(1,N)

ENDDO 14 CONTINUE

NDELMAX=NDELMAX+2*(NTY-2) C EQUAZIONE SULL'ANGOLO R=0 Z=H IEQ=NDELMAX+1

DO N=l,NTX

IVAR=NTX*NTY+NTX+(N-1)*NTY+NTY PD(IEQ,IVAR)=AL1X(1,N) ENDDO NDELMAX=NDELMAX+1

C B.C. ON Z=H Dl VRZ=0

DO 16 1=1, NTX IEQ=I+NDELMAX DO M=l,NTY IVAR=M+(I-1)*NTY

PD(IEQ,ΓVAR)=ALIY(NTY,M)

ENDDO 16 CONTINUE

NDELMAX=NDELMAX+NTX C B.C. SU R=l

DO 18 J=2,NTYM1 IEQ=NDELMAX+(J-1) DO N=l,NTX IVAR=J+(N-1)*NTY PD(IEQ,IVAR)=AL1X(NTX,N)

ENDDO 18 CONTINUE

NDELMAX=NDELMAX+NTYM2

C B.C. SU D1CZ=0. Z=0. DO 1=1, NTX

IEQ=NDELMAX+I DO M=l,NTY

IVAR=2*NTX*NT Y+NTX^(I- 1 )*NTY+M PD(IEQ,IV AR)=AL 1 Y( 1 ,M) ENDDO

ENDDO NDELMAX=NDELMAX+NTX

C B.C. R=0. D1CR=0.

DO J=2.NTYM1 IEQ=NDELMAX÷J-1 DO N=l,NTX

IVAR=2*NTX*NTY+NTX+(N-1)*NTY+J PD(IEQ_;IVAR)=AL1X(1,N) ENDDO ENDDO

NDELMAX=NDELMAX+NTYM2

C B.C. DlCZ=Kc DIFNH3*HMAX*H*DELC LN Z=l

DO I=1,NTX IEQ=NDELMAX+I rVAR=2*NTX*NTY+NTX+(I- 1)*NTY+NTY

PD(IEQ,IVAR)=XKC DIFNH3*HMAX*H(I) IVAR=NTX*NTY+I

PD(IEQ,IVAR)=XKC/DIFNH3*HMAX*CONC(I,NTY) ENDDO c CONTINUITY EQUATION

DO 20 I=2,NTX

R=ROOTX(I)

DO 22 J=2,NTY

LEQ=NDELMAX+(J-2)*NTXM1 +1- 1 C DER. WITH RESPECT TO VR(I,J) rVAR=J+(I-l)*NTY

PD(IEQ,IVAR)=VRRIF RMAX*(l R+ALlX(I,r))

DO N=l,NTX

LF (N.NE.I) THEN IVAR=J+(N-1)*NTY

C DER. WITH RESPECT TO VR(N,J)

PD(IEQ,rVAR)=VRRLF RMAX*ALlX(I,N)

ENDIF

ENDDO C DER. WITH RESPECT TO H(I)

ΓVAR=NTX*NTY+I

D1VZZ=0. DO M=l,NTY

D1VZZ=D1VZZ+AL1Y(J,M)*VZ(LM) ENDDO

PD(IEQ,IVAR)=-(VZRLF/(H(I)*H(I))/HMAX*D 1 VZZ) C DER. WITH RESPECT TO VZ(I,M)

DO M=l,NTY

IV AR=NTX*NT Y+NTX+M+(I- 1 )*NTY PD(IEQ,IVAR)=VZRTF/H(I)/HMAX*AL1Y(J,M)

ENDDO 22 CONTINUE 20 CONTINUE

NDELMAX=NDELMAX+NTYM1 *NTXM1 C NAVIER-STOKES EQUATION

DO 24 I=2,NTXM1 DO 26 J=2,NTYM1

IEQ=NDELMAX+( J-2)*(NTX-2)+(I- 1 ) ΓVAR=(I-I)*NTY+J

C DER. WITH RESPECT TOVRL(I )

PD(IEQ,IVAR)=-CJ*VRJΛIE+XMU0/RHO*VRRIF/((HMAX*HMAX)*(H(I)*H(I)))*

1 AL2Y(J,J) C DER. WITH RESPECT TO VR(I,M) DO M=l,NTY

IVAR=(I-1)*NTY+M IF (M .NE. J) THEN

PD(IEQ,IVAR)=XMU0/PJiO*VRRIF/(HMAX**2*H(I)**2)*AL2Y(J,M) ENDIF ENDDO C DER. WITH RESPECT TO H(I) IVAR=NTX*NTY+I D2VRZ=0. DO M=l,NTY D2VRZ=D2VRZ+AL2Y(J,M)*VR(I,M)

ENDDO PD(IEQ,IVAR)=-2*XMU0/RHO*VRRIF/(HMAX**2*H(I)**3)*D2VRZ

26 CONTINUE 24 CONTINUE NDELMAX=NDELMAX+NTXM2*NTYM2

C

DO 28 1=1, NTX

IEQ=NDELMAX+I ΓVAR=NTX*NTY+I C DER. WITH RESPECT TO H(I)

PD(IEQ,IVAR)=-AL1X(I,I)*VR(I,NTY)*VRRIF*HMAX RMAX+CJ*(-HMAX)

DO N=1,NTX IF (N.NE.I) THEN IVAR=NTX*NTY+N C DER. WITH RESPECT TO H(N)

PD(IEQ,IVAR)=-AL 1 X(I,N)* VR(I,NTY)* VRRIF *HMAX RM AX

ENDIF ENDDO C DER. WITH RESPECT TO VR(I,NTY) IVAR=(I-1)*NTY+NTY

D 1HR=0. DO N=1,NTX D1HR=D1HR^AL1X(LN)*H(N) ENDDO

PD(IEQ,IVAR)=-D 1 HR* VRRIF*HMAX RMAX C DR. WITH RESPECT TO VZ(I,NTY)

ΓVAR=NTX*NTY+NTX+(I-I)*NTY+NTY PD(IEQ,IVAR)=VZRLF

28 CONTINUE

DO IEQ=l,NEQ

WRITE(30,*) (PD(IEQ,rVAR),IVAR=l,NEQ) ENDDO RETURN

END

FUNCTION P0_H2O(TEMP) ! PO IN ATM P0_H2O=(EXP(18.3036-3816.44/(TEMP-46.13)))/760. END

FUNCTION PO SO(TEM) ! PO IN ATM C 1=76.964 C2=7623.8 C3=7.4924

C4=5.9436E-18

P0_ISO=EXP(C 1 -C2/TEM-C3 *LOG(TEM)+C4*TEM* *6.)/101500.

END FUNCTION XMU(TEMP,CONCEN)

! VISC. IN Pa*S INCLUDE 'SPLNCOM.F' C CONCEN=CONCEN+100.

A=-8.430086958099 B=6.1228655 C=38.007371

XMU=EXP(A+B/(CONCEN/RHO))*EXP(C/(TEMP-273.15))*l.E-3 END FUNCTION FERMI(X) !PROFILO INIZIALE

INCLUDE 'SPINCOM.F

FERMI=l.E-3+HFERMI/(EXP((X-AFERMI) BFERMI)+l) END FUNCTION OMEGA(TIME)

INCLUDE 'SPLNCOM.F' C IF (TIME.LT.TSPLN) THEN

IF (TIME.LT.TACC) THEN

OMEGAOME GAI+(OME GAF-OMEGAI)/TAC C *TLME ELSE OMEGA=OMEGAF ENDIF C ELSE

C OMEGA=OMEGAF*EXP(-0.5*(TLME-TSPLN-TACC))+0.001

C OMEGA=0.001

C ENDIF

END c

FUNCTION TSAT(TLME,NI) INCLUDE 'SPLNCOM.F

TSAT=(TSATI-T0(NI))*TLME/TLMETERM+T0(NI) END

C

FUNCTION GAMMA(XCOMP) C

A=0.299753

B=-13.224378

GAMMA=EXP((l-XCOMP)**2*(A+2*(B-A)*XCOMP))

END

Bibliography

Acrivos A., Shah M.J., Petersen E.E., J.Appl.Phys 31, 963, (1960)

Bird R.B., Stewart W.E., Lightfoot E.N., Transport phenomena, John Wiley & Sons, Inc, (1960)

Bomside D.E., Macosko C.W., Scriven L.E., J. imaging Technol. 13, 122, (1987)

Britten J.A., Thomas I. , J. Appl. Phys. 71(2), 972, (1992)

P.N. Brown, AC. Hindmarsh, and L.R. Petzold, Using Krylov Methods in the Solution of Large- Scale Differential-Algebraic Systems, SIAM J. Sci. Comp., 15 (1994), pp. 1467-1488

P.N. Brown, AC. Hindmarsh, and L.R. Petzold, Consistent Initial Condition Calculation for Differential-Algebraic Systems, LLNL Report UCRL-JC-122175, August 1995; submitted to SIAM J. Sci. Comp., 1995

Dente M., Ranzi E., Principi di Ingegneria Chimica, CLUP, Milano, (1978)

Emslie A.G., Bonner F.T., Peck L.G., J.Appl.Phys., 29, 858, (1958)

Flack W.W., Soong D.S., Bell AT., Hess D.W., J.Appl.Phys 56, 1 99, (1984)

Givens F.L., Daughton W.L, J. Electrochem. Soc. 126, 269, (1979)

Givens F.L., Daughton W.L., J. Electrochem. Soc. 129, 173, (1982)

Higgins B.G., Phys. Fluids 29, 3522, (1986)

Jenecke S.A., Ind. Eng. Chem. Fundam. 23, 425, (1984)

Jenecke S.A., Polym. Eng. Sci. 23, 830, (1983)

Jenecke S.A., Schuldt S.B., Ind. Eng. Chem. Fundam. 23, 432, (1984) Kreith F„ Taylor J.H., Chong J.P., J. Heat Transfer, 81 , 95, (1959)

Meyerhofer D., J. Appl. Phys. 49, 3993, (1978)

Monti S., Proprieta reologiche e volumetriche di serie omologhe oligomehche, Tesi di Laurea in Ingegneria Chimica del Politecnico di Milano (1996-1997)

Ohara T., Matsumoto Y., Ohashi H., Phys. Fluids 1, 1949, (1989)

Perry R.H., Green D., Perry's Chemical Engineering Handbook - 50^th Ed., McGraw-Hill, (1984)

Petzold L.R., Brown P.N., Hindmarsh A.C., Ulrich C.W., SDASPK Routine, (1990)

Regh T.J., Higgins B., AlChe J. 38 (4), 489, (1992)

Regh T.J., Higgins B.G., Phys. Fluids 31, 1360, (1988)

Reid R.C., Prausnitz J.M., Poling B.E., The Properties of Gases and Liquids, McGraw-Hill, (1988)

Schlichting H., Boundary-Layer Theory, McGraw-Hill, (1979)

Sparrow E.M., Gregg J.L, J. Heat transfer 81 , 294, (1960)

Sparrow E.M., Gregg J.L., J. Heat transfer, ASME, Series C, 81 , 249, (1959)

Sukanek P.C., J. Imaging Tβchnol. 11 , 184, (1985)

Tanner R.I., Engineering rheology, Oxford, (1985) Villadsen J., Michelsen M.L., Solution of Differential Equation Models by Polynomial Approximation, Prentice-Hall, (1978)

von Karman Th, Z. Angew.Math.Mech. 1, 233, (1921)

In view of the above, it will be seen that the several objects of the invention are achieved and other advantageous results attained.

When introducing elements of the present invention or the preferred embodiment(s) thereof, the articles "a", "an", "the" and "said" are intended to mean that there are one or more of the elements. The terms "comprising", "including" and "having" are intended to be inclusive and mean that there may be additional elements other than the listed elements.

As various changes could be made in the above constructions without departing from the scope of the invention, it is intended that all matter contained in the above description or shown in the accompanying drawings shall be interpreted as illustrative and not in a limiting sense.

Claims

ClaimsWHAT IS CLAIMED IS:

1. A method of forming thin layers of wax on surfaces of mounting blocks for mounting semiconductor wafers on the mounting blocks so that each wax layer has a desired wax layer shape and thickness selected to facilitate shaping of each wafer during subsequent polishing of the wafer, the method comprising the steps of: a. selecting values of at least some variables having an effect on the shape and thickness of the wax layer, b. modeling the wax layer numerically multiple times while changing at least ^u one of said values to determine values of said variables for said desired wax layer shape and thickness, c. forming wax layers on the surfaces of the mounting blocks by using said determined values of said variables to achieve said desired wax layer shape and for producing a consistent, desired shape of the wafers after polishing.

2. A method as set forth in claim 1 wherein the step of forming includes depositing wax on the surface and rotating the block.

3. A method as set forth in claim 2 wherein said at least one variable comprises final angular velocity of the mounting block, a change in the value of the final angular velocity being capable of affecting the predicted wax layer shape and thickness, said step of modeling the wax layer comprising changing the final angular velocity to determine the value of final angular velocity for the desired shape and thickness, and wherein the step of forming wax layers comprises controlling final angular velocity at said determined value to achieve the desired wax layer shapes and thicknesses.

4. A method as set forth in claim 2 wherein said at least one variable comprises time for spinning the block, a change in said spinning time being capable of affecting the predicted wax layer shape, said step of modeling the wax layer comprising changing the spinning time to determine the value of acceleration time for the desired shape and thickness, and wherein the step of forming wax layers comprises selecting the spinning time to achieve the desired wax layer shapes.

5. A method as set forth in claim 1 wherein said modelling step produces a predicted roll-off value of the predicted wax layer shape.

6. A method as set forth in claim 1 wherein said at least one variable ^ϋ includes a thermal gradient in the block measured from a center of the block to an edge of the block, a change in the thermal gradient being capable of affecting the predicted wax layer shape, said step of modeling the wax layer comprising changing the thermal gradient of the block to to determine the value of thermal gradient for the desired shape and thickness, and wherein said step of forming wax layers on the mounting blocks includes 5 controlling the thermal gradient to achieve the desired wax layer shapes.

7. A method as set forth in claim 1 wherein said at least one variable includes a block temperature, a change in the block temperature being capable of affecting the predicted wax layer shape, said step of modeling the wax layer comprising changing the block temperature while holding the other of said variables at the determined ⁰ values to to determine the value of block temperature for the desired shape and thickness, and wherein the step of forming the wax layers comprises controlling the block temperature value to achieve the desired wax layer shapes.

8. A method as set forth in claim 1 wherein the step of modeling the wax layer is performed using a computer processor.

5 9. A method as set forth in claim 1 wherein the step of modeling the wax layer integrates a mass conservation equation in a two dimensional domain using orthogonal collocation.

10. A method of predicting the shape and thickness of a wax layer formed on a mounting block for mounting a semiconductor wafer on the block, the method comprising the steps of: mounting the wafer on the mounting block, coating a surface of the wafer with a lubricating fluid, selecting values of at least some variables having an effect on the shape and thickness of the wax layer, depositing wax on a surface of the wafer and spinning the wafer using said selected values to form a predictive wax layer on the wafer, - demounting the silicon wafer from the mounting block, and measuring the shape and thickness of the predictive wax layer to predict the shape and thickness of said wax layer formed on the mounting block at said selected values.

11. The method of claim 10 further comprising the steps of removing the wax layer from the wafer, measuring the shape and thickness of the wafer and comparing the shape and thickness of the wax layer and the wafer.