Sign up
Title:
COMPUTERIZED PROCESS CONTROL SYSTEM FOR THE GROWTH OF SYNTHETIC QUARTZ CRYSTALS
United States Patent 3805044
Abstract:
A computer-controlled, hydrothermal, crystal-growing process, for example, for growing synthetic quartz crystals of uniformly high Q. The disclosure includes a mathematical model of the crystal-growing autoclave and an algorithm by which a digital, process-control computer can make adaptive changes to the currents supplied to heating elements associated with the upper and lower zones of the autoclave, during the 25 day growth cycle, to yield a superior crystal of high Q.


Inventors:
Bhattacharyya, Ranendra K. (Kendall Park, NJ)
Fiore Jr., Angelo Ralph (Westville, NH)
Rudd, David William (E. Derry, NH)
Application Number:
05/131916
Publication Date:
04/16/1974
Filing Date:
04/07/1971
Assignee:
Western Electric Company Incorporated (New York, NY)
Primary Class:
Other Classes:
23/301, 117/71, 117/72, 117/202, 117/224, 117/926, 117/943, 700/29, 700/89, 700/117, 700/209
International Classes:
B01D9/02; (IPC1-7): G06F15/46
Field of Search:
235/151
View Patent Images:
Primary Examiner:
Botz, Eugene G.
Attorney, Agent or Firm:
Kain, Miller Winter W. M. R. P. R. C.
Claims:
1. A hydrothermal method for growing synthetic crystals which exhibit predetermined mechanical properties, comprising the steps of:

2. The method according to claim 1 wherein said heat applying step comprises:

3. The method according to claim 2 wherein said selective temperature altering step comprises:

4. The method according to claim 1 wherein said synthetic crystal comprises

5. A hydrothermal method for growing synthetic crystals which exhibit substantially uniform mechanical properties throughout all crystallographic regions of the crystal, comprising the steps of:

6. The method according to claim 5 wherein said heat applying step comprises:

7. The method according to claim 6 wherein said selective temperature differential increasing step comprises:

8. The method according to claim 5 wherein said crystal comprises quartz.

9. In a hydrothermal crystal-growing process of the type wherein synthetic crystals are grown in an autoclave, a method of selectively controlling the rate of crystal growth comprising the steps of:

10. The process according to claim 9 including the further steps of:

11. The process according to claim 10 wherein said selected physical parameters include at least the temperature of the upper and lower zones

12. A method of controlling the rate of growth of a synthetic crystal grown from a seed crystal suspended in a nutrient solution within a substantially vertically disposed autoclave, comprising the step of:

13. A method of maintaining substantially constant the rate of growth of a quartz crystal which is grown from a seed crystal suspended in a nutrient solution of solvated quartz within a substantially vertically disposed autoclave, comprising the step of:

14. Apparatus for hydrothermally growing synthetic crystals comprising:

15. Apparatus according to claim 14 wherein said computer means is programmed to include a mathematical model descriptive of the thermodynamic behavior of said autoclave, and said apparatus further comprises:

16. Apparatus according to claim 14 wherein:

17. Apparatus according to claim 14 wherein:

18. Apparatus for hydrothermally growing synthetic quartz crystals exhibiting substantially uniform mechanical Q throughout all crystallographic regions thereof, which comprises:

19. Apparatus according to claim 18 wherein said digital computer includes a memory storage device priorly programmed to include a mathematical model descriptive of the thermodynamic behavior of said autoclave to step changes in the heater currents supplied to the heating elements associated with said upper and lower zones, said memory storage device having further been priorly programmed to include the temperature profile that it is desired said upper and lower autoclave zones should attain, as time proceeds, said apparatus further comprising:

20. Apparatus according to claim 19 wherein said digital computer includes at least a central processing unit, a clock circuit, a control circuit and a process control input circuit and a process control output circuit, said apparatus further comprising:

21. The apparatus according to claim 26 wherein said current controlling means comprises first and second magnetic amplifiers respectively associated with the heating elements associated with the upper and lower

22. The apparatus according to claim 20 wherein said current controlling means comprises first and second silicon control rectifier circuits respectively associated with the heating elements associated with the

23. The apparatus according to claim 20, further comprising:

24. A method of controlling, by computer, a hydrothermal crystal-growing process, of the type wherein synthetic crystals are grown in an autoclave functionally divided into upper and lower zones and having electrically operated heating elements associated therewith, said crystal-growing process being divided into at least a warm-up phase and a run phase, said warm-up phase being subdivided into at least three regions, the method comprising the steps of:

25. A method of controlling, by computer, a hydrothermal crystal growing process, of the type wherein synthetic crystals are grown in an autoclave functionally divided into upper and lower zones and having electrically operated heating elements associated therewith, said crystal-growing process being divided into at least a warm-up phase and a run phase, the method comprising the steps of:

Description:
BACKGROUND OF THE INVENTION

1. Field of the Invention

Broadly speaking, this invention relates to the growth of synthetic crystals. More particularly, in a preferred embodiment, this invention relates to a computer-controlled process for growing synthetic crystals having substantially uniform mechanical properties throughout all crystallographic regions of the crystal.

2. Discussion of the Prior Art

In 1880, Pierre and Jacques Curie discovered that when mechanical stress was applied to opposite faces of a natural crystal, such as Rochelle Salt, an electric potential was developed across some other pair of faces. They also discovered that the application of an electric potential to opposite faces of the crystal resulted in mechanical strain, and a change in dimensions, between two other opposite faces of the crystal.

This phenomenon is called piezoelectricity and was regarded as little more than a scientific curiosity until the 1930's, when it was realized that the inherent mechanical resonance of a piezoelectric crystal could be utilized to control the frequency of an electronic oscillator, or to improve the band-pass characteristics of a tuned audio-frequency, or radio-frequency filter.

Of the many known crystals exhibiting this phenomenon, silicon dioxide, or quartz, has proved to be most suitable for use in electronic circuitry. This is due to the superior physical properties possessed by quartz, for example, high mechanical strength and good stability at temperatures in excess of 100°C.

An important factor in analyzing the performance of a vibrating piezoelectric crystal is its internal friction. This parameter may be defined as the ratio of the energy which is converted into heat within the crystal to the total energy supplied to the crystal when resonating at its natural frequency. Clearly, the internal friction is related to imperfections in the crystal structure. Because this energy loss is low, the numerical value for the internal friction of most crystals is a very small number. It is, thus, more convenient to define a quality factor, Q, equal to the reciprocal of the internal friction. This mechanical Q is analogous to the Q of a tuned L,R,C circuit defined as

Q = (ωL/R) = (1/ωCR)

where ω = 2π × the frequency of resonance.

The excellent circuit properties of quartz are due, in part, to the very high ratio of mass to elastance (this is equivalent to a high L/C ratio in a conventional tuned circuit) and to the very high ratio of mass to damping (this is equivalent to a high Q in a conventional tuned circuits). The Q of a natural quartz crystal is in the order of from 10,000 to 30,000 with values ranging up to 500,000 or more for specially treated crystals mounted in a vacuum.

The telephone industry is the largest single user of quartz crystals. These crystals are employed, for example, in the channel filters associated with frequency multiplex carrier systems, and as frequency control elements in the oscillators associated with coaxial cable carrier systems, and the like.

As the demand for telecommunications channels continues to increase, the telephone industry has been forced to use carrier systems having carrier frequencies and sidebands which fall higher and higher in the frequency spectrum. Over the last few decades, there has, thus, been a corresponding increase in the Q requirements for crystal control elements. For example, the Q requirements of a typical quartz resonator have been found to be as follows: from 0 to 100 Khz a Q of about 100,000 is required; from 100 to 300 Khz a Q of about 300,000 is required; from 300 Khz to about 1 Mhz a Q of from 400,000 to 800,000 is required; and from 1 to 10 Mhz crystals having Q's up to 1,000,000 are required.

Silicon dioxide, or natural quartz, is believed to constitute approximately 1/10 of the earth's crust. Although large deposits of low quality quartz are found in the United States and in Madagascar, this low quality quartz is of no direct use for the critical requirements of the electronic and telecommunication industries. Large crystals, of a quality sufficiently high for telecommunications use, are found only in Brazil. Although, in processed form, quite literally worth more than their weight in gold, Brazilian quartz is obtained, not from any organized mining operation, but rather from entrepreneurs who dig the crystals from the ground or find them in the beds of rivers. Obviously, this is a very unstable source of supply and requires the maintenance of large and expensive stockpiles. In any event, the supply of natural quartz crystals, of appropriate size and quality, is no longer capable of fulfilling the demands of the electronics and telecommunication industries.

Attempts to grow synthetic piezoelectric crystals date back well over 100 years, but because piezoelectricity itself was little more than a laboratory phenomenon, and there was no economic incentive, these attempts were generally unsuccessful. It was not until World War II that any serious effort was undertaken to find a way to create synthetic piezoelectric crystals of useful size and Q.

In the United States, the use of the synthetic piezoelectric material EDT (ethylene diamine tartrate) was considered, but artificial crystals of EDT did not prove entirely satisfactory as the very solubility of the crystals contributing to the ease with which they could be artificially grown, was a handicap during the processing of the crystals, as was their fragility. Also, EDT crystals were electronically inferior to natural quartz.

After World War II, a review of captured documents revealed that German scientists had achieved some success in growing synthetic quartz crystals using a hydrothermal technique with silica as the nutrient. The results of these experiments, however, were unpredictable and continuity was not achieved; only a porous mass of fine quartz needles resulted and these fine needles were, of course, totally unsuited for any commercial electronic application.

Spurred on by the German research, and foreseeing a shortage of natural quartz in the years to come, scientists in the United States continued research into the problem of growing synthetic piezoelectric crystals, of reasonable size and Q, which would be suitable for use in electronic applications.

United States Pat. No. 2,785,058, which issued on Mar. 12, 1957 in the name of E. Buehler, and which is assigned to a subsidiary of the assignee of the instant invention, disclosed the first really practical method of growing synthetic quartz crystals suitable for electronic use. As disclosed in that patent, a long cylindrical autoclave is divided by a baffle plate into an upper zone, called the crystallization zone, and a lower zone, called the dissolving zone. The lower, or dissolving zone, comprises about one half the total volume of the autoclave and is filled with small pieces of natural quartz, referred to as the nutrient. In the upper, or crystallization zone, small plates of quartz, of known orientation, called the seed crystals, are suspended from a rack. The remaining volume of the autoclave is filled to about 80 percent capacity with a solvent which is capable of dissolving quartz, for example, a week aqueous, alkaline solution of sodium hydroxide (NaOH). As discussed in U.S. Pat. No. 3,356,463, dopants, such as lithium nitrate (LiNO2) may be added to the solution to improve the characteristics of the crystals. The vessel is then sealed and heat is applied to the upper and lower zones of the autoclave. Typically, the upper zone is heated to about 630°F and the lower zone to about 700°F. The action of the heat causes the nutrient quartz to dissolve in the alkaline solution and, because the upper zone is relatively cooler than the lower zone, the solvated quartz flows by convection to the upper end of the autoclave. The baffle plate which separates the two halves of the autoclave contains a plurality of small holes, and serves to maintain two essentially isothermal zones within the vessel. In addition, the baffle plate helps to channel the convected flow of solvated quartz from the lower zone to the upper zone. The action of the heat which is applied to the autoclave causes the aqueous solution of sodium hydroxide within the autoclave to expand, generating pressures in the order of 25,000 pounds per square inch. Because the upper, or crystallization zone, is some 70°F cooler than the lower zone, when the aqueous solution arrives at the upper zone, it cools and quartz precipitates out of solution and nucleates onto the seed crystals in the form of a single crystal. The process varies in time, depending upon the size of the crystals desired, but nominal growing times are in the order of from three to four weeks.

It has been observed that the Q of crystals grown by the process described in U. S. Pat. No. 2,785,058, varies throughout all crystallographic regions of the crystal. This has been attributed to impurities in the crystal which are segregated in response to variations in growth rate during the growth process. These variations in growth rate are believed to be a consequence of the design of the autoclave and its relation to the growing crystal, in addition to the natural anisotropy of the crystal. For example, during the first few days of synthesis, the size of the crystal is not sufficiently large to inhibit the convective flow of solvated quartz from the lower zone of the autoclave to the upper zone of the autoclave. The crystals, therefore, grow at a rapid rate. This effect, however, decreases with time, because as the crystals grow, they occupy an increasingly larger portion of the upper autoclave zone, thereby inhibiting the convective flow of solvated quartz, and reducing the rate of crystal growth. The decrease in growth rate leads, in turn, to a rejection of impurities, with a subsequent enhancement of the crystal Q. Thus, the Q of the crystal in the vicinity of the seed is lower than the Q at the extremities of the crystal.

In the above-described process, the seed crystals grow in X, Y, and Z directions. However, because of the crystallographic orientation of the seed crystal, and the physical arrangement within the autoclave, by far the greatest amount of growth occurs in the Z direction. For this reason, the above-described crystal growing process is referred to as a fast Z-growth, or fast basal growth process.

At the end of the growth cycle, the synthetic crystals are removed from the autoclave for further processing. Crystal units are then cut from the synthetic crystals and, as is well known, these cuts may be made in any of several orientations, for example, X cut, Y cut, BT cut, CT cut, etc. For telecommunication purposes, however, most quartz crystal devices operate in the lower and middle frequency ranges. Thus, the majority of crystal cuts are made in the basal plane, parallel to the seed crystal. Unfortunately, the Z growth plane is particularly susceptible to impurity segregation and, as previously discussed, these impurities tend to decrease the value of the crystal Q by increasing the internal friction of the crystal. Thus, crystal units which are cut from slices parallel to the seed crystal are more likely to have larger Q variations than crystal units which are cut from slices making an angle with the seed.

As previously mentioned, because the rate of crystal growth falls continuously during the growth cycle, the number of impurities which succeed in attaching themselves to the crystal also falls, thus, steadily raising the Q of the crystal in the direction of the crystal surfaces. Accordingly, the Q of a basal cut slice will vary throughout the slice and the average Q, from slice to slice, will also vary.

The actual difference in Q between a basal cut slice taken close to the seed and one taken close to the Z surface is considerable. For example, a slice taken close to the seed may have a Q of only 120,000, whereas a slice cut from near the outer Z surface of the same crystal may have a Q of as high as 260,000, more than a 2:1 ratio. Clearly, this is a most undesirable situation, for if the Q of the near slice is considered to be the limiting factor, then slices taken further from the seed will have a Q far in excess of the requirements of the circuit with which they are intended to operate, resulting in wasteful economic loss. On the other hand, if the Q of slices taken some distance from the seed crystal is considered to be the limiting factor, then the entire portion of the crystal which is close to the seed must be scrapped, resulting again, in wasteful economic loss.

The problem, then, is to find a technique for growing synthetic crystals whereby the crystals may be made to exhibit predetermined mechanical properties throughout selected crystallographic regions thereof.

This problem has been solved by the instant invention, in which, a substantially vertically disposed autoclave is charged with a crystal growing nutrient, the autoclave being functionally divided into two essentially isothermal zones. The autoclave is then partially filled with a solvent capable of dissolving said nutrient, and at least one seed crystal suspended, proximate the upper zone of the autoclave. The autoclave then is sealed and heat is applied to the upper and lower zones thereof to increase the pressure and temperature above the levels required to dissolve the nutrient in the solvent. The temperature of the lower zone is established at a higher value than that of the upper zone to permit a convective flow of solvated nutrient over the surfaces of the seed crystal. Thereafter, the temperature differential between the upper and lower zones is altered selectively as time proceeds so as to change the rate of flow of solvated nutrient over the seed crystal. In this fashion, the growth rate of the crystal is selectively altered thereby causing the crystal to exhibit the desired predetermined mechanical properties.

One illustrative apparatus for practicing the above method comprises a substantially vertically disposed autoclave which is functionally divided into an upper and a lower zone. The autoclave is adapted to receive a charge of crystal growing nutrient in the lower zone and of being substantially filled with a solvent capable of dissolving the nutrient. The apparatus further includes means for suspending at least one seed crystal within the autoclave, proximate the upper zone thereof as well as first and second heating means which are respectively associated with the upper and lower zones.

A source of energy for the first and second heating means is also provided as well as computer means for respectively adjusting the amount of energy supplied to the first and second heating means from said source, whereby the temperature differential between the upper and lower zones of the autoclave is selectively altered during the crystal growing process.

The invention, and its mode of operation, will be more fully understood by reference to the following detailed description and to the drawings, in which:

DESCRIPTION OF THE DRAWINGS

FIG. 1 is a partially cut-away, isometric view of an hydrothermal autoclave of a type suitable for use with the present invention;

FIG. 2 is a graph depicting the relationship between crystal growth and growth rate, as a function of time, for a crystal grown in the autoclave of FIG. 1;

FIG. 3 is a graph depicting the relationship between the Q of a crystal and the elapsed time for a crystal grown in the autoclave of FIG. 1;

FIG. 4 is a drawing, prepared from an actual photography of a shadowgram, showing how the Q of a crystal is assigned to different physical locations within the body of the crystal;

FIG. 5 is a graph depicting the relationship between the Q of a crystal and its growth rate for a crystal grown in the autoclave of FIG. 1;

FIG. 6 depicts an illustrative profile of the desired temperatures for the upper and lower zones of a computer-controlled autoclave according to the present invention;

FIG. 7 is a graph which depicts the autocorrelation between samples in the mathematical model of the autoclave of the present invention;

FIG. 8 is a graph depicting the root mean square error, as a function of time, for the mathematical model of the autoclave of the present invention;

FIG. 9 depicts a graph of the estimated value of the parameter w, as a function of γ, for region 4 of the upper zone in the mathematical model of the autoclave of the present invention;

FIG. 10 is a graph depicting the relationship between the upper and lower zone temperatures, as a function of time, during the warm-up phase of the autoclave of FIG. 1;

FIG. 11 depicts the formula used to predict the temperature of the autoclave of FIG. 1;

FIG. 12 depicts the formulae used to compute the heater currents for the heating elements associated with the autoclave of FIG. 1;

FIG. 13 depicts the formulae used to update the value of the parameter K in the mathematical model of the autoclave of FIG. 1;

FIG. 14 depicts the values of the estimated off-line parameters for the mathematical model of an experimental hydrothermal autoclave actually built and tested;

FIG. 15 depicts one illustrative apparatus configuration for practicing the present invention;

FIGS. 16 - 22 depict an illustrative logical flow chart for implementing the algorithm of the present invention;

FIG. 23 illustrates the manner in which FIGS. 16 - 22 should be assembled;

FIG. 24 depicts a graph of the estimated value of the parameter w, as a function of γ, for region 4 of the lower zone in the mathematical model of the autoclave of FIG. 1; and

FIG. 25 depicts a graph of the estimated value of the parameter w, as a function of γ, for region 1 of the upper zone in the mathematical model of the autoclave of FIG. 1.

DETAILED DESCRIPTION OF THE INVENTION

This invention is based upon the discovery that the uniformity of synthetically-grown crystals may be improved by controllably increasing the temperature differential between the upper and lower zones of an autoclave, throughout the growth cycle, to maintain a relatively constant rate of flow of solvated nutrient over the seed crystals.

The effect of increasing the temperature differential between the upper and lower zones of the autoclave, considered by itself, is to increase the rate of flow of solvated nutrient over the seed crystals. However, as previously noted, as the crystals expand they fill an ever-increasing proportion of the autoclave. Thus, there is a tendency for the rate of flow of solvated nutrient to steadily decrease with time. By carefully adjusting the temperature differential as a function of time, these two effects can be balanced, one against the other, so that the net effect is a relatively constant flow of solvated nutrient over the seed crystals, throughout the entire growth cycle.

Implementation of the above-described technique, however, is by no means a simple matter. First, the crystal growing process must be mathematically characterized so that the temperature differential necessary to maintain a constant rate of convective flow can be calculated for any instant of time during the growth cycle. Next, a real-time, direct digital control algorithm must be derived so that appropriate changes may be made to the current which is supplied to the heating coils associated with the upper and lower zones of the autoclave to produce the desired temperature differential for the time period of interest.

The illustrative embodiment of the invention, to be discussed below, pertains to the growth of synthetic quartz crystals. One skilled in the art, however, will appreciate that the invention is not so limited, and may be used with equal facility in any hydrothermal crystal growing process. The use of the invention to grow crystals other than quartz merely requires modification of the parameters of the control algorithm and the crystal growing equations.

Before discussing, in detail, the mathematical characterization of the crystal growing process and the control algorithm, the operating environment will briefly be reviewed.

Hydrothermal techniques for the growth of synthetic quartz are, of course, well known and have been widely reported in the literature (see, for example, the article entitled, "Growing Quartz Crystals," by R. A. Sullivan, Western Electric Engineer, Vol. III, No. 2, (April 1959), pp. 2-10). Thus, the autoclave to be described below, per se, forms no part of this invention, but is described merely to facilitate an understanding of the development of the process characterization equations and the control algorithm.

As shown in FIG. 1, crystal growing apparatus 100 comprises a hollow, cylindrical autoclave 101 vertically disposed on a bed of insulating fire brick 102 set in a concrete support pad 103. Typically, autoclave 101 is fabricated from a chromium-molybdenum steel alloy and measures 9 feet long, with an inside diameter of 6 inches and a wall thickness of 31/2 inches. A perforated baffle 104, positioned about halfway up the central aperture, divides the autoclave into an upper, or crystallization zone, and a lower, or dissolving zone. A wire basket 106, filled with crushed nutrient quartz 107 is positioned in the lower zone of the autoclave. A wire rack 108 supporting a plurality of seed crystals 109 is positioned in the upper zone of the autoclave above the crushed quartz. The remaining portion of the autoclave aperture is filled to about 80 percent capacity with a weak, aqueous, alkaline solution of sodium hydroxide (NaOH) and one or more dopants, such as Lithium Nitrate (LiNO2). A threaded seal 111, typically a modified Bridgeman seal, is screwed into the upper end of autoclave 101 to form a sealed enclosure. A first plurality of heater bands 112 are positioned about the upper zone of the autoclave and connected, in parallel, to a first source of electric energy (not shown). A second plurality of heater bands 113 are positioned about the lower zone of the autoclave and connected, in parallel, to a second source of electric energy (not shown). Typically, each heater band comprises two 440-volt, 500-watt heater elements, connected in parallel, and clamped end to end around the circumference of the autoclave. The current which is supplied to the heater bands associated with the upper and lower zones of the autoclave may be individually adjusted to any desired value. A pressure sensing device 114 passes down through an aperture in threaded seal 111 into the body of the autoclave to permit the pressure within the sealed portion to be monitored throughout the crystal growing process. A first thermocouple 116 is placed within a well formed in the wall of the autoclave to permit the temperature in the lower zone of the autoclave to be continuously monitored. A second thermocouple 117 is placed within a well formed in threaded seal 111, to similarly permit the temperature in the upper zone of the autoclave to be continuously monitored. One or more additional thermocouples 118 and 119 may be positioned in contact with the upper and lower heater bands, respectively, to permit the temperature at the outer surfaces of the autoclave to be monitored. A rupture disc 121 is connected to threaded seal 111 to minimize the risk of a catastrophic explosion if the pressure within the autoclave, typically from 25,000 to 30,000 pounds per square inch, accidentally rises above the safe limits of the autoclave. The entire autoclave is surrounded by a jacket 112 of thermal insulation, such 122 asbestos, to reduce heat loss.

In the operation of prior art apparatus, such as above-described, current is supplied to heater bands 113, which are associated with the lower zone of the autoclave, to rapidly increase the temperature thereof to a value of approximately 700°F. At the same time, current is supplied to heater bands 112, which are associated with the upper zone of the autoclave, to increase the temperature thereof to some value approximately 70°F less than the temperature of the lower zone of the autoclave. As a result of the application of heat to the autoclave, the aqueous solution of sodium hydroxide expands to completely fill the central aperture of the autoclave and then continues to expand raising the pressure within the autoclave to approximately 25,000 pounds per square inch. Under these conditions, the crushed nutrient quartz 104 in basket 106 dissolves in the sodium hydroxide solution and rises, by convection, to the upper half of the autoclave, which it will be recalled, is some 70°F cooler than the lower half. Because of the difference in temperature, the solvent cools and quartz precipates out of solution and deposits on seed crystals 109 steadily increasing their size, particularly in the Z direction.

In the prior art, appropriate heater current adjustments are made manually, or by means of an analog control device, to maintain the temperature of the upper and lower zones of the autoclave substantially constant or, in any event, to provide a constant temperature differential therebetween. At the end of approximately 25 days, when crystals 109 have grown to the desired size, the supply of current to the heater bands is terminated and the autoclave permitted to cool. Seal 111 is then broken and rack 108 with the fully grown crystals thereon removed from the autoclave, followed by basket 106 and what little nutrient remains in undissolved form. Crystals grown by the above process will have a Q factor which varies throughout all crystallographic regions of the crystal.

CHARACTERIZATION OF THE CRYSTAL GROWING PROCESS

As previously mentioned, before one can calculate the precise manner in which the temperature of the upper and lower zones of an autoclave must be varied, as a function of time, to grow crystals of uniform Q, it is first necessary to characterize the crystal growing process. An exact mathematical characterization of the crystal growing process is extremely difficult, due to the large number of process variables present. For example, among the parameters affecting the growth rate and the Q of a synthetically grown crystal, are the percentage of autoclave fill, the concentration and chemical composition of the dopants, the concentration and the chemical composition of the solvant, the amount and nature of the nutrient, the nature and crystallographic orientation of the seed crystals, etc., etc. Fortunately, because of experience gained during many years of crystal growth, we realized that the characterization process could be vastly simplified by choosing specific operating conditions known to produce acceptable crystals and building around these conditions process equations of utility, yet which are sufficiently simple to be easily programmed on a general purpose digital computer. The problem, then, is reduced to that of characterizing the growth process in terms of growth rate and Q, since these two variables are the ones of greatest interest, when considering the individual merit of a grown crystal. As a first step in the assembly of the data required to characterize the growth process, a curve was constructed which depicts the thickness of a synthetic crystal as a function of the growing time. This curve, FIG. 2 of the drawing, was constructed by performing a number of experimental syntheses under identical, known operating conditions. The individual runs were terminated at frequent intervals, for example, every three days, and the thickness of the resultant crystals was measured. Typical operating conditions for these experimental syntheses were temperatures of 675°F for the crystallization zone and 745°F for the dissolving zone, with about 83.5 percent fill. The measurements obtained by these series of experiments were then plotted and a curve was fitted to the plotted data, with the aid of a computer, using the method of least squares. Thus, at a specific, nominal crystallizing temperature and a specific temperature differential between the dissolving and the crystallizing zones, the following growth-time relationship was obtained: ##SPC1##

where G is the crystal growth in mils; t equals time in days; ai is an experimentally defined coefficient. The derivative of Equation 1, dG/dt, which expresses the rate of crystal growth as a function of time, is also illustrated in FIG. 2. It will be noted that Equation 1 is not defined for values of t less than 0.5 days. The reason for this is that, during the first half-day of the growth cycle, the seed crystals are etching because, at this time, the solvent is essentially pure sodium hydroxide (NaOH).

As the next step in the characterization process, the spectrophotometric technique disclosed in U. S. Pat. No. 3,351,757 was used to measure the Q at different thicknesses within each crystal grown by the experimental synthesis. The applicable growth equation was then used to relate the average Q's, so obtained, to the growing time. For the conditions illustrated in FIG. 2, a graph of the relationship between the average Q and the growing time is shown in FIG. 3. More specifically, a crystal which was grown under the conditions establishing Equation 1 was cut to produce a Y-cut optical flat about 0.25 inches thick. This flat was then scanned using infrared light at a wave number of 3500 cm-1. Next, a shadowgram of the Y-cut optical flat was prepared and matched to the infrared scan so that positions within the growth regions could be identified with the Q of the crystal. This analysis is shown in FIG. 4 Because the magnification effect of the camera is known, by utilizng certain landmarks in the optical flat (principally the location of the seed crystal) the infrared scan was divided into segments with numbers ascending in chronological order on either side of the seed crystal, which segments are proportional to thickness segments in the original crystal. Each of these thickness segments may be identified by a Q value expressed in terms of α, where α is defined as the extinction coefficient at room temperature using light with a wave number of 3500 cm-1. The value of α is obtained by noting the absorption difference between an infrared scan of the optical flat made with infrared light of a wave number of 3800 cm-1 and a second scan made with infrared light of a wave number of 3500 cm-1. Thus, FIG. 4 relates the Q description of a crystal grown under the conditions described by Equation 1 to various thickness segments in the growth region. With this analysis the Q, as a function of growth rate, throughout the crystal may be determined. This information is obtained by solving Equation 1 for the time, t, and d(thickness)/dt for each of the thickness segments shown in FIG. 4. After performing this operation, the Q for each thickness segment is identified with its growth rate so that the Q as a function of growth rate is obtained for a crystal grown under the operational conditions defined by Equation 1. These results are given in Table A, below, which lists the average growth rate, the final thickness of the crystal, the Q description throughout the crystal, the average Q and the range in Q for each region of the crystal. FIG. 5 of the drawing is a graph illustrating the data set forth in Table A.

TABLE A

Q VS GROWTH RATE; Tupper =675°F, Tlower =745°F, % --------------------------------------------------------------------------- FILL=83.5

Total Time Growth Region Growth in Rate Q×10-6 QRANGE Y 10- 6 mils/day __________________________________________________________________________ 1 122 1.46 130 0.357 0.357-0.356 2 244 2.44 119 0.432 0.469-0.401 3 366 3.51 109 0.452 0.411-0.403 4 488 4.68 99.3 0.517 0.574-0.461 5 610 5.97 90.2 0.583 0.639-0.527 6 730 7.39 81.7 0.612 0.659-0.564 7 854 8.96 74.0 0.645 0.721-0.569 8 976 10.70 67.5 0.706 0.797-0.616 9 1098 12.60 61.4 0.812 0.943-0.681 10 1220 14.70 55.9 0.853 0.958-0.748 11 1342 17.00 50.6 1.07 1.27-0.879 12 1464 19.50 44.5 2.03 2.09-1.98 13 1586 22.50 36.2 >2.2 >2.2 __________________________________________________________________________

From Table A, it will be observed that the average growth rate is 78 mils/day and that the Q at this average growth rate is 645,000. It will alwo be noted that the range in Q is from 350,000 to 2,200,000.

If the assumption is made that the rate of crystal growth is proportional to the rate at which solvated quartz flows over the seed crystals, then from the results of the experimental syntheses, illustrated in FIGS. 2, 3, and 5, it is a relatively easy matter to compute what the autoclave temperature differential should be at any given instant of time, so that the overall growth rate will remain relatively constant, thereby yielding a crystal of uniform Q.

When the growth rate is so adjusted, the Q in the vicinity of the seed will be elevated and the Q in the outer Z extremity will be diminished, thereby raising the average Q of the crystal and the uniformity of Q throughout the crustal. Thus, crystal plates may be cut from the entire crystal, and the economic waste of the prior art is avoided.

FIG. 6 illustrates the general shape of the temperature versus time curves for the upper and lower zones of a typical hydrothermal autoclave. Strictly speaking, the assumption that the growth rate is proportional to the rate at which the solvated quartz flows over the seed crystals is incorrect, as it ignores the fact that the geometry of the autoclave may affect the growth rate. Thus, the temperature versus time curves will be slightly different for different autoclave geometries. However, the curves shown in FIG. 6 are typical.

It will be noted that, initially, the temperature differential between the upper and lower autoclave zones is about 70°F, as in the prior art. As time progresses, however, it will be noted that the temperature differential is steadily increased; for example, to about 125°F by the tenth day, to about 170°F by the 15th day, and close to 200°F by the 25th day. For practical reasons it has been found to be more expedient to only slightly raise the temperature of the lower zone and to supply most of the increase in temperature differential by dropping the temperature of the upper zone. It will be apparent, however, that the increase in temperature differential could be aportioned between the temperatures of the upper and lower zones; for example, by holding the temperature of the upper zone constant and raising the temperature of the lower zone. It should be born in mind, however, that if the temperature of either zone is raised too high, the hydrothermal process may be affected and the limit of the physical strength of the autoclave may be reached.

An analytical formulation of all the physical phenomena that relate the current in the heater coils to the temperature at various points within the autoclave is extremely difficult. Among the factors to be considered are:

1. the convection of the viscous fluid within the autoclave;

2. the conduction of heat along the wall of the autoclave; and

3. the heat radiation between the autoclave and its surroundings.

Let the height of the autoclave be h and the radius be a. Since the autoclave is symmetrical about its vertical axis, a two-dimensional co-ordinate system can be erected where z denotes the direction of the axis, and r denotes the direction along a radius. Let the density of the fluid be denoted by ρ, its kinematic viscosity by ν, and its thermal diffusivity by λ. Let v and u be the vertical and the horizontal components of the convection velocity of a fluid particle. Now define the following dimensionless quantities:

τ = λt/h2 where t is time

Z = z/h, R = r/a, U = hu/λ, V = hv/λ

θ = normalized temperature

The energy equation for convection is given by the equation ##SPC2##

No closed-form analytic solution of Equation 2 is known to exist to satisfy arbitrary boundary conditions.

The heat conduction along the wall of the autoclave is essentially a diffusion phenomenon given by the equation

(δθ/δτ) = G∇2 θ (3)

where,

2 is the Laplacean operator, and G is a constant whose value depends on the thermal conductivity, the specific heat, and the density of the wall material.

The radiation phenomenon is governed by Boltzman's Law, according to which the rate of energy radiation from the autoclave is:

φ[θ4 - Ω] where Ω is a constant dependent on the room temperature (assumed to be constant), and φ is another constant. If a direct current of fixed magnitude x is fed into the autoclave heater coil, the rate at which energy is fed into the autoclave, after the coil has reached a steady temperature is ψx2, where ψ is a constant. Thus, the energy balance equation becomes:

ψx2 - φ[θ4 - Ω] + N = ηdθ/dt (4)

where,

η is a constant dependent on the specific heat and the mass of the autoclave, and N is a random term representing the stochastic nature of the process.

The over-all relationship between the autoclave heater current and the temperature at various points in the autoclave is a complex combination of Equations 2, 3, and 4. But the problem associated with developing a complicated model is that, first of all, no matter how complex the model is, it still remains an inexact structural description of the process. Secondly, with the increase in complexity of the model, the number of parameters to be estimated increases. Because the number of observation samples is limited, the estimation error increases, due to the lack of redundancy. Thirdly, the model becomes very sensitive to errors in estimation, as can be expected due to the presence of the second and the fourth degree terms in Equation 4. Finally, a complicated equation makes mathematical manipulations rather unwieldy.

From a series of experiments run on an actual autoclave, a study of current and temperature data reveals two factors that are pertinent to the temperature control of the autoclave. The first is that, for a step input of heater current, the temperature rise has the general shape of the response of a first order differential equation with a transport delay associated therewith. The second factor is that if the heater current is chosen to be a rectangular wave whose width is T2 - T1 (i.e., the rectangle exists between the two instants of time T1 and T2), and, if the resulting temperature starts rising at T1 + τ(where τ is the transport delay), the temperature does not necessarily start falling at T2 + τ. Depending on the magnitude of the rectangular current, the temperature may still go on rising at T2 + τ for a few minutes to come. This second factor is taken care of, in an approximate way, in the design of the control algorithm, to be discussed below. The first factor can be described by the following equation

dy(t + τ)/dt + Ay(t + τ) = Bx(t) + C + N (5)

where,

y(t + τ) -- is the temperature at the instant t + τ;

x(t) -- is the heater current at the instant t;

τ -- is a parameter representing the transport delay (i.e., the amount of time required for the effect of a sudden change in current to be detected at the thermocouple);

N -- is the stochastic noise term that constitutes the random element in the response of the autoclave temperature;

A, b, and C are the three parameters of the model. The parameter A determines the time constant of the system (i.e., the amount of time during which the temperature goes through a transient state until it attains a new steady-state temperature for a new setting of the current); B/A is the steady-state differential gain (i.e., the amount of change in the steady-state temperature for a unit change in current); and C/A is the zero-set of the temperature (i.e., the steady-state temperature of a zone of the autoclave when there is no current in the corresponding heater). This latter term takes care of the room temperature as well as the "temperature coupling" between the upper and the lower zones. The temperature coupling can be defined as the effect on the temperature of the upper zone due to the current in the heater of the lower zone, and vice-versa.

Note, that there is no valid mathematical reason to allow the two constants C and N to exist separately, because their summation can be replaced by a single variable. However, they are kept separate to explicitly express the different physical factors involved in the process. Equation 5 has proved to be a satisfactory model that approximates the process on a piece-wise basis, within an acceptable bound on structural errors. Since Equation 5 is a piece-wise approximation, the values of the four parameters of this equation are dependent on the operating point of the process.

Provisions may be made to measure sample values of x(t) and y(t). However, the measurements are not free from errors. Let the measured values of the samples xi and yi be denoted as ui and vi where,

ui = xi + ζi and vi = yi + ξi (6)

and ζi and ξi are the measurement errors. The current and temperature transducers are such that ζi is negligible, and ξi is within a band of ±1°F. Since ζi is negligible, ui is identically equal to xi for all practical purposes, and henceforth will be used synonimously.

Equation 5 can be written in a discrete form, that is,

vi+k+l - vi+k + A/2 [vi+k+l + vi+k ] Δt

= B ui Δt + C Δt + εi (7)

where, εi is a function of ζi, ξi+k, ξi+k+l and Ni. The interval between two consecutive samples is Δt, and k is an index such that kΔt = τ. Equation 7 can be rewritten as:

αi + Aβi + Bγi + K = εi (8)

where the definitions of αi, βi and γi are obvious by direct comparison with Equation 7. A, B and K are the parameters of the model to be estimated. The index k in Equation 7 is another term to be evaluated.

The model described by Equation 8 is only an approximate description of the heating process, with A, B, K and the index k (i.e., τ) as the unknowns, and with εi as the stochastic term. The question arises as to how many parameters shuld be estimated, how often the parameters should be updated, and what type of statistical estimation scheme should be followed. Since the process, in reality, is nonlinear, and Equation 8 is only a structural approximation therefor, the well-known properties of a least square estimate for a linear equation (i.e., efficiency, consistency and unbiasedness) cannot be claimed. Nevertheless, because of its simplicity and acceptable overall results, the method of least square estimates was used.

Some initial experiments were run with x(t) as a step function of different magnitudes, and v(t) was measured in the transient, as well as in the steady state. The transient regions reveal that τ can be considered to be a constant, independent of the values of x and y. Typically, the value of this constant is 12 minutes for the bottom zone of the autoclave and 15 minutes for the upper zone of the autoclave. The autocorrelation of the steady-state regions has the general shape shown in FIG. 7. One plausible interpretation of this figure is that consecutive samples of ε begin to be correlated as the interval between two consecutive samples becomes smaller than three minutes. Since samples of ε taken at intervals not less than three minutes are uncorrelated, no control action can be taken on x(t) to diminish ε(t), especially in view of the fact that the transport delay between x(t) and y(t) is of the order of 12 to 15 minutes. Hence, it is unnecessary to make estimates of ε, in the sense that this information does not serve any useful purpose. Thus, the parameters A, B, and K in Equation 8 may be estimated, while εi constitutes the residual error.

The steady-state value of v(t) for a certain step function x(t), would be constant, for all practical purposes, if ε(t) were zero and if A, B and K remained constant. The presence of ε(t), however, gives rise to high frequency fluctuations. It should also be noted that there are slowly varying trends in the steady-state region of v(t), proving thereby that at least one, if not all, of the parameters A, B, and K are slowly time-varying. To establish this slow time-variation, a certain observation interval T was chosen in the steady-state region of the temperature response for a step current input. A straight line, parallel to the abscissa was fitted through v(t) over the interval T, and the mean square error of the fit calculated. As shown in FIG. 8, this mean square error, plotted as a function of T, as T varies from 15 minutes to 5 hours, shows that for all practical purposes, A, B, and K can be considered to be constant over a period of approximately one and a half hours. In other words, the time-varying system can be "frozen" over a period of approximately one and a half hours. In a real-time, on-line, control situation, this "freeze" time must be divided into two parts; the observation interval, and the prediction interval. The observation interval is the time during which input-output data are collected to estimate the parameters of the model, and the prediction interval is the time over which the estimated model performance is extrapolated and assumed to be the predicted value of the output of the process, on the basis of which the control apparatus has to adjust itself. It is believed, in this particular case, that the one and a half hour "freeze" time can be effectively divided into a half-hour observation interval and a one hour prediction interval. A question arises as to which of the three parameters A, B, and K should be estimated on an on-line basis and which of them should be estimated on an off-line basis. A study of the transient regions of v(t) for x(t) as a step function reveals that an estimate of A obtained from data collected over half an hour is too sensitive to ζ, ξ and N. This sensitivity is due to the fact that half an hour is so short a time that the small segment of the exponential transient of v(t), over the half-hour interval, is almost a straight line. Hence, the presence of ζ, ξ and N introduces a significant error in the estimate of A. This estimation error, in turn, produces unacceptable prediction errors. These considerations have led to an off-line estimation scheme for A. Values of this estimate are obtained at different regions within the operating range and are stored in the computer controlling the crystal growing process. When the process is in operation, the computer selects that value of A which corresponds to the operating region of the process at that moment.

To determine whether both B and K should be estimated on an on-line real-time basis, a new parameter W is defined as follows:

W = Bγ + K (9)

since, as previously mentioned, ζi is negligible, ui will remain, more or less, constant if xi is held at a known constant value for all values of i. Thus, if the current is kept constant, the subscript i can be dropped from γi in Equation 8. The estimated value of W is plotted as a function of x in FIG. 9. The linear nature of this plot suggests that B can be reasonably approximated as a constant. Hence, it can be estimated on an off-line basis for different operating regions of the process, and the corresponding values can be fed to the computer when the process is in operation. Thus, the only parameter that is held responsible for the slow time variation of the process is K. The estimation scheme consists, therefore, of an off-line, least square estimation of A and B, from previously collectd data, for various regions of the process operation, and an on-line, real time least square estimation of K while the process is in operation.

Let A and B be the off-line estimates of A and B, respectively. Let Ko be the initial least square estimate of K that is fed to the computer, along with A and B, before enough on-line data are collected to update the estimate of K on a real-time basis. The value of Ko is also evaluated, on an off-line basis, along with A and B. The observation time (i.e., the estimation time) to estimate K is approximately half an hour, as mentioned before. The rate of sampling is determined by the fact that the consecutive samples of ε begin to be correlated if the sample interval is less than three minutes, as shown in FIG. 7. Accordingly, the sample interval was chosen to be three minutes. This allows ten samples per observation period. Since only one parameter, that is K, is estimated on an on-line basis, the ratio between the number of samples and the number of parameters to be estimated is sufficiently large to make the estimate statistically meaningful. The expression for K is given by ##SPC3##

where n is the number of samples collected during the half-hour observation period. The above equation is used every three minutes (i.e., the sample interval) to update K on a real-time basis. As mentioned before, A and B are off-time estimates. Depending on the operating point, the corresponding values of A and B are fed to the computer.

The purpose of the off-line estimation of the model parameters is to evaluate the value of the parameters A, B and Ko that must be stored in the computer for ready access when the process is in operation. The value of these parameters is not constant over the entire operating region of the process. Therefore, the range of operation has to be divided into smaller regions during each of which the parameters have to be estimated separately, and stored in the computer. In each region, for off-line estimation purposes, K can be replaced by Ko in Equation 8, and the model equation becomes,

αi + Aβi + Bγi + Ko = εi (11)

The entire range of operation of the process can be divided into two phases; the warm-up phase, and the run or continuous phase. During the warm-up phase, the maximum possible heater current is applied to the lower zone to bring the temperature of the lower zone from room temperature to the initial operating nutrient-dissolving temperature (approximately 750°F) in the minimum possible time. The heater current applied to the upper zone is controlled during this phase to maintain a steady difference in temperature (approximately 70°F) between the lower and upper zones, the upper zone having the lower temperature. During the run phase, the heater current supplied to the upper and lower zones is controlled to attain the respective temperature profiles shown in FIG. 6.

A typical plot of the temperatures at the upper and the lower zones, against time, during the warm-up phase, is shown in FIG. 10. The current supplied to the lower zone is fixed at its maximum value, for example 4 mA at the input of a current amplifier whose output feeds the heater coils, during the entire warm-up phase. The temperature of the lower zone rises, in a manner which is analogous to the response of a single pole transfer function, to about 400°F. At this point a change of state occurs for the fluid in the autoclave. The material properties of the fluid in this new state are difficult to describe exactly. However, under the extreme temperature and pressure conditions which exist within the autoclave, the well-defined liquid surface of the solvent ceases to exist, and the entire volume within the autoclave becomes filled with a sodium hydroxide -- quartz -- air system. During this transition state, the temperature of the lower zone rises, almost linearly, to about 500°F. After that, the temperature of the lower zone starts rising again, in a manner analogous to the response of a single pole transfer function (of the lower zone) with a different gain and time-constant, until the lower zone attains the initial operating temperature, approximately 750°F. This brings the process to the end of the warm-up phase, at which point the run phase begins. Thus, a natural breakdown of the entire process is as follows: ##SPC4##

The model parameters must be estimated for each of these regions separately. First the lower zone will be discussed and then the upper zone.

Lower Zone

Since the heater current supplied to the lower zone is held constant throughout the entire warm-up phase (i.e., regions, 1, 2 and 3 of FIG. 10), the model described by Equation 11 may be further simplified for the warm-up phase of the lower zone, for reasons described below.

For all practical purposes, γi in Equation 11 remains constant for each value of i since xi remains the same for each value i, and the value of ζi is negligible. Thus, the two terms Bγi and Ko in Equation 11 may be combined into one term, thus reducing the number of parameters to be estimated. If this is not done, there will be no unique solutions for B and Ko. As proof of this, assume the subscript i is dropped from γi in Equation 11, then solve Equation 11 for A, B and Ko to minimize the error ##SPC5##

(where n is the total number of samples). Let A, B and Ko be the least square estimates given by the equation, ##SPC6##

D' = the matrix transpose of D

and n = the number of sample points

Since the second and the third column vectors of the matrix D in Equation 13 have identical directions, matrix [D'D] must be singular. Hence, its inverse, as it appears in Equation 12, does not exist. Since unique solutions do not exist for B and Ko, and since γ is a non-zero constant, then Ko in Equation 11 can be arbitrarily equated to zero. It should be noted that if zero is a feasible value for γ, Ko cannot be assumed to be zero because Bγ + Ko can assume a non-zero value. Thus, Equation 11 becomes,

αi + Aβi + Bγ = εi (14)

The least square estimate of A and B is given by: ##SPC7##

Inserting experimentally derived values of αi, βi and γ in Equation 15, it is found that:

A = 0.0034 for Region 1 of B = 0.0142755 the lower zone (17) and, A = 0.00172 for Region 3 of B = 0.013207 the lower zone (18)

A and B are not evaluated for Region 2 of the lower zone because, as will become evident, the characteristics of the lower zone in Region 2 are not critical for the overall control of the process. Region 2 is assumed to be a continuation of Region 1. The structural errors caused by this assumption were found to be well within acceptable limits.

In Region 4 of the lower zone, γ is no longer constant because the current supplied to the lower zone is controlled so that the lower temperature can track a prescribed profile. Thus, instead of Equation 14, Equation 11 must be used as the model equation. The operating range of the lower temperature in Region 4 is in the neighborhood of 750°F. Accordingly, some test data were collected to study the bottom temperature with an initial value 750°F, while the unamplified current supplied to the lower zone was held fixed at various constant values (for example, 0 ma, 1 ma, 2 ma and 3 ma, respectively). These four current settings can be denoted as uij, where j goes from 1 to 4, and i is the running index of the samples. For any particular j, uij remains constant for every value of i, because the test current signals are step fucntions. Hence, the subscript i can be dropped from uij. The corresponding γj, when inserted in Equation 9, defines the variables Wj where K is replaced by Ko. Equation 11 thus becomes:

αij + Aβij + Wj = εij (19)

where j goes from 1 to 4.

The estimates of A and Wj, j = 1, 2, 3, 4, are obtained so as to minimize E where ##SPC8##

where nj is the number of sample points for the jth current setting.

Let these estimates be denoted by A and Wj, respectively. Differentiating E with respect to A, one obtains ##SPC9##

Differentiating E with respect to Wj, one obtains

Wj = - (Sj + A Rj)/nj (22)

Where, ##SPC10##

Solving Equations 21 and 22 for A, we get: ##SPC11##

Equation 23 is used to evaluate A. The values of Wj, for j = 1, 2, 3 and 4, are obtained by inserting this value of A in Equation 22. FIG. 9 shows the plot of Wj against γj as j goes from 1 to 4. The plot is approximately a straight line. Referring to the definition of Wj, the slope of this straight line is defined as B, i.e., the estimate of B. The ordinate of the point of intersection of the straight line with the Wj axis is defined as Ko, i.e., the estimate of Ko. Based on experimental data, the numerical values of the parameters are given by:

A = 0.00258717 For Region 4 of B = 0.006867175 (24) the lower zone Ko = -0.12539601

Upper Zone

The heater current supplied to the zone is controlled such that, during warm-up the upper temperature is some 70°F below the lower temperature, and that, during the run phase, the upper temperature tracks the prescribed temperature profile. Thus, unlike the current supplied to the lower zone, the current supplied to the upper zone is not held constant in any of the four regions. Hence, instead of Equation 14, Equation 11 must be used as the model for the upper zone, and the parameters have to be estimated in each region.

In Region 1, test data were collected in the same manner as in Region 4 of the lower zone. With room temperature as the initial value of the upper temperature, experiments were run wherein the heater currents, at the input of the power amplifier for the upper zone, were chosen to be step functions of magnitudes 0 ma., 1 ma., 2 ma., and 3 ma., respectively, while the current supplied to the lower zone was held at its usual operating value. In each experiment, the transient value of the upper temperature was recorded. The upper temperature in Region 1 was found to increase linearly with time with a slope which is a linear function of a dummy variable H where H(t) = x2 (t) + x(t). The linear relationship implies that A = 0. Hence, Equation 19 becomes:

αij + Wj = εij

where j goes from 1 to 4 (corresponding to currents 0, 1, 2 and 3 ma.), and Wj = βγj + Ko, where γj has the following new definition:

γj = -Hj Δt, i.e., uj is replaced by Hj.

The least square estimate of Wj is given by Equation 22 with zero substituted for A. The values of B and Ko are read from FIG. 12 in the same fashion as FIG. 10 is read for region 4 of the bottom zone. These values are given by:

A = 0 For region 1 of the B = 0.001382634 (25) upper zone Ko = -0.05991414

In Region 2, when the fluid in the autoclave is undergoing a transition from one state to another, the temperature differential between the upper and lower zones assumes values less than 70°F, see FIG. 10. In other words, the temperature of the upper zone becomes greater than it should be. Nothing much can be done to implement better control of the temperature of the upper zone, because the current supplied to the upper zone in region 2 is invariably zero. It is for this reason that, as previously mentioned, the characteristics of the bottom zone in region 2 are not critical. Inspite of the fact that zero current is supplied to the upper heater, the temperature at the upper zone increases because Ko of Equation 11 changes abruptly as soon as the fluid within the autoclave undergoes a change of state from Region 1 to Region 2. Since the current supplied to the upper zone in Region 2 is always zero, Equation 11 can be simplified as:

αi + Aβi + Ko = εi (26)

The parameter B is irrelevant. The estimates of A and Ko that minimize ##SPC12##

are given by: ##SPC13##

Inserting experimentally obtained values for αi and βi in Equation 27, we obtain

A = 0.006764 For Region 2 of (28) Ko = -0.29436 the upper zone

For Region 3, it was experimentally found that if the values of A and Ko from Equation 28 were inserted in Equation 11, and if B was estimated on a least square basis, the residual error was well within acceptable limits. Hence, the estimate of B is given by: ##SPC14##

where A and Ko are given by Equation 28. Thus, the parameters for the upper zone for Region 3 are given by:

A = 0.006764 For Region 3 of B = 0.02 (30) the upper zone Ko = -0.29436

The parameters for the upper zone for Region 4 are obtained in the same manner as those for Region 1 of the upper zone. The model is described by Equation 11. A typical operating temperature for the upper zone in Region 4 is about 680°F. Accordingly, with 680°F as the initial value of the upper zone the current supplied to the upper zone was held constant at 0 mA, 0.5 mA, and 1 mA, respectively. The transient top temperature data was collected for each case. Equation 23 was then used to evaluate A with the minor alteration that j goes from 1 to 3 instead of going from 1 to 4. FIG. 13 was then used to evaluate B and Ko. The numerical values are given by:

A = 0.00499285 For Region 4 of B = 0.01392232 (31) the upper zone Ko = -0.28624648

The model equations and typical parameter estimations for all regions, for both the upper and lower zones, are set forth in Table B, below. ##SPC15##

It will be appreciated that the numerical values of A, B and Ko, as they appear in Equations 17, 18, 24, 25, 28, 30 and 31, depend on the scale factors of the temperature or current transducers used to measure ui and vi. Thus, there is nothing unique about the particular values stated. In the illustrative embodiment of the invention, the above-mentioned numerical values were obtained using the temperature values as read by a thermocouple where 1°F corresponds to 0.0307252 mv. The heater currents were measured at the input of the power amplifiers whose outputs feed the actual heater coils. The gain of these amplifiers will depend on the number and capacity of the heater elements used for a given autoclave.

Summarizing, the control algorithm can be divided into two phases: (1) the warm-up phase, and (2) the run phase. The objectives to be accomplished during these two phases have already been discussed. During the warm-up phase, the current supplied to the lower zone is held constant at its maximum allowable value, and the current supplied to the upper zone is controlled so as to maintain the upper temperature lower than the lower zone temperature by about 70°F. During the run phase, the heater current supplied to both the top and the bottom zones is controlled so that the temperature profiles shown in FIG. 6 are obtained.

The model representing the on-line operation of the upper zone in each of the three regions of the warm-up phase is given by Equation 8. Depending on the region of operation, the values of A and B are read from Equation 25, 28 or 30, and inserted for A and B in Equation 8. Note that the value of B in Region 2 of the upper zone is immaterial, so long as it remains positive. Hence, the value of B from Equation 25 is allowed to remain in the memory region of the computer for Region 2 also. The initial value of K for Equation 8 in each region of operation, is assumed to be Ko of the corresponding region as given by Equation 25, 28 or 30. As the process operation continues, the value of K is updated by equating it to K of Equation 10. The observation time for this updating is chosen to be thirty minutes. Since the interval between two consecutive samples is three minutes, the value of n in Equation 10 is 10. A direct comparison of coefficients of Equations 7 and 8 reveals that:

αi = vi+k+1 - vi+k ; βi = 1/2[ vi+k+1 + vi+k ]Δ t

γi = -ui Δt ; K = -CΔt (32)

The index k is equal to 5 for the upper zone, since the transport delay in the top zone is 15 minutes. Thus, ui for the upper zone (i.e., the required value of the current supplied to the upper zone at the present instant) can be obtained by solving Equation 7, provided proper predictions are made for yi+k+1 and yi+k of the upper zone. To obtain these predicted values, the lower zone temperature was predicted using the model for the lower zone. Once yi+k+1 and yi+k for the lower zone were found, the desired temperature difference between the lower and the upper zones can be subtracted from the predicted temperatures of the lower zone to obtain yi+k +1 and yi+k for the upper zone. However, this control scheme was found to be not entirely satisfactory due to the fact that the value of ui became highly sensitive to the nondeterministic terms, that is to ζi, as defined in Equation 6 and to εi. For a small positive value of (αi - εi), in the order of 1°F, the required ui was found to be higher than the maximum allowable value. Similarly, for a small negative value of αi, ui must be set equal to zero. Thus, the control algorithm resembles a "bang-bang" situation. As previously mentioned, if a rectangular current wave, of width T, is supplied to the heating coils, the resulting temperature function will be spread out and have a width much greater than T. This phenomenon, coupled with the sensitivity of ui with respect to the nondeterministic terms, gave rise to oscillations in the process operation.

One pragmatic way to solve this problem is to decrease the sensitivity of ui with respect to the nondeterministic terms. This was done by predicting the temperatures of the lower zone by: yi+k+1, yi+k+2, . . . , yi+k+r where r is a positive integer whose value must be empirically determined. A typical value of r is 15 . From these lower zone temperature values, the corresponding values for the desired upper zone temperature are determined. Let these upper zone temperature values be denoted by di+k+1, di+k+2, . . . , di+k+r. Now, ui for the upper zone is determined so as to attain di+k+1, di+k+2, . . . , di+k+r in a minimum mean square sense. The mathematical formulae for implementing this scheme are developed below.

The on-line model for the bottom zone is given by Equation 8. Combining Equations 6, 7 and 8, we obtain:

yi+1 = [ 2/2 + AΔt] [ yi { 1 - AΔt/2}- Bγi-k - K] + . . .

a stochastic term (33)

where the subscript i corresponds to the present instant of time. Inserting the estimated values of the parameters in Equation 33, and neglecting the stochastic term which constitutes the error term, we obtain:

yi+1 = [ 2/2 + AΔt] [yi {1 - AΔt/2} - Bγi-k - K] (34)

where yi+1 is the predicted value of yi+1 for the lower zone, A and B are read from Equation 17 or 18, depending on the region of operation, and K is evaluated on an on-line basis using Equation 10 . The initial value of K (i.e., Ko) is equated to zero as previously discussed. The index k corresponds to the transport delay in the lower zone. Equation 34 is used in a recursive fashion to evaluate the required predicted values of the lower zone temperature. To start the recursion, yi is equated to vi of he lower zone (i.e., the measured value of the lower zone temperature at the present instant of time). The desired upper zone temperatures di+k+1, di+k+2, . . . , di+k+r are evaluated directly from the predicted bottom temperatures yi+k+1, yi+k+2, . . . , yi+k+r.

The on-line model for the upper zone, as mentioned earlier, is given by Equation 8 which, combined with Equations 6 and 7, results in Equation 34, where yi+1 is the predicted value of the upper zone temperature and A, B, and K correspond to the upper zone.

Define three variables, μ, ν, ω as:

μ = [2/2 + AΔt] (1 - AΔt/2) (35) ν = 2/2 + AΔt (36)

ωi-k = - Bγi-k - K (37)

thus, Equation 34 for the top zone can be rewritten as:

yi+1 = μyi = νωi-k (38)

Assume that ωi = ωi+1 = . . . ωi+r-1, where r is the previously defined positive integer. This assumption is a reasonable approximation, because the value of ω has not been, in practice, found to vary widely over r sampling intervals. The problem is to find that value of ωi which minimizes: ##SPC16##

where yi+j is the predicted value of the upper zone temperature. Differentiating J with respect to ωi and equating to zero, we obtain ##SPC17##

where j and m are running indexes for summation, yi+k is the predicted upper zone temperature, and the superscripts indicate the power (in the sense that 2 would imply a square). Equation 38 is used to evaluate yi+k starting the recursion with yi = vi. From Equation 37 it follows that the required value of the upper current at the instant i is given by:

ui = (ωi + K)/(BΔt) (41)

where ωi is given by Equation 40.

During the run phase, both the currents supplied to the upper and lower zones are independently controlled so that the temperatures at the upper and at the lower zone may respectively follow their prescribed profiles. Thus, Equation 40 is used for the upper as well as for the lower zone, where the term di+k+j denotes the sample value of the corresponding prescribed profile. μ and ν are evaluated separately for the lower and for the upper zone using Equations 35 and 36, where the value of A is read from Equation 24 for the lower zone and from Equation 31 for the upper zone. Equation 38 is used to evaluate yi+k where the values of μ, ν and k for the upper zone are used to calculate the upper zone current, and their values for the lower zone are used to calculate the lower zone current. Knowing the values of ωi for the upper and the lower zones, Equation 41 is used to determine the required upper and lower zone currents. The values of B and K in Equation 41 are used accordingly, i.e., Equation 24 is used for the lower zone, and Equation 31 is used for the upper zone. A flow diagram of the entire control algorithm is shown in FIG. 16.

Summary of the Algorithm

First the algorithm performs a series of tests to determine whether the crystal-growing process is operating in one of the three regions of the warm-up phase, or in the run phase, and based upon that determination, selects the appropriate values of the estimated off-line parameters, recalculating the parameter K, as need be, as time proceeds. The algorithm also directs that continuous samples of the autoclave temperatures and of the heater currents supplied thereto be taken and, in addition, based upon the results of these samples, continuously predicts the desired temperatures for the upper and lower zones of the autoclave, using either the fixed differential DELT, or the difference mandated by the temperature profile stored in the memory region of the computer. Finally, with appropriate consideration for the thermal lag in the system, the algorithm computes the necessary currents which must be supplied to the heater elements to produce these desired temperatures, at some interval in the future.

Detailed Description of the Algorithm

As shown in FIG. 16, the first step in the logical flow chart is to define the parameters used in the flow chart. The abbreviation INT represents the time interval between two successive samples of temperature and current, for both the upper and lower zones of the autoclave. Typically, this interval will be three minutes.

With regard to the parameter R, it will be recalled that there is a certain transport delay associated with the mathematical model for the autoclave. This is a measure of the thermal inertia of the autoclave, which has substantial physical mass. As discussed previously, it is necessary to anticipate the reaction of the autoclave to changes in heater current. For example, a step change in current will not effect a change in autoclave temperature until some finite interval of time, typically twelve minutes has elapsed. This means, therefore, that the algorithm must predict temperatures for quite some time in the future so that the desired currents which will establish those temperatures can be computed.

Unfortunately, during some portions of the warm-up or run phases, the desired temperature approaches very close to the predicted temperature and a direct application of the mathematical model would suggest extremely high heater currents which, for many reasons, are not practical. Thus, the algorithm is structured to assume that, because the autoclave temperature profile passes through the desired temperature profile, it is possible to calculate those autoclave heater currents which will cause the actual autoclave temperature to approach, in a least square sense, the desired temperature profile. Thus, the parameter R represents the number of temperature points on the profile, which must be considered, in the future, to calculate a present current value.

The abbreviation DELT represents the temperature differential between the upper and lower autoclave zones, during the three regions of the warm-up phase.

With regard to the abbreviation NOBS, it will be recalled that the parameters A and B in the mathematical model for the autoclave are assumed to be constant, for both the upper and lower zones of the autoclave, during each of the three regions of the warm-up phase and the run phase, but that the parameter K must be continually updated. It will also be recalled that the current supplied to the heating coils, and the temperature of the autoclave, are continuously sampled during the growth cycle, for example, every three minutes. As explained in the preceding mathematical analysis, the parameter K must be updated approximately every half hour. Thus, approximately ten discrete samples of temperature and current are required before the algorithm has enough data to accurately update the parameter K. The factor NOBS, then, simply represents the number of samples of temperature and current which must be accumulated by the computer before K may be safely updated.

As previously discussed, the crystal-growing process is divided into a warm-up phase and a run phase. The warm-up phase is, in turn, subdivided into three regions, defined essentially by the temperature attained by the lower autoclave zone. For example, as shown in FIG. 10, the first region extends from room temperature to approximately 400°F, the second region extends from approximately 400° to approximately 500°, the third region extends from approximately 500° to approximately 750°, and the run phase extends from approximately 750° to the temperature required by the desired temperature profile. The parameters S1, S2, and S3, therefore, are included in the algorithm as logical switches which permit the algorithm to determine whether it is operating in any one of the three regions of the warm-up phase or in the run phase. Thus, typically S1 will have a value of 400°, S2 will have a value of 500°, and S3 will have a value of 750°. It should be born in mind, of course, that the actual temperature value at which these switches are set will vary according to the type of crystal being grown and the physical arrangement within the autoclave.

NT and NB represent logical flags and are representative of the transport delay in the upper and lower autoclave zones. For example, if the interval INT between successive samples is three minutes, NT = 5 means that the transport delay associated with the upper autoclave zone must be 15 minutes, and if NT = 4, the transport delay associated with the lower zone must be 12 minutes.

T Data (l), and B Data (l) represent the desired temperature profile for the upper and lower autoclave zones, respectively, over the entire 25-day run phase, the data being defined for every INT, i.e., for each three-minute interval of the entire run phase.

As shown in FIG. 16, the next step in the logical flow chart is to enter into the computer the particular formula which is to be used to predict the temperature of the autoclave, as a function of the current supplied to the heating elements associated therewith. This formula is set forth in FIG. 11. In like fashion, the formulae which are used to compute the heater currents, as a function of the desired temperature and the appropriate parameters A, B, and K for the phase and autoclave zone under consideration, are set forth in FIG. 12. Further, the formulae used to update the parameter K are set forth in FIG. 13, and are also entered into the computer at this time.

Next, the actual values for the estimated off-line parameters A, B, and K, for the upper and lower autoclave zones, for all three regions of the warm-up phase and the run phase, are now entered into the computer. The actual values for the illustrative quartz-growing process disclosed, are set forth in FIG. 14.

It will be appreciated that the actual values of INT, R, DELT, NOBS, S1, S2, S3, NT, NB, T Data (l), B Data (l) and the values of A, B, and K are functions of the type of crystal being grown and the physical characteristics of the particular autoclave used.

The next step in the logical flow chart is to initialize the values of the parameters M1, M2, and M3 at -1 at l = 1. The parameters M1, M2, and M3 are flags which are included as tests so that the algorithm, as it proceeds, can determine whether or not it is operating in the first, second or third region of the warm-up phase, or in the run phase. For example, the algorithm is so structured that when M1 changes from its initialized value of -1 to 0, the algorithm knows that the process has just crossed the boundary between region 1 of the warm-up phase and region 2 of the warm-up phase. Similarly, when M2 changes from -1 to 0 the algorithm knows that the process has just crossed the boundary between region 2 and region 3 of the warm-up phase, and when M3 changes from -1 to 0, the algorithm knows that the process has just crossed the boundary between region 3 of the warm-up phase and the run phase. The parameter l is, of course, the running index of the given, desired temperature profile data, as defined for T Data () and B Data (l), above.

The next step in the flow chart is to define a continuously shifting time index, i, which is shorter in duration than l. This is done, as a matter of expediency, to economize the use of the memory section of the computer. The reason that a continuously shifting time index, shorter than the cumulative time index l, is required is that only approximately 30 sampled values of current and temperature need be stored in the computer, at any given time, as the parameters of the process are continually being updated. Thus, storage of more than 30 samples of temperature and current is redundant, thus, as the shifting time index i reaches the value 30, the earliest stored temperature and current samples are erased. Clearly then, the first step in the algorithm is to initialize i at 1.

The factors Tx and Bx are the two principle variables in the process, and represent the currents which are supplied to the upper and lower zones of the autoclave, respectively. Since, at this stage of the algorithm, the growth process has just begun, Txx(i) is, in fact, Tx(l), and Bx(i) is Bx(l), and both are set equal to an initial value of zero.

Now that these housekeeping chores have been attended to, the algorithm begins to control the actual crystal-growing process. The first step in the algorithm is to read the first value of Bv(i), that is, Bv(l), where Bv represents the actual temperature of the lower zone of the autoclave. After this temperature has been read, the computer performs a logical calculation asking the question, "Is Bv(l) the actual temperature of the lower zone, less than the value of the switch S1 ?" which, as we have previously seen, typically has a value of 400°F. If the answer to this question is "yes," this means that the algorithm is now operating in the first region of the warm-up phase.

Accordingly, the next step in the algorithm is to set a flag, called NFLAG equal to 1. As will be seen later, NFLAG is used to select the appropriate values of A, B and Ko, the estimated off-line parameters for the various regions and phases of the process, for use in the temperature predicting and current calculating formulae. The parameter NFLAG differs from the parameters M1, M2, M3 in that these latter parameters indicate when a transition between regions or zones has occurred whereas NFLAG is merely an indicator telling which region or zone the process is operating.

Thus, if the answer to the preceding question is "no," the algorithm then asks another question, i.e., "Is the actual value of the temperature of the lower zone less than S2 ?" where S2 is the second switch, which typically has a value of 500°F. If the answer to this second question is "yes," this means that the algorithm is operating in the second region of the warm-up phase. Accordingly, NFLAG is set equal to 2. If the answer to the second question is "no," the algorithm then asks a third question, namely, "Is the actual value of the temperature of the lower zone less than the switch S3 ?" where S3, the third switch typically has a value of 750°F. If the answer to this third question is "yes," then the algorithm is operating in the third and final region of the warm-up phase. Accordingly, the algorithm sets NFLAG equal to 3. On the other hand, if the answer to the third question is "no," then obviously the algorithm is no longer operating in any of the three regions of the warm-up phase, but must now be operating in the run phase. Accordingly, NFLAG is set equal to 4 so that the appropriate values of estimated off-line parameters will be employed whenever any of the formulae set forth in FIGS. 11 through 13 are employed.

Since it has now been determined in which region or phase the algorithm is operating, the next step in the flow chart is to read in the measured values of the temperature of the upper autoclave zone and the lower autoclave zone, Tv(i) and Bv(i), respectively. It will be noted that the algorithm is re-entrant at this point, as designated by the circle letter A. This will be referred to later.

The next step is to ask the question, "Is NFLAG = 4?" which condition, if true, would indicate that the algorithm is operating in the run phase. If the answer to this question is "no," as it will be during any of the three regions of the warm-up phase, the algorithm then directs that the heater current for the lower zone of the autoclave, Bx(i), be set to its maximum safe value.

It will be recalled that during the warm-up phase, the heater current supplied to the lower zone of the autoclave is set at its maximum value, to minimize warm-up time, and only the upper zone current is controlled, to maintain a constant temperature differential between the upper and lower zones of the autoclave.

After this step has been taken, or in the alternative, if the question previously asked indicates that NFLAG is indeed 4, the next step in the algorithm is to read the measured values of heater current, Tu(i) and Bu(i), for the upper and lower zones of the autoclave, respectively. The formulae used to predict temperature, compute the heater current and update the value of the parameter K are generalized for use with both the upper and lower zones of the autoclave. This permits the equations to be written as subroutines which may be stored in the computer memory and called at will. Thus, the next step in the algorithm is to set the parameter k (FIGS. 11, 12, and 13) equal to NB, where NB, it will be recalled, is a number indicative of the number of intervals INT corresponding to the transport delay of the lower zone of the autoclave.

It will also be recalled that the temperature predicting formula set forth in FIG. 11 predicts the temperature of the autoclave one interval INT in advance. It is, of course, necessary to keep a record of how often the temperature prediction formula has been called, since it is only necessary to predict NB temperature points in advance for the lower zone, and NT temperature points in advance for the upper zone. Thus, j represents a running index of the iteration of the temperature predicting formula, and as soon as j = k, the iteration is terminated. The parameter k will, of course, equal NB or NT, depending on the current autoclave zone of interest. Here, since the algorithm is concerned with the lower zone of the autoclave, the algorithm sets k equal to NB and j is initialized at 0.

The parameter y(i) also appears in the generalized temperature predicting formula and the generalized current calculating formula. Thus, the algorithm sets y(i) equal to Bv(i), the observed temperature of the lower autoclave zone.

The next step in the flow chart is to enter the estimated off-line parameters for the lower zone of the autoclave which correspond to the current value at which NFLAG is set. These parameters are, of course, used in both the temperature predicting formula and the current calculating formulae. Once these parameters have been called from the storage area of the computer into the central processing unit, the algorithm asks the question, "Is i = NOBS?" where NOBS represents the number of samples between two consecutive updatings of the process parameter K. If the answer to this question is "no," as will be the case initially, the next step dictated by the flow chart is to compute the temperature of the lower zone of the autoclave using the temperature predicting formula of FIG. 11, with the appropriate values for the parameters A, B and K and with u set equal to Bu, that is to say, with the measured value of the current equal to the measured value of the current in the heaters associated with the lower zone of the autoclave.

If the answer to the previous question, "Is i = NOBS?" is "yes," this means that it is now appropriate to update the value of the parameter K. Accordingly, the flow chart directs that K be updated using the formula set forth in FIG. 13 with v equal to the measured value of the lower zone temperature, and u equal to the measured value of the lower zone heater current. Then, using this updated K, the algorithm proceeds to compute the temperature of the lower zone using the temperature predicting formula set forth in FIG. 11, as previously discussed.

The next step in the algorithm is to compute the desired upper zone temperature profile, d(i + j + 1), which is equal to the predicted lower zone temperature profile y(i + j + 1) minus DELT, the desired temperature differential between the upper zone and the lower zone, typically 70°F.

As shown in FIG. 18, the next step in the algorithm is to ask the question, "Is j<k - 1?" where j is the running index of iteration of the temperature predicting formula, and k is representative of the transport delay corresponding to the lower zone of the autoclave, which was priorly set equal to NB. If the answer to the preceding question is "yes," the present value of j is increased by a factor of one, and the flow chart becomes re-entrant, as shown in FIG. 17. The algorithm, therefore, continues to predict the temperature of the lower zone until the condition j = k occurs, at which time the answer to the question, "Is j<k - 1?" is "no," causing the algorithm to pose the question, "Is NFLAG = 4?" The purpose in asking this latter question is to determine whether the algorithm is operating in the warm-up phase or in the run phase. If the result of this question indicates that the algorithm is now operating in the run phase, it is no longer proper to predict the upper zone temperature by subtracting the fixed quantity DELT from the lower zone temperature. Rather, it is necessary to refer to the desired temperature profiles, as stored in the computer memory, as will be explained later.

Assume now that the answer to the question, "Is NFLAG = 4?" is "no." This means that the algorithm must still be operating in the warm-up phase; accordingly, the algorithm then asks the question, "Is y(i + j + 1)<S1 ?" that is to say, is the predicted temperature y of the lower zone less than the value of the first switch S1. In other words, is the algorithm still operating in the first region of the warm-up phase? If the answer to this question is "yes," and the algorithm is still operating in the first region, the flow chart then poses the question, "Is M1. M2 = 0?"

It will be recalled, that initially M1, M2 and M3 were set equal to -1, and that M1 serves as a test to see if the boundary between the first and second regions of the warm-up phase has been crossed. The reason that the algorithm provides an additional test to determine which region of the warm-up phase the process is operating in is that the continuously shifting time index i, must be reset immediately that the parameters of the process cross the border between one warm-up region and another. This is so because i is used to update the parameter K every NOBS interval, and if i were not reset each time a border between one region and another were crossed, the algorithm would incorrectly attempt to update the value of K at the end of the next ith samples.

Thus, the expression, "Is M1.M2 = 0?" is merely an expedient way to test to see if either M1 or M2 is 0, which would be the case if a border between warm-up regions had been crossed. If the answer to this test is "yes," i is initialized again at 1. If the answer is "no," i is not re-initialized, but retains its present value and k, which, as we have previously seen, is one of the parameters in the generalized temperature predicting formula and the generalized current computing formula, is set equal to NT, the transport delay corresponding to the upper zone of the autoclave.

The algorithm then asks the question, "Is j<k + R - 1?" It will be recalled, that j is the running index of reiteration of the temperature predicting formula. This current test, therefore, determines whether or not the temperature has been predicted, k + R - 1 times. If the answer to this question is "yes," this means that the algorithm has not yet computed enough predicted temperature points for the lower autoclave zone, and the algorithm becomes re-entrant at the instruction "Set j = j + 1." This instruction, in turn, re-enters the algorithm at the top of the temperature predicting loop, i.e., at the instruction, "Call estimated off-line parameters corresponding to NFLAG."

Thus, the algorithm continues to predict the lower zone temperatures until the condition occurs where j becomes equal to k + R, which means, in effect, that the algorithm now has sufficient predicted temperatures for the lower zone, and that the algorithm may now proceed to compute the temperatures for the upper zone, by subtracting the constant DELT if the algorithm is operating in the warm-up phase, or by reference to the stored upper and lower zone temperature profile if in the run phase. In any event, the algorithm next directs that j be set to 0, and that y(i), the predicted temperature be set equal to the measured temperature for the upper zone, Tv(i).

Let us return now to the question that the algorithm posed earlier, namely, "Is y(i + j + 1)<S1 ?" If the answer to this question is "no," this can only mean that the algorithm is no longer operating in the first region of the warm-up phase. Accordingly, the algorithm directs that M1 be increased in value by 1, and that NFLAG be set equal to 2. The algorithm then asks the further question, "Is y(i + j + 1)<S2 ?", that is, is the predicted temperature for the lower zone, y(i + j + 1) less than the switch S2, which delimits the second region of the warm-up phase. If the answer to this further question is "yes," then the algorithm knows that it is operating in the second region of the warm-up phase and the algorithm re-enters at the previously discussed instruction which serves as a test to determine whether or not i should be reset to 1; namely, at the question, "Is M1. M2 = 0?"

On the other hand, if the answer to the question, "Is y(i + j + 1)<S2 ?" is "no," this must mean that the algorithm is no longer operating in the second region of the warm-up phase, but has shifted to the third region of the warm-up phase. Accordingly, M2 is increased in value by 1, and the value of NFLAG set equal to 3. The algorithm then asks the question, "Is y(i + j + 1)<S3 ?" If the answer to this latest question is "yes," the algorithm must again be operating in the third region of the warm-up phase and the algorithm once again re-enters at the instruction which serves as a test to determine whether or not i should be reinitialized, i.e., "Is M1. M2 = 0?"

On the other hand, if the answer to the question, "Is y(i + j + 1)<S3 ?" is "no," this must mean that the algorithm is no longer operating in the third and final region of the warm-up phase. Accordingly, the algorithm branches and sets M3 equal to M3 + 1, as shown in FIG. 21. The question is now asked, "Is M3 = 0?" This again is a test to determine whether or not the boundary between the third region of the warm-up phase and the beginning of the run phase has been crossed. If the answer to that question is "no," and the algorithm is, therefore, not operating at the boundary between the third region of the warm-up phase and the run phase, j is set equal to 1. The algorithm then proceeds to compute the desired temperature for the lower zone D(i + k + j) by equating this quantity to B Data (l + j - 1), i.e., the corresponding value from the desired temperature profile priorly stored in the memory region of the computer. On the other hand, if the answer to the question, "Is M3 = 0?" is "yes," this indicates that the algorithm has just crossed the boundary between region 3 of the warm-up phase and the run phase. Accordingly, the continuously shifting index i must be reset to 1 and NFLAG set equal to 4. Further, the estimated off-line parameters corresponding to the current value of NFLAG are entered from the memory area of the computer into the central processing unit.

The algorithm then re-enters at the instruction "Set j = 1" and continues to compute the desired temperatures for the lower zone of the autoclave in the run phase. Next, the algorithm asks the question, "Is j<R?", that is, is the running index of iteration of the temperature predicting formula i less than R, where R represents the total number of samples which must be taken in the future to fit, in a least square sense, the desired temperature profile. If the answer to this question is "yes," i.e., if j is less than R, the value of j is increased by one and another desired temperature point is computed. This process is re-iterated until j becomes equal to R, at which point the algorithm poses yet another question, "Is i = NOBS?", i.e., is i equal to the total number of desired observation points between two successive updatings of the parameter K. If the answer to this question is "no," the algorithm then proceeds to compute the heater currents for the lower zone of the autoclave which correspond to the desired temperatures previously calculated. After these heater currents are computed, signals are transmitted to the magnetic amplifiers associated with the heaters of the lower zone, as discussed with reference to the description of the overall apparatus.

If the answer to the question, "Is i = NOBS?" is "yes," this means that the required number of desired temperatures have been calculated, and it is now necessary to update the value of the parameter K, which is, of course, done by means of the formula set forth in FIG. 13, with v set equal to Bv, the lower zone temperature, and u equal to Bu, the lower zone current. The algorithm then re-enters at the current computing instruction, previously discussed, and continues by setting k = NT, the transport delay associated with the upper zone, and by re-initializing the value of j, the running index of reiteration for the temperature predicting formula, at 1. Then, the algorithm begins to compute the desired temperatures for the upper zone of the autoclave according to the formula d(i + k + j) = T Data (l + j + NT - NB - 1), where this latter term represents the desired temperature points for the upper zone of the autoclave as set forth in the temperature profile for the upper zone, priorly stored in the memory area of the computer. The factor, l, is, of course, the running index for the entire run phase; the quantity j is the running index of reiteration for the temperature predicting formula; NT and NB represent the transport delay of the upper and lower zones of the autoclave, respectively; and the expression, (l + j + NT - NB - 1) merely establishes a matching relationship between the indices of the upper zone and lower zone sample values which are, of course, offset from one another by the differing indices and delay times.

Next, the algorithm again asks the question, "Is j<R?" If the answer to this question is "yes," j is increased by 1 and the algorithm re-enters at the instruction at which the desired upper zone temperatures are computed. This process is reiterated until j becomes equal to R, at which point a sufficient number of upper zone temperatures have been calculated, so the algorithm proceeds by increasing the value of l, the running index for the entire run phase, by 1, and re-entering the loop at point B of the algorithm as shown in FIG. 18. Thus, when j becomes equal to R, the algorithm becomes re-entrant at the instruction "Set j = 0, y(i) = Tv(i)."

Next, as shown in FIG. 19, the estimated off-line parameters corresponding to NFLAG + 4, i.e., the values of A, B, and K for the upper zone of the autoclave for the warm-up and run phase, as the case may be, are entered from the memory region of the computer into the central processing unit.

The algorithm then asks the question, "Is NFLAG = 1?" If the answer to this question is "yes," i.e., if the algorithm is operating in the first region of the warm-up phase, a fictitious variable, Tu(i) is defined equal to Tu2 (i) + Tu(i), and this fictitious variable is employed as if it were the actual measured current supplied to the upper zone of the autoclave. As explained previously, this fiction is resorted to only in the first region of the warm-up phase, and for all other regions and phases, the actual value of Tu(i) is employed. On the other hand, if the answer to the question, "Is NFLAG = 1?" is "no," the algorithm directs no change in the value assigned to Tu(i) and then asks the question, "Is i = NOBS?" If the answer to this question is "yes," it is now necessary to update the value of the parameter K using the formula set forth in FIG. 13, with v equal to Tv and u equal to Tu. If, however, the answer to this question is "no," then K is not updated, at this time, and the algorithm proceeds to compute the desired temperatures for the upper zone of the autoclave, using the temperature predicting formula in FIG. 11 with u equal to Tu and the appropriate estimated off-line parameters A, B, and K.

Next, the algorithm asks the question, "Is j<k - 1?" that is, is j, the running index of iteration for the temperature predicting formula, equal to k? It will be recalled, that k was priorly set equal to NT, the transport delay of the upper zone. If this test indicates that j is smaller than k, then the value of j is increased by one, and the algorithm re-enters at the instruction which computes the desired temperatures for the upper zone of the autoclave. This procedure is reiterated, as often as need be, until j becomes equal to k, at which point the algorithm enters the formulae used to compute the heater currents (FIG. 12) from the memory area of the computer and proceeds to compute the appropriate values for these heater currents. Recalling that it was previously found necessary to use a fictitious value for the measured heater current for the upper zone during region 1 of the warm-up phase only, the algorithm asks the question, "Is NFLAG = 1?" If the answer to this question is "yes," and the algorithm is operating in region 1 of the warm-up phase, the calculated value for the heater current, x(i) will be much higher than it should be and must, therefore, be converted back to it original form by solving the inverse equation to the one previously used to define the fictitious value of u(i).

After this has been done, the algorithm directs that the appropriate signals be sent to the D to A converter and to the magnetic amplifier associated with the upper zone of the autoclave. The same is true if the answer to the question, "Is NFLAG = 1?" is "no," indicating that the algorithm is not operating in region 1 of the warm-up phase, in which case it is not necessary to use a fictitious value for the sampled upper zone heater current, u(i) and the computed heater current x(i) needs no conversion.

Finally, the algorithm asks the question, "Is l>lmax ?" i.e., is the cumulative index for the run phase greater than the maximum index value for the entire crystal-growing process? If the answer is "yes," then the process may now be terminated, as obviously the crystals are fully grown. If the answer is "no," however, the algorithm continues by again asking the question, "Is i < NOBS?" If the answer to this question is "yes," this indicates that the continuously shifting time index, i, which is shorter than l, must be increased by a value of 1, and the algorithm becomes re-entrant at point A in FIG. 17, that is, at the instruction to read the measured temperature values Bv(i) and Tv(i) for the lower and upper zones of the autoclave, respectively. At this point the algorithm recycles, predicting temperatures and calculating currents, as previously described, until the time that, with l still smaller than lmax, the question, "Is i<NOBS?" is answered "no." When this condition occurs the algorithm goes through a procedure of establishing "temporary" values for the current and temperature samples for the upper and lower zones of the autoclave. Basically, this procedure merely alters the value of the index associated with each particular recorded sample to account for the fact that i is a continuously shifting index. For example, the measured lower zone current Bu(30) is changed to Bu(29) so that the next incoming sample of heater current may be assigned the index Bu(30). Next, as part of the same procedure, the algorithm changes the index i so that it becomes equal to i - 1, and asks the question, "Is i = 1?" If the answer to this question is "no," the algorithm again re-enters the temporary index shifting mode and reiterates around this loop until i becomes 1, at which point the "temporary" values for the current and temperature samples, for both the upper and lower zones, become the "permanent" values.

The algorithm then asks the question, "Is i = NOBS - 1?" If the answer to this question is "no," then i is set equal to i + 1, and the procedure whereby the "temporary" current and temperature samples become the "permanent" values is reiterated until i = NOBS - 1, at which point the new "permanent" values are entered in the memory area of the computer, and the program re-enters at the instruction where i is set equal to i + 1, the algorithm proceeding to point A, as previously discussed.

For the sake of clarity, various subroutines which print out error conditions, warning signals, and the like for the operator, or which shut down the system to avoid catastrophic failure of the autoclave, if temperatures or pressures become excessive, have not been shown, as these are well known in the art.

The preferred embodiment of the invention, discussed above, proceeded on the assumption that it was desired to produce a crystal having a uniform Q throughout all crystallographic regions of the crystal. To that end, FIG. 6 shows a temperature profile which is calculated to yield a constant flow of solvated quartz over the seed crystal and hence, to yield a uniform growth rate in the crystal. One skilled in the art, however, will appreciate that once the growth process for a particular type of crystal has been characterized, it is a relatively easy matter to substitute some other desired temperature profile to grow a crystal having a non-uniform Q throughout all crystallographic regions of the crystal. One could, if one wished, deliberately set out to manufacture a crystal having a Q which steadily decreases towards the extremities of the crystal or one in which the Q is initially low, then rises to a maximum and finally falls off again to the initial low value. Such crystals could be made by merely modifying the overall shape of the temperature profiles shown in FIG. 10. The algorithm shown in FIGS. 16 to 23 would require essentially no modification to handle such crystals.

FIG. 15 depicts an illustrative arrangement of an apparatus which may be used, according to this invention, to grow synthetic quartz crystals of substantially uniform Q. It must again be emphasized, however, that the invention is not limited to the hydrothermal growth of synthetic quartz, but may advantageously be used to grow many other types of crystalline materials. In that event, it may prove necessary to make modifications to the apparatus shown, but these modifications would be well within the competence of one skilled in the art.

As shown in FIG. 15, heater bands 112 and 113, which are associated with the upper and lower zones of autoclave 101, are connected to the output of a first and a second magnetic amplifier 131 and 132, respectively. Magnetic amplifiers 131 and 132 are connected, in turn, via leads 133 and 134, respectively, to a power supply 136, which supplies the extremely high current needed to raise the temperature of autoclave 101 to well over 600°F. Magnetic amplifiers 131 and 132 are, of course, required because the digital control apparatus, to be described below, is not ordinarily capable of directly handling the extremely high power levels required to drive heating elements 112 and 113.

The input of magnetic amplifier 131, which controls the heating current which is supplied to the upper zone of the autoclave, is connected, via a lead 137, to a first digital-to-analog converter 138, thence, via a lead 139, to a process control output circuit 141, which is associated with a digital computer 142. In a similar manner, the input of magnetic amplifier 132 is connected, via a lead 143, to a second digital-to-analog converter 144, thence, via a lead 146, to output circuit 141. As will be apparent from later discussions, low level digital signals, representative of the heater currents that it is desired to feed to heater elements 112 and 113, are converted by digital-to-analog converters 138 and 144, into analog form and are amplified by magnetic amplifiers 131 and 132, respectively, to the levels necessary to drive the heater elements.

Digital computer 142, typically comprises a central processing unit 151, which is associated with, and coupled to, a memory storage device 152, a control circuit 153, a clock circuit 154, an input device 156, an output device 157, a process control input circuit 158, and as previously discussed, a process control output circuit 141.

Digital computer 142 may be any of several well known commercially available computers, and thus, per se, forms no part of this invention. See, for example, the digital computer disclosed in U. S. Pat. No. 3,400,371 and the book entitled, "Computers: Their Operation and Application," by Berkeley and Wainwright, Rheinhold Publishing Company, New York, New York.

Pressure transducer 114, as previously discussed, is connected through seal 111 into the body of autoclave 101 and is connected, via a lead 161, to the input of a first analog-to-digital converter 162, thence, via a lead 163, to process control input circuit 158. Thus, as will be discussed below, an analog signal indicative of the pressure within the autoclave is converted into digital form for use within computer 132. In a similar manner, temperature transducers 116 and 117, which may typically comprise a thermocouple and associated reference junction are connected via leads 164, 166 to the inputs of a second and a third analog-to-digital converter 167, 168, thence, via leads 169 and 171, respectively, to process control input circuit 158. Thus, analog signal indicative of the temperatures of the upper and lower zones of the autoclave are converted into digital form for use in the computer.

Memory storage device 152 may comprise, for example, an array of magnetic cords, a magnetic disc or drum or a magnetic tape device or any combination thereof. Input device 156 may comprise, for example, a Hollerith card reading device, a punched paper tape reader or the keyboard of a teletypewriter. Output device 157 may comprise, for example, a Hollerith card punch, a paper tape punch or the printer portion of a teletypewriter. A clock circuit is connected to central processing unit 151 via control circuit 153 and supplies the necessary timing pulses to computer 142 to insure synchronous operation of the various components of computer 142 and, in addition, maintains a running record of the crystal-growing process.

In operation a computer program comprising essentially the logic flow chart shown in FIG. 16 is entered into the memory storage device 152 by means of input device 156. In addition, the estimated off-line parameters for the upper and lower zones of the autoclave for the three warm-up phases and the run phase are also entered into the memory end of the computer. This information may be in any suitable form, for example, in machine language or in one of the higher level languages, such as FORTRAN or COBOL, depending upon the type of computer employed.

The central processing unit 151 operating under control of the program stored in memory device 152 performs the necessary calculations corresponding to the particular phase of the process, and forwards the resulting digital signals via lines 139 and 146 to digital-to-analog converters 138 and 144. These converters convert the digital signals to analog form and forwards them to magnetic amplifiers 131 and 132, respectively. Greatly amplified versions of the analog signals and then applied to the heating elements 112 and 113, respectively, to control the heating of the upper and lower zones of the autoclave. The temperature in the upper and lower zones of the autoclave and the pressure within the autoclave are continuously monitored by transducers 114, 116 and 117 and the analog outputs from these transducers are converted by analog-to-digital converters 162, 167 and 168 and applied in digital form to process control input circuit 158. As dictated by the stored program in computer 142, these parameters are sampled from time to time and are compared with the desired temperature and pressure parameters as stored in storage device 152, computer 142 making the necessary computations to alter the heater currents supplied to the autoclave if the temperatures depart from their desired values. Clock circuit 154 determines the interval between successive samples of temperature and pressure and also maintains a continuous running record of the crystal-growing process so that the appropriate estimated off-line parameters may be called from the memory storage device 152 into the central processing unit 151 to perform the appropriate computations. After clock circuit 154 indicates that the run phase has been completed, computer 142 will disconnect all power to the heaters and signal the operator via output device 157 teat the growth cycle has been completed.

Advantageously, the program which is inserted into memory storage device 152 may include appropriate error signals and to alert an operator if for some unexpected reason conditions within the autoclave depart from the desired conditions or to shut down the entire system if the pressure within the autoclave approaches some predetermined safety limit.

One skilled in the art, would appreciate that the arrangement shown in FIG. 15 is merely illustrative of the many arrangements which could be employed to implement the program and flow chart shown in FIGS. 16-22. It must again be emphasized that the invention is not limited to the growth of quartz crystals, but has general applicability to the growth of any crystal which may grow hydrothermally. Such crystals include, but are not limited to, emerald, corundum (Al2 O3), in both its ruby and white sapphire form, berlinite (AlPO4) calcite (CaCO3), zincite (ZnO), tourmaline, magnetite, asbestos, fluorite, scheelite, garnets and zircons. See K. Nassau, "Growing Synthetic Crystals," Lapidary Journal, Vol. 18, Numbers 1-6, (1964).

Thus, for example, in the jewelry trade, synthetic rubys, sapphires and emeralds could be grown in a predetermined manner, so that the jewels exhibited unusual coloring or optical effects. In the communications and laser industry, laser crystals and optical modulators, beam splitters, etc. could be grown to exhibit predetermined optical and mechanical effects.

Indeed, the invention is of use in any industry or trade where crystals, of the type discussed, are employed.