WO2002069391A1  Method of forming a thin wax layer on a mounting block for mounting a semiconductor wafer  Google Patents
Method of forming a thin wax layer on a mounting block for mounting a semiconductor wafer Download PDFInfo
 Publication number
 WO2002069391A1 WO2002069391A1 PCT/IT2001/000094 IT0100094W WO02069391A1 WO 2002069391 A1 WO2002069391 A1 WO 2002069391A1 IT 0100094 W IT0100094 W IT 0100094W WO 02069391 A1 WO02069391 A1 WO 02069391A1
 Authority
 WO
 Grant status
 Application
 Patent type
 Prior art keywords
 wax layer
 shape
 wax
 block
 thickness
 Prior art date
Links
Classifications

 B—PERFORMING OPERATIONS; TRANSPORTING
 B24—GRINDING; POLISHING
 B24B—MACHINES, DEVICES, OR PROCESSES FOR GRINDING OR POLISHING; DRESSING OR CONDITIONING OF ABRADING SURFACES; FEEDING OF GRINDING, POLISHING, OR LAPPING AGENTS
 B24B37/00—Lapping machines or devices; Accessories
 B24B37/34—Accessories
 B24B37/345—Feeding, loading or unloading work specially adapted to lapping
Abstract
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 singleside polishing.
A conventional process for polishing a single side of a semiconductor wafer includes "wax mounting" the wafer on a diskshaped 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 singleside 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_{0}, ω, 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 NonNewtonian fluids Acrivos, Shah, and Petersen (1960) extended the analysis of Emslie et al. (1958) to the flow of a nonNewtonian 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 nonNewtonian fluids. Janecke and Shuldt studied the coating flow of nonNewtonian 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 Newtonianlike behavior at low shear rates while the Bingham and powerlaw 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 masstransfer 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 filmaveraged equations. Flack et al. (1984) numerically simulated the process taking into account the nonNewtonian 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 nonNewtonian behavior. However, many simplifications, such as using a constant evaporation rate, a constant solute concentration at the surface, and filmaveraged 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 gasphase 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 onedimensionai 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. Highspeed 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_{0}
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. 47 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 A1 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 noncontact 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 nonuniformity 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 ith component
Jι mass diffusive flux
0, diffusivity of the ith component in the liquid mixture radial wax velocity
Vz axial wax velocity _{3} azimuthal (tangential) wax velocity r radial coordinate z axial coordinate s azimuthal (tangential) coordinate ω angular velocity k_{d}is_{k} thermal conductivity of the disk δ_{d}isk half of the disk's thickness i_{h}erm heat transfer coefficient λ_{ev}, latent heat of evaporation of the ith component
T_{b} ^{s} temperature in the gas phase bulk k_{c},ι mass transfer coefficient
Co, '^{s} mass fraction of the component ith in the gas phase at the gasliquid interface ω,^{*}'' mass fraction of the component ith in the liquid phase at the gasliquid interface ω,^{ά,s} mass fraction of the component ith in the gas phase bulk
ω,^{w} mass fraction of the component ith in the liquid phase bulk
P°(T) vapor pressure at temperature T MM_{T}j average molecular weight of the liquid mixture
MM_{T},g average molecular weight of the gas mixture
MM, molecular weight of the component ith
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
* = ω ^{2}r + (1.1) dt p dz^{2}
^{"}_ 1 a f δT a^{2}τ ^{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 liquidvapour 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 gasliquid interface:
concentration axialsymmetry 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_{2}, 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 waxblock 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^{2}/S Dc3H70H, air= 1.02641 x 10^{'5} m^{2}/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^{"} (Tamb280)) 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^{5} 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^{4}. 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_{2}, ..., 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://wwwusers.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 ith 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 rolloff (a measure of shape) was studied by defining sensitivity coefficients as follows:
dlnH
^{3}thkJ lnx_{j}
_ d\n(roll off)
' rolloff, 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 singleside polishing machine. Preferably, the shape of the wafer W is measured and data from such measurement, such as the rolloff 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 207208, the processor calculates the pointtopoint difference between the two sets of data and generates a threedimensional 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 redistribute the threedimensional 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 printout 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_{I}OMEGAL,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(AH,0Z)
EXTERNAL RESJAC.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)_{1}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 algebraicdifferential equations DELTA(i)=G(T,Y/YPRIME)
C NEQ=IPAR
C
IMPLICIT REAL(AH,0Z)
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=(J1)+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+(NTY2)
DELTA(IEQ)=D1VZR 11 CONTINUE
NDELMAX=NDELMAX+2*(NTY2)
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+J1
DELTA(IEQ)=D1CR
IEQ=NDELMAX+NTYM2+J1
DELTA(IEQ)=DlCISOR
ENDDO NDELMAX=NDELMAX÷2*NTYM2
C B.C. D1CZ=BI0T*RH0GAS/CRIF*(FMH20GFMH20B) 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 FMISOLCISO(I,NTY)*CRTF/RHO
1 XLIQI=FMH20L/PMH20/(FMH20L/PMH20+FMISOL/PMISO+(l FMH20LFMISOL)
1 /PMSEC)
XISOI=FMISOL/PMISO/(FMH20L/PMH20+FMISOL/PMISO+(lFMH20LFMISOL) 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+(lYVAPIYISOI) 1 *PMARIA)
FMISOG=YISOI*PMISO/(YNAPI*PMH20+YISOI*PMISO+(lYVAPIYISOI) 20 1 *PMARIA)
C
ΓEQ=ΝDELMAX+I
BIOT=XKC/DLFH20*HMAX*H(I)
DELTA(IEQ)=D 1 CZBI0T*RH0GAS/CRIF*(FMH20GFMH20B)
25 C
IEQ=NDELMAX+NTX+I
BIOTISO=XKCISO/DIFISO*HMAX*H(I)
DELTA(IEQ)=D1CISOZBIOTISO*RHOGAS/CRTF*(FMISOG0.)
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=(TSATIT0(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 1TZQ*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+J1
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 1TZBI0TT*(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+(J2)*NTXM1+I1
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 — NAVIERSTOKES 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+(J2)*(NTX1)+(I1)
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+(lFMH20LFMISOL) 20 1 PMSEC)
XISOI=FMISOL/PMIS0/(FMH2OL/PMH2O+FMISOL/PMIS0+(lFMH2OLFMISOL) 1 PMSEC)
TI=T(I,NTY)*TRLF
YVAPI=P0_H2O(TI)/l.*XLIQI*GAMMA(XLIQI) ^{25} YISOI=P0_ISO(TI)/l.*XISOI
FMH20G=YVAPI*PMH20/(YVAPI*PMH20+YISOI*PMISO+(lYVAPIYISOI) 1 *PMARIA)
FMISOG=YISOI*PMISO/(YVAPI*PMH20+YISOI*PMISO+(lYVAPIYISOI) 1 *PMARIA) 30 C
IEQ=NDELMAX+(I1) 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.*(FMH2OGFMH2OB)+XKCISO/789.*(FMISOG0.)) _{ή 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+(J2)*(NTX1)+(I1) 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+(J2)*(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+(J2)*(NTX1)+(I1) 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(AH,0Z) DIMENSION Y(NEQ),YPRIME(NEQ)
INCLUDE 'spincom.f c
DO 1 JX=1,NTX D0 2 JY=1,NTY JS=(JX1)*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+(JX1)*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+(JX1)*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+(JX1)*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=(JX1)*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+(JX1)*NTY+JY YPRLME(JS)=CONC 1 (JX, JY) 41 CONTINUE
40 CONTINUE D0 45 JX=1_{;}NTX
D0 46 JY=1,NTY IS=NN4+(JX1)*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+(JX1)*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=(JX1)*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+(JX1)*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+(JX1)*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+(JX1)*NTY+JY
T(JX,JY)=Y(JS) CONTINUE CONTINUE ENDIF TIME DERIVATIVES
D0 6 JX=1,NTX DO 7 JY=l,NTY JS=(JX1)*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+(JX1)*NTY+JY
CONC 1 (JX, JY)=YPRIME(JS) 25 CONTINUE
20 CONTINUE
DO 45 JX=1,NTX DO 46 JY=1,NTY
JS=NN4+(JX1)*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+(JX1)*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 VRZD 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+(lFMH20LFMISOL) 1 /PMSEC)
XISOI=FMISOL/PMISO/(FMH20L/PMH20+FMISOL/PMISO+(1FMH20LFMISOL) 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+(lYVAPIYISOI) 1 *PMARIA)
FMISOG=YISOI*PMISO/(YVAPI*PMH20+YISOI*PMISO+(1YVAPIYISOI) 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.*(FMH2OGFMH2OB)+XKCTSO/J89.*(FMISOG0.))
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+(J2)*(NTX1)+(I1)
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+(J2)*(NTX1 )+(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(AH,0Z) 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 280300 K
XMUAΓR=(0.175+5.E4*(TAMB280))* 1.E4 !Pa*s
XKAIR=(0.0247+8.E5*(TAMB280))* l.E3 !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.15957E7 !M^{Λ}2/S DISALR=0.01013*TAMB**lJ5*4J48938E8
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(AH,0Z) INCLUDE 'spincom.f DTΓME=TFΓN/NINTER c c
NTX=NTX+N0X+N1X NTY=NJY+N0Y+N1Y
NTXM1=NTX1
NTYM1=NTY1
NTXM2=NTX2
NTYM2=NTY2 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(AH,0Z)
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)*TRIF273.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.E3
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. CHIMICAFISICA 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(AH,0Z) 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
_{30} 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.E20 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(NDNTY,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)_{I}ROOT(ND)
C AB=AL+BE
AD=BEAL 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=I1 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/(ZI) CONTINUE X=0. L=l
DO 20 1=1, N XD=0.
XN=1. XD1=0. XN1=0.
DO 30 J=1,N XP=(XDIF 1 (J)) *XNDJJ2(J)*XD XP 1 =(XDTF 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=ZCZ/(XROOT(Jl)) Z=Z/ZC X=XZ
LF (L.GT.500) GO TO 1 L=L+1
IF (ABS(Z/X).GT.l.E5) 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+1I 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=XROOT(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)_{5}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(AH,0Z) 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(ID2) 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(AH,0Z) DIMENSION Y(NEQ),YPRIME(NEQ) RETURN
END
C c
SUBROUTINE SDHELP(IDID,IFIL) IMPLICIT REAL(AH,0Z)
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(AH,0Z)
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(AH,0Z) 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(*,*)'illconditioned jacobian, restart' IF(IDID.EQ.8) WRITE(*,*)PD is singular, redundant equations^{1} IF(IDID.EQ.9) WRITE(*,*)'iHposed 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(AH,0Z)
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=(J1)+NDELMAX IVAR=J PD(IEQ,IVAR)=1. IEQ=(J l)+NDELMAX+(NTY2) DO N=l,NTX
IN R=NTX*NTY+NTX (N1)*NTY+J PD(IEQ,IVAR)=AL1X(1,N)
ENDDO 14 CONTINUE
NDELMAX=NDELMAX+2*(NTY2) C EQUAZIONE SULL'ANGOLO R=0 Z=H IEQ=NDELMAX+1
DO N=l,NTX
IVAR=NTX*NTY+NTX+(N1)*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+(I1)*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+(J1) DO N=l,NTX IVAR=J+(N1)*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÷J1 DO N=l,NTX
IVAR=2*NTX*NTY+NTX+(N1)*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+(J2)*NTXM1 +1 1 C DER. WITH RESPECT TO VR(I,J) rVAR=J+(Il)*NTY
PD(IEQ,IVAR)=VRRIF RMAX*(l R+ALlX(I,r))
DO N=l,NTX
LF (N.NE.I) THEN IVAR=J+(N1)*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 NAVIERSTOKES EQUATION
DO 24 I=2,NTXM1 DO 26 J=2,NTYM1
IEQ=NDELMAX+( J2)*(NTX2)+(I 1 ) ΓVAR=(II)*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=(I1)*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=(I1)*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+(II)*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.30363816.44/(TEMP46.13)))/760. END
FUNCTION PO SO(TEM) ! PO IN ATM C 1=76.964 C2=7623.8 C3=7.4924
C4=5.9436E18
P0_ISO=EXP(C 1 C2/TEMC3 *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/(TEMP273.15))*l.E3 END FUNCTION FERMI(X) !PROFILO INIZIALE
INCLUDE 'SPINCOM.F
FERMI=l.E3+HFERMI/(EXP((XAFERMI) BFERMI)+l) END FUNCTION OMEGA(TIME)
INCLUDE 'SPLNCOM.F' C IF (TIME.LT.TSPLN) THEN
IF (TIME.LT.TACC) THEN
OMEGAOME GAI+(OME GAFOMEGAI)/TAC C *TLME ELSE OMEGA=OMEGAF ENDIF C ELSE
C OMEGA=OMEGAF*EXP(0.5*(TLMETSPLNTACC))+0.001
C OMEGA=0.001
C ENDIF
END c
FUNCTION TSAT(TLME,NI) INCLUDE 'SPLNCOM.F
TSAT=(TSATIT0(NI))*TLME/TLMETERM+T0(NI) END
C
FUNCTION GAMMA(XCOMP) C
A=0.299753
B=13.224378
GAMMA=EXP((lXCOMP)**2*(A+2*(BA)*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 DifferentialAlgebraic Systems, SIAM J. Sci. Comp., 15 (1994), pp. 14671488
P.N. Brown, AC. Hindmarsh, and L.R. Petzold, Consistent Initial Condition Calculation for DifferentialAlgebraic Systems, LLNL Report UCRLJC122175, 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 (19961997)
Ohara T., Matsumoto Y., Ohashi H., Phys. Fluids 1, 1949, (1989)
Perry R.H., Green D., Perry's Chemical Engineering Handbook  50^{th} Ed., McGrawHill, (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, McGrawHill, (1988)
Schlichting H., BoundaryLayer Theory, McGrawHill, (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, PrenticeHall, (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
Priority Applications (1)
Application Number  Priority Date  Filing Date  Title 

PCT/IT2001/000094 WO2002069391A1 (en)  20010227  20010227  Method of forming a thin wax layer on a mounting block for mounting a semiconductor wafer 
Applications Claiming Priority (2)
Application Number  Priority Date  Filing Date  Title 

PCT/IT2001/000094 WO2002069391A1 (en)  20010227  20010227  Method of forming a thin wax layer on a mounting block for mounting a semiconductor wafer 
TW90110001A TW492062B (en)  20010227  20010426  Method of forming a thin wax layer on a mounting block for mounting a semiconductor wafer 
Publications (1)
Publication Number  Publication Date 

WO2002069391A1 true true WO2002069391A1 (en)  20020906 
Family
ID=11133626
Family Applications (1)
Application Number  Title  Priority Date  Filing Date 

PCT/IT2001/000094 WO2002069391A1 (en)  20010227  20010227  Method of forming a thin wax layer on a mounting block for mounting a semiconductor wafer 
Country Status (1)
Country  Link 

WO (1)  WO2002069391A1 (en) 
Citations (2)
Publication number  Priority date  Publication date  Assignee  Title 

US3475867A (en) *  19661220  19691104  Monsanto Co  Processing of semiconductor wafers 
US5256599A (en) *  19920601  19931026  Motorola, Inc.  Semiconductor wafer wax mounting and thinning process 
Patent Citations (2)
Publication number  Priority date  Publication date  Assignee  Title 

US3475867A (en) *  19661220  19691104  Monsanto Co  Processing of semiconductor wafers 
US5256599A (en) *  19920601  19931026  Motorola, Inc.  Semiconductor wafer wax mounting and thinning process 
NonPatent Citations (1)
Title 

GIVENS F. L., DAUGHTON W. J.: "ON THE UNIFORMITY OF THIN FILMS: A NEW TECHNIQUE APPLIED TO POLYIMIDES", JOURNAL OF THE ELECTROCHEMICAL SOCIETY., vol. 126, no. 2, February 1979 (19790201), ELECTROCHEMICAL SOCIETY. MANCHESTER, NEW HAMPSHIRE., US, pages 269  272, XP001031006, ISSN: 00134651 * 
Similar Documents
Publication  Publication Date  Title 

Salbu  Compressible squeeze films and squeeze bearings  
Wilkes et al.  Computational and experimental analysis of dynamics of drop formation  
Meunier et al.  Analysis and treatment of errors due to high velocity gradients in particle image velocimetry  
Yarin et al.  Evaporation of acoustically levitated droplets  
Cazabat et al.  Evaporation of macroscopic sessile droplets  
Koşar et al.  Thermalhydraulic performance of MEMSbased pin fin heat sink  
Hoffman  A study of the advancing interface. I. Interface shape in liquid—gas systems  
Davidchack et al.  The anisotropic hardsphere crystalmelt interfacial free energy from fluctuations  
Wildman et al.  Granular temperature profiles in threedimensional vibrofluidized granular beds  
US7082371B2 (en)  Fundamental mistuning model for determining system properties and predicting vibratory response of bladed disks  
Mackaplow et al.  A numerical study of the rheological properties of suspensions of rigid, nonBrownian fibres  
Ngan  On the nature of the dynamic contact angle: an experimental study  
Roy et al.  A lubrication model of coating flows over a curved substrate in space  
US6331075B1 (en)  Device and method for measuring thermal conductivity of thin films  
Cain et al.  Dynamic contact angles on smooth and rough surfaces  
Oligschleger et al.  Simulation of thermal conductivity and heat transport in solids  
Shen et al.  Stresses, curvatures, and shape changes arising from patterned lines on silicon wafers  
Hu et al.  Investigation of the natural convection boundary condition in microfabricated structures  
Escher et al.  Experimental investigation of an ultrathin manifold microchannel heat sink for liquidcooled chips  
Fouras et al.  Targetfree stereo PIV: a novel technique with inherent error estimation and improved accuracy  
Marsh et al.  Dynamic contact angles and hydrodynamics near a moving contact line  
Wantke et al.  Measurements of the surface elasticity in medium frequency range using the oscillating bubble method  
Charonko et al.  Assessment of pressure field calculations from particle image velocimetry measurements  
Francois et al.  Computations of drop dynamics with the immersed boundary method, part 1: numerical algorithm and buoyancyinduced effect  
Hodges et al.  The motion of a viscous drop through a cylindrical tube 
Legal Events
Date  Code  Title  Description 

AK  Designated states 
Kind code of ref document: A1 Designated state(s): CN JP KR SG US 

AL  Designated countries for regional patents 
Kind code of ref document: A1 Designated state(s): AT BE CH CY DE DK ES FI FR GB GR IE IT LU MC NL PT SE TR 

121  Ep: the epo has been informed by wipo that ep was designated in this application  
DFPE  Request for preliminary examination filed prior to expiration of 19th month from priority date (pct application filed before 20040101)  
122  Ep: pct application nonentry in european phase  
NENP  Nonentry into the national phase in: 
Ref country code: JP 