JP4449347B2 - シミュレーションによるｏｓｆリングの分布予測方法 - Google Patents

Description

【０００１】
【発明の属する技術分野】
本発明は、チョクラルスキー（以下、ＣＺという。）法にて引上げられるシリコン等の単結晶内のグローイン欠陥の密度分布及びサイズ分布をコンピュータシミュレーションしてＯＳＦ（Oxidation Induced Stacking Fault：酸化誘起積層欠陥）リングの分布を予測する方法に関するものである。ここで言うグローイン欠陥とは、ボイドと内壁酸化膜からなるボイド欠陥のことであり、例えばウェーハ表面に観察される結晶起因のパーティクル（Crystal Originated Particle、以下、ＣＯＰという。）や、結晶内部に観察されるＬＳＴＤ（Laser Scattering Tomograph Defects）に対応する。
【０００２】
【従来の技術】
従来、この種のシミュレーション方法として、図８に示すように、総合伝熱シミュレータを用いてＣＺ法によるシリコン単結晶４引上げ時の引上げ機１内のホットゾーン構造及びそのシリコン単結晶４の引上げ速度に基づいて、シリコン融液２の熱伝導率を操作することによりシリコン融液２の内部温度分布を予測し、この内部温度分布からシリコン単結晶４のメッシュの座標及び温度をそれぞれ求め、更にシリコン単結晶４内の格子間シリコン及び空孔の拡散係数及び境界条件に基づいて拡散方程式を解くことにより、上記格子間シリコン及び空孔の密度分布をコンピュータを用いて求める方法が知られている。このシミュレーション方法では、ホットゾーンの各部材がメッシュ分割されてモデル化される。特にシリコン融液２のメッシュは計算時間を短くするために１０ｍｍ程度と比較的粗く設定される。
【０００３】
一方、ＣＺ法によるシリコン単結晶育成時の操業条件及び炉内の温度条件をコンピュータにデータとして入力し、単結晶内部に導入される空孔起因のグローンイン欠陥のサイズ及び密度と単結晶の成長速度との関係をコンピュータを用いたシミュレーションにより求め、このシミュレーション結果に基づいて単結晶内部に導入される空孔起因のグローンイン欠陥のサイズ及び密度が所定値になるように成長速度を選択して単結晶を育成するシリコン単結晶の製造方法が知られている（例えば、特許文献１参照）。このシリコン単結晶の製造方法では、総合伝熱解析プログラムを用いて炉内の熱分布のシミュレーションを行い、単結晶の育成時の単結晶中心の鉛直方向における温度勾配を求める。
【０００４】
また、抵抗率が０．０２Ωｃｍ以下であって、熱酸化処理をした際にリング状に発生するＯＳＦがウェーハ中心部で消滅したシリコンウェーハの表面にエピタキシャル層が形成されたエピタキシャルウェーハが知られている（例えば、特許文献２参照）。このエピタキシャルウェーハでは、ボロンのようなｐ型不純物を４×１０¹⁸ａｔｏｍｓ／ｃｍ³以上ドープし、熱酸化処理を行った際にＯＳＦリングがウェーハ中心部で消滅するようにＣＺ法によりシリコンウェーハを作製し、更にこのウェーハの表面にエピタキシャル層をＣＶＤ法にて形成することにより製造される。このように製造されたエピタキシャルウェーハでは、エピタキシャル層を形成した後の熱処理によりウェーハ面内で均一に高密度のＢＭＤを発生させてウェーハ面内で均一なＩＧ効果を得ることができるので、エピタキシャル層を形成したときにそのエピタキシャル層にＬ／Ｄの転写痕を生じさせないようになっている。
【０００５】
【特許文献１】
特開２００２−１４５６９６号公報
【特許文献２】
特開２０００−３３１９３３号公報
【０００６】
【発明が解決しようとする課題】
しかし、上記従来の格子間シリコン及び空孔の密度分布のシミュレーション方法では、実際の引上げ機においては発生するシリコン融液の対流を考慮しておらず、またシリコン融液のメッシュが比較的粗いため、固液界面形状の再現性が悪く、精度の良い単結晶内温度分布を提供できない。そのため、シリコン単結晶内の格子間シリコン及び空孔の密度分布が実測値と大幅に相違する問題点があった。また、上記従来のシミュレーション方法では、結晶冷却過程に比較的高温で生成されるシリコン単結晶内の欠陥、例えばボイドや酸素析出物の種類や、密度分布及びサイズ分布は判らなかった。
【０００７】
また、上記従来の特許文献１に示されたシリコン単結晶の製造方法では、シリコン単結晶の鉛直方向のグローンイン欠陥のサイズ分布及び密度分布を予測できるけれども、シリコン単結晶の半径方向のグローイン欠陥のサイズ分布及び密度分布を予測できず、このためＯＳＦリングの分布を予測できない問題点があった。
更に、上記従来の特許文献２に示されたエピタキシャルウェーハでは、不純物を高密度で添加した場合、添加密度に応じてＯＳＦリングの位置がシフトすることを確認できず、ＯＳＦリングの分布を予測することもできない問題点があった。
【０００８】
本発明の目的は、単結晶内でのＯＳＦリングの分布を正確に予測できる、シミュレーションによるＯＳＦリングの分布予測方法を提供することにある。
本発明の別の目的は、不純物（ドーパント）を添加した場合、この不純物の添加密度に応じてＯＳＦリングの分布を正確の予測できる、シミュレーションによるＯＳＦリングの分布予測方法を提供することにある。
【０００９】
【課題を解決するための手段】
請求項１に係る発明は、図１〜図７に示すように、引上げ機１１による単結晶１４の融液１２からの引上げ開始時から単結晶１４の冷却完了時までの引上げ機１１のホットゾーンをメッシュ構造でモデル化する第１ステップと、ホットゾーンの各部材毎にメッシュをまとめかつこのまとめられたメッシュに対する各部材の物性値とともに単結晶１４の引上げ長及びこの引上げ長に対応する単結晶１４の引上げ速度をそれぞれコンピュータに入力する第２ステップと、各部材の表面温度分布をヒータの発熱量及び各部材の輻射率に基づいて求める第３ステップと、各部材の表面温度分布及び熱伝導率に基づいて熱伝導方程式を解くことにより各部材の内部温度分布を求めた後に融液１２が乱流であると仮定して得られた乱流モデル式及びナビエ・ストークスの方程式を連結して解くことにより対流を考慮した融液１２の内部温度分布を更に求める第４ステップと、単結晶１４及び融液１２の固液界面形状を単結晶の三重点Ｓを含む等温線に合せて求める第５ステップと、第３ステップから第５ステップを三重点Ｓが単結晶１４の融点になるまで繰返し引上げ機１１内の温度分布を計算して単結晶１４のメッシュの座標及び温度を求めこれらのデータをそれぞれコンピュータに入力する第６ステップと、単結晶１４の引上げ長及び引上げ高さを段階的に変えて第１ステップから第６ステップまでを繰返し引上げ機１１内の温度分布を計算して単結晶１４のメッシュの座標及び温度を求めこれらのデータをそれぞれコンピュータに入力する第７ステップと、単結晶１４の融液１２からの引上げ開始時から単結晶１４の冷却完了時までの時間を所定の間隔毎に区切りこの区切られた時間間隔毎に第７ステップで求めた単結晶１４のメッシュの座標及び温度のデータから単結晶１４の引上げ長及び引上げ高さと単結晶１４内の温度分布とを求める第８ステップと、単結晶１４内の空孔及び格子間原子の拡散係数及び境界条件に基づいて拡散方程式を解くことにより所定の時間間隔の経過した後の空孔及び格子間原子の密度分布を求める第９ステップと、空孔の密度分布に基づいてボイド２１の形成開始温度を求める第１０ステップと、単結晶１４内のそれぞれのメッシュの格子点における温度が次第に低下してボイド２１の形成開始温度になったときのボイド２１の密度を求める第１１ステップと、単結晶１４内のそれぞれのメッシュの格子点における温度がボイド２１の形成開始温度より低いときのボイド２１の半径を、このボイド２１の成長速度が空孔の拡散速度に依存する拡散律速でありかつボイド２１の形状が球状であるとして求める第１２ステップと、ボイド２１の周囲に成長する内壁酸化膜２２の形成開始温度を求める第１３ステップと、単結晶１４内のそれぞれのメッシュの格子点における温度が次第に低下して内壁酸化膜２２の形成開始温度より低くなったときのボイド２１の半径と内壁酸化膜２２の厚さとを互いに関連させて求める第１４ステップと、第８ステップから第１４ステップを単結晶１４の冷却が完了するまで繰返す第１５ステップと、ボイド２１の密度分布及び半径と内壁酸化膜２２の厚さに基づいてＯＳＦリング２３のリング幅及び位置を求める第１６ステップとを含む、コンピュータを用いてシミュレーションによりＯＳＦリングの分布を予測する方法である。
【００１０】
この請求項１に記載されたＯＳＦリングの分布予測方法では、先ず融液１２の対流を考慮して融液１２から成長する単結晶１４内の温度分布を求めるだけでなく、更に冷却過程における単結晶１４内の温度分布までも求める、即ち融液１２から切離された単結晶１４の冷却過程における単結晶１４の徐冷及び急冷の効果を考慮して解析する。次に単結晶１４内のボイド２１の密度分布及び半径と内壁酸化膜２２の厚さを求めた後に、ボイド２１の密度分布及び半径と内壁酸化膜２２の厚さに基づいてＯＳＦリング２３のリング幅及び位置を求める。これによりＯＳＦリング２３の単結晶１４内での分布を正確に予測することができる。
【００１１】
請求項２に係る発明は、請求項１に係る発明であって、更に融液に添加する不純物の密度に基づき融液と単結晶の固液界面における空孔又は格子間原子の密度を変更して、請求項１に記載された第１〜第１６ステップを行う、コンピュータを用いてシミュレーションによりＯＳＦリングの分布を予測する方法である。
この請求項２に記載されたＯＳＦリングの分布予測方法では、融液に不純物を添加した場合、この不純物の添加密度に応じてＯＳＦリングの分布を正確に予測することができる。
【００１２】
【発明の実施の形態】
次に本発明の実施の形態を図面に基づいて説明する。
図５に示すように、シリコン単結晶引上げ機１１のチャンバ内には、シリコン融液１２を貯留する石英るつぼ１３が設けられる。この石英るつぼ１３は図示しないが黒鉛サセプタ及び支軸を介してるつぼ駆動手段に接続され、るつぼ駆動手段は石英るつぼ１３を回転させるとともに昇降させるように構成される。また石英るつぼ１３の外周面は石英るつぼ１３から所定の間隔をあけてヒータ（図示せず）により包囲され、このヒータは保温筒（図示せず）により包囲される。ヒータは石英るつぼ１３に投入された高純度のシリコン多結晶体を加熱・溶融してシリコン融液１２にする。またチャンバの上端には図示しないが円筒状のケーシングが接続され、このケーシングには引上げ手段が設けられる。引上げ手段はシリコン単結晶１４を回転させながら引上げるように構成される。
【００１３】
このように構成されたシリコン単結晶引上げ機１１におけるシリコン単結晶１４内のボイド欠陥の密度分布及びサイズ分布のシミュレーション方法を図１〜図７に基づいて説明する。
先ず第１ステップとして、シリコン単結晶１４を所定長さＬ₁（例えば１００ｍｍ）まで引上げた状態におけるシリコン単結晶引上げ機１１のホットゾーンの各部材、即ちチャンバ，石英るつぼ１３，シリコン融液１２，シリコン単結晶１４，黒鉛サセプタ，保温筒等をメッシュ分割してモデル化する。具体的には上記ホットゾーンの各部材のメッシュ点の座標データをコンピュータに入力する。このときシリコン融液１２のメッシュのうちシリコン単結晶１４の径方向のメッシュであってかつシリコン融液１２のシリコン単結晶１４直下の一部又は全部のメッシュ（以下、径方向メッシュという。）を０．０１〜５．００ｍｍ、好ましくは０．２５〜１．００ｍｍに設定する。またシリコン融液１２のメッシュのうちシリコン単結晶１４の長手方向のメッシュであってかつシリコン融液１２の一部又は全部のメッシュ（以下、長手方向メッシュという。）を０．０１〜５．００ｍｍ、好ましくは０．１〜０．５ｍｍに設定する。
【００１４】
径方向メッシュを０．０１〜５．００ｍｍの範囲に限定したのは、０．０１ｍｍ未満では計算時間が極めて長くなり、５．００ｍｍを越えると計算が不安定になり、繰返し計算を行っても固液界面形状が一定に定まらなくなるからである。また長手方向メッシュを０．０１〜５．００ｍｍの範囲に限定したのは、０．０１ｍｍ未満では計算時間が極めて長くなり、５．００ｍｍを越えると固液界面形状の計算値が実測値と一致しなくなるからである。なお、径方向メッシュの一部を０．０１〜５．００の範囲に限定する場合には、シリコン単結晶１４直下のシリコン融液１２のうちシリコン単結晶１４外周縁近傍のシリコン融液１２を上記範囲に限定することが好ましく、長手方向メッシュの一部を０．０１〜５．００の範囲に限定する場合には、シリコン融液１２の液面近傍及び底近傍を上記範囲に限定することが好ましい。
【００１５】
第２ステップとして上記ホットゾーンの各部材毎にメッシュをまとめ、かつこのまとめられたメッシュに対して各部材の物性値をそれぞれコンピュータに入力する。例えば、チャンバがステンレス鋼にて形成されていれば、そのステンレス鋼の熱伝導率，輻射率，粘性率，体積膨張係数，密度及び比熱がコンピュータに入力される。またシリコン単結晶１４の引上げ長及びこの引上げ長に対応するシリコン単結晶１４の引上げ速度と、後述する乱流モデル式（１）の乱流パラメータＣとをコンピュータに入力する。
【００１６】
第３ステップとして、ホットゾーンの各部材の表面温度分布をヒータの発熱量及び各部材の輻射率に基づいてコンピュータを用いて求める。即ち、ヒータの発熱量を任意に設定してコンピュータに入力するとともに、各部材の輻射率から各部材の表面温度分布をコンピュータを用いて求める。次に第４ステップとしてホットゾーンの各部材の表面温度分布及び熱伝導率に基づいて熱伝導方程式（１）をコンピュータを用いて解くことにより各部材の内部温度分布を求める。ここでは、記述を簡単にするためｘｙｚ直交座標系を用いたが、実際の計算では円筒座標系を用いる。
【数１】

【００１７】
上記式（１）において、ρは各部材の密度、ｃは各部材の比熱、Ｔは各部材の各メッシュ点での絶対温度、ｔは時間である。またλ_x，λ_y及びλ_zは各部材の熱伝導率のｘ，ｙ及びｚ方向成分であり、ｑはヒータの発熱量である。
【００１８】
一方、シリコン融液１２に関しては、上記熱伝導方程式（１）でシリコン融液１２の内部温度分布を求めた後に、このシリコン融液１２の内部温度分布に基づき、シリコン融液１２が乱流であると仮定して得られた乱流モデル式（２）及びナビエ・ストークスの方程式（３）〜（５）を連結して、シリコン融液１２の内部流速分布をコンピュータを用いて求める。
【数２】

【００１９】
上記式（２）において、κ_tはシリコン融液１２の乱流熱伝導率であり、ｃはシリコン融液１２の比熱であり、Ｐｒ_tはプラントル数であり、ρはシリコン融液１２の密度であり、Ｃは乱流パラメータであり、ｄはシリコン融液１２を貯留する石英るつぼ１３壁からの距離であり、ｋはシリコン融液１２の平均流速に対する変動成分の二乗和である。
【数３】

【００２０】
上記式（３）〜式（５）において、ｕ，ｖ及びｗはシリコン融液１２の各メッシュ点での流速のｘ，ｙ及びｚ方向成分であり、ν_lはシリコン融液１２の分子動粘性係数（物性値）であり、ν_tはシリコン融液１２の乱流の効果による動粘性係数であり、Ｆ_x，Ｆ_y及びＦ_zはシリコン融液１２に作用する体積力のｘ，ｙ及びｚ方向成分である。
【００２１】
上記乱流モデル式（２）はｋｌ（ケイエル）−モデル式と呼ばれ、このモデル式の乱流パラメータＣは０．４〜０．６の範囲内の任意の値が用いられることが好ましい。乱流パラメータＣを０．４〜０．６の範囲に限定したのは、０．４未満又は０．６を越えると計算により求めた界面形状が実測値と一致しないという不具合があるからである。また上記ナビエ・ストークスの方程式（３）〜（５）はシリコン融液１２が非圧縮性であって粘度が一定である流体としたときの運動方程式である。
【００２２】
上記求められたシリコン融液１２の内部流速分布に基づいて熱エネルギ方程式（６）を解くことにより、シリコン融液１２の対流を考慮したシリコン融液１２の内部温度分布をコンピュータを用いて更に求める。
【数４】

【００２３】
上記式（６）において、ｕ，ｖ及びｗはシリコン融液１２の各メッシュ点での流速のｘ，ｙ及びｚ方向成分であり、Ｔはシリコン融液１２の各メッシュ点での絶対温度であり、ρはシリコン融液１２の密度であり、ｃはシリコン融液１２の比熱であり、κ_lは分子熱伝導率（物性値）であり、κ_tは式（１）を用いて計算される乱流熱伝導率である。
【００２４】
次いで第５ステップとして、シリコン単結晶１４及びシリコン融液１２の固液界面形状を図５の点Ｓで示すシリコンの三重点Ｓ（固体と液体と気体の三重点(tri-junction)）を含む等温線に合せてコンピュータを用いて求める。第６ステップとして、コンピュータに入力するヒータの発熱量を変更し（次第に増大し）、上記第３ステップから第５ステップを三重点がシリコン単結晶１４の融点になるまで繰返した後に、引上げ機１１内の温度分布を計算してシリコン単結晶のメッシュの座標及び温度を求め、これらのデータをコンピュータに記憶させる。
【００２５】
次に第７ステップとして、シリコン単結晶１４の引上げ長Ｌ₁にδ（例えば５０ｍｍ）だけ加えて上記第１ステップから第６ステップまでを繰返した後に、引上げ機１１内の温度分布を計算してシリコン単結晶１４のメッシュの座標及び温度を求め、これらのデータをコンピュータに記憶させる。この第７ステップは、シリコン単結晶１４の引上げ長Ｌ₁が長さＬ₂（Ｌ₂はシリコン融液１２から切離されたときのシリコン単結晶１４の長さ（成長完了時の結晶長）である。）に達してシリコン単結晶１４がシリコン融液１２から切離された後、更にシリコン単結晶１４が引上げられてその高さＨ₁（Ｈ₁はシリコン単結晶１４の直胴開始部からシリコン融液１２の液面までの距離である（図５）。）がＨ₂（Ｈ₂は冷却完了時のシリコン単結晶１４の直胴開始部からシリコン融液１２の液面までの距離である。）に達するまで、即ちシリコン単結晶１４の冷却が完了するまで行われる。なお、シリコン単結晶１４がシリコン融液１２から切離された後は、シリコン単結晶１４の引上げ高さＨ₁にδ（例えば５０ｍｍ）だけ加え、上記と同様に上記第１ステップから第６ステップまでを繰返す。
【００２６】
シリコン単結晶１４の引上げ高さＨ₁がＨ₂に達すると、第８ステップに移行する。第８ステップでは、シリコン単結晶１４をシリコン融液１２から成長させて引上げ始めたときｔ₀から、シリコン単結晶１４をシリコン融液１２から切離して更にシリコン単結晶１４を引上げ、その冷却が完了したときｔ₁までの時間を、所定の間隔Δｔ秒（微小時間間隔）毎に区切る。このときシリコン単結晶１４内の格子間シリコン及び空孔の拡散係数及び境界条件のみならず、後述するボイド２１及び内壁酸化膜２２（図６（ａ））の密度分布及びサイズ分布を求めるための式に用いられる定数をそれぞれコンピュータに入力する。上記区切られた時間間隔Δｔ秒毎に、第７ステップで求めたシリコン単結晶１４のメッシュの座標及び温度のデータから、シリコン単結晶１４の引上げ長Ｌ₁及び引上げ高さＨ₁と、シリコン単結晶１４内の温度分布とを求める。
【００２７】
即ち、第１〜第７ステップでシリコン単結晶のメッシュの座標及び温度を引上げ長δ毎に求め、シリコン単結晶を例えば５０ｍｍ引上げるのに数十分要するため、この数十分間でのシリコン単結晶のメッシュの温度変化を時間の関数として微分することにより、時刻ｔ₀からΔｔ秒後におけるシリコン単結晶１４の引上げ長Ｌ₁及び引上げ高さＨ₁とシリコン単結晶１４内の温度分布を算出する。次にシリコン単結晶１４内の空孔及び格子間シリコンの拡散係数及び境界条件に基づいて拡散方程式を解くことにより、Δｔ秒経過後の空孔及び格子間シリコンの密度分布を求める（第９ステップ）。
【００２８】
具体的には、空孔の密度Ｃ_vの計算式が次の式（７）で示され、格子間シリコンの密度Ｃ_iの計算式が次の式（８）で示される。式（７）及び式（８）において、密度Ｃ_v及び密度Ｃ_iの経時的進展を計算するために、空孔と格子間シリコンの熱平衡がシリコン単結晶の全表面で維持されると仮定する。
【数５】

【００２９】
上記式（７）及び式（８）において、Ｋ₁及びＫ₂は定数であり、Ｅ_i及びＥ_vはそれぞれ格子間シリコン及び空孔の形成エネルギーであり、Ｃ_ve及びＣ_ieは空孔の平衡密度及び格子間シリコンの平衡密度である。またｋ_Bはボルツマン定数、Ｔは絶対温度を意味する。
【００３０】
上記平衡式は時間で微分され、空孔及び格子間シリコンに対してそれぞれ次の式（９）及び式（１０）になる。
【数６】

【００３１】
上記式（９）及び式（１０）において、Θ(ｘ)はヘビサイド関数（Heaviside function）である。即ち、ｘ＜０のときΘ(ｘ)＝０であり、かつｘ＞０のときΘ(ｘ)＝１である。またＴ_pは内壁酸化膜２２（ボイド２１表面に形成されたシリコンの酸化膜（ＳｉＯｘ膜））の形成開始温度であり、Ｔ_vはボイド２１の形成開始温度である。更に式（９）及び（１０）のそれぞれ右側の第１項はフィックの拡散式であり、右側の第１項中のＤ_v及びＤ_iは、次の式（１１）及び（１２）で表される空孔及び格子間シリコンの拡散係数である。
【数７】

【００３２】
上記式（１１）及び式（１２）において、△Ｅ_v及び△Ｅ_iはそれぞれ空孔及び格子間シリコンの活性化エネルギーであり、ｄ_v及びｄ_iはそれぞれ定数である。また式（９）及び式（１０）のそれぞれ右側の第２項中のＥ_v ^t及びＥ_i ^tは熱拡散による空孔及び格子間シリコンの活性化エネルギーであり、ｋ_Bはボルツマン定数である。式（９）及び式（１０）のそれぞれ右側の第３項のｋ_ivは空孔及び格子間シリコンペアの再結合定数である。式（９）及び式（１０）のそれぞれ右側の第４項のＮ_vはボイド２１の密度であり、ｒ_vはボイド２１の半径であり、更に式（９）の右側の第５項のＮ_pは内壁酸化膜２２の密度であり、Ｒ_pは内壁酸化膜２２の外半径である。
【００３３】
上記式（９）が成立つのは、空孔が析出するための空孔の流束が十分に大きく、シリコン単結晶１４を構成するＳｉマトリックスと内壁酸化膜を構成するＳｉＯｘとの単位質量当りの体積差を埋められる場合、即ちγＤ_OＣ_O≦Ｄ_v（Ｃ_v−Ｃ_ve）の場合である。上記以外の場合、即ちγＤ_OＣ_O＞Ｄ_v（Ｃ_v−Ｃ_ve）の場合には、次式（１３）が成立つ。ここで、γとは酸素１原子に対する空孔の消費割合であり、Ｄ_Oは酸素の拡散係数であり、Ｃ_Oは酸素密度である。
【数８】

【００３４】
上記式（１３）において、Θ(ｘ)はヘビサイド関数（Heaviside function）である。即ち、ｘ＜０のときΘ(ｘ)＝０であり、かつｘ＞０のときΘ(ｘ)＝１である。またＴ_pは内壁酸化膜２２の形成開始温度であり、Ｔ_vはボイド２１の形成開始温度である。更に式（１３）の右側の第１項はフィックの拡散式であり、右側の第１項中のＤ_vは、式（１１）で表される空孔の拡散係数である。
【００３５】
次に第１０ステップとして、上記拡散方程式を解くことにより求めた空孔の密度Ｃ_v分布に基づいて、ボイド２１の形成開始温度Ｔ_vを次の式（１４）から求める。
【数９】

【００３６】
上記式（１４）において、Ｃ_vはシリコン単結晶１４中の空孔密度であり、Ｃ_v0はシリコン単結晶１４の融点Ｔ_mでの空孔平衡密度であり、Ｅ_vは空孔形成エネルギーである。またσ_vはシリコン単結晶１４の結晶面(１１１)における界面エネルギーであり、ρはシリコン単結晶１４の密度であり、ｋ_Bはボルツマン定数である。
【００３７】
第１１ステップとして、シリコン単結晶１４内のそれぞれのメッシュの格子点における温度が次第に低下してボイド２１の形成開始温度Ｔ_vになったときに、次の近似式（１５）を用いてボイド２１の密度Ｎ_vを求める。
【数１０】

【００３８】
上記式（１５）において、ｎ_vはボイド２１の密度であり、（dT/dt）はシリコン単結晶１４の冷却速度である。上記ボイド２１の密度は空孔過飽和度に依存するけれども、計算するシリコン単結晶１４での核形成温度のような狭い温度領域では一定とみなしてよい。またＤ_vはボイド２１の拡散係数であり、ｋ_Bはボルツマン定数であり、Ｃ_vはシリコン単結晶１４中の空孔密度である。
【００３９】
第１２ステップとして、シリコン単結晶１４内のそれぞれのメッシュの格子点における温度がボイド２１の形成開始温度Ｔ_vより低いときのボイドの半径ｒ_vを、次の式（１６）から求める。ここでボイド２１の成長速度は空孔の拡散速度に依存する拡散律速である。またボイド２１の形状は実際には八面体であるけれども、ここでは計算の効率化から球状として扱う。
【数１１】

【００４０】
上記式（１６）において、ｔ１はシリコン単結晶１４のメッシュの格子点における温度がボイド２１の形成開始温度Ｔ_vまで低下したときの時刻である。またｒ_vcrはボイド２１の臨界径であり、この臨界径ｒ_vcrは上記時刻ｔ１での値とする。更にＣ_ve及びＣ_ieはそれぞれの状態での空孔及び格子間シリコンの平衡密度である。
【００４１】
次に第１３ステップとして、ボイド２１の周囲に成長する内壁酸化膜２２の形成開始温度Ｔ_pを求める。内壁酸化膜２２の形成開始温度Ｔ_pがボイド２１の形成開始温度Ｔ_vより小さい場合には、酸素と空孔が結合し、ボイド２１表面に内壁酸化膜２２が成長する（図６（ａ））。本発明のモデルでは、上記内壁酸化膜２２はボイド２１が発生するとすぐに成長するものとして扱った。ここで、「すぐに」とは、「ボイドが発生した次の繰返し計算から」という意味である。従って、内壁酸化膜２２の形成開始温度Ｔ_pはボイド２１の形成開始温度Ｔ_vと同一となる。
【００４２】
更に第１４ステップとして、シリコン単結晶１４内のそれぞれのメッシュの格子点における温度が次第に低下して内壁酸化膜２２の形成開始温度Ｔ_pより低くなったときのボイド２１の半径ｒ_vと内壁酸化膜２２の厚さｄとを互いに関連させて求める。一般的に、酸素はその拡散とシリコンとの反応によって内壁酸化膜２２が形成される。酸素が内壁酸化膜２２外からこの酸化膜の外周面に流入する流束Ｊ_Oは、ボイド２１の中心を原点とする内壁酸化膜２２の外半径をＲ_pとするとき、次の式（１７）から求められる。
【数１２】

【００４３】
上記式（１７）において、Ｄ_Oは酸素の拡散係数であり、Ｃ_Oは酸素密度であり、Ｃ_Oeifは球状のボイド２１の中心を原点とする半径Ｒ_pにおける酸素平衡密度であり、ｇ_OはＳｉＯｘ及びＳｉの界面、即ち内壁酸化膜２２の外周面における酸素原子の反応割合である。
【００４４】
また空孔が内壁酸化膜２２外からこの酸化膜の外周面（半径Ｒ_pの球面）に流入する流束Ｊ_vは次の式（１８）から求められる。
【数１３】

【００４５】
上記式（１８）において、ｇ_vは内壁酸化膜２２の外周面での空孔の反応割合であり、Ｃ_veifは内壁酸化膜２２の外周面での空孔の平衡密度である。シリコンの内壁酸化膜２２が成長すると同時に空孔が消費され、空孔の消費割合は酸素１原子に対してγである。このため内壁酸化膜２２が歪みなく成長するための条件は、酸素の流束Ｊ_Oに対して空孔の流束はγＪ_Oであるので、空孔の流束γＪ_Oは次の式（１９）から求められる。
【数１４】

【００４６】
従って、(Ｊ_v−γＪ_O)は内壁酸化膜２２に吸収される空孔の流束となる。ここで空孔の内壁酸化膜２２外からこの酸化膜への流束が増大する場合((Ｊ_v−γＪ_O)＞０)、或いは空孔のボイド２１内から酸化膜への流束が時間の経過に対して増大する場合((Ｊ_v−γＪ_O)＜０)の上記流束の時間に対する変化の割合をαとすると、このαは内壁酸化膜２２の厚さに依存すると考えられる。このため内壁酸化膜２２の厚さｄ、即ち(Ｒ_p−ｒ_v)がｄ₀未満であるならば、内壁酸化膜２２は空孔を酸化膜外及びボイド２１内の双方向から部分浸透させると考えられる。換言すれば、ｄ＜ｄ₀であってα≠０である場合には、空孔はある割合で内壁酸化膜２２を貫通してボイド２１の成長に使われる。但し、格子間シリコンは酸素析出層が核形成されると同時にボイド２１の成長に使われなくなる。
【００４７】
また空孔の内壁酸化膜２２外からこの酸化膜への流束量が少ない場合、即ち(Ｊ_v−γＪ_O)＜０である場合には、歪みの無い内壁酸化膜２２を成長させることができないので、空孔はボイド２１内から酸化膜に流入し、ボイド２１の外径は小さくなる。この結果、ｄ≧ｄ₀であってα＝０の場合には、内壁酸化膜２２は空孔を全く通過させなくなる。そして、n_vをボイド２１中の空孔の数とすると、dn_v/dtは空孔の流入量となり、dn_v/dt＝α(Ｊ_v−γＪ_O)であるならば、次の式（２０）が成り立つ。なお、式（２０）において、ｒ_vは式（１６）で求めたボイド２１の半径であり、Ｒ_pは内壁酸化膜２２の外半径である。またＤ_Oは酸素の拡散係数であり、Ｃ_Oは酸素密度である。更にＣ_Oeifは球状のボイド２１の中心を原点とする半径Ｒ_pにおける酸素平衡密度であり、ｇ_Oは内壁酸化膜２２の外周面における酸素原子の反応割合である。
【数１５】

【００４８】
同時に内壁酸化膜２２を構成するＳｉＯｘ分子の数は式（２０）の酸素が内壁酸化膜２２外からこの酸化膜の外周面に流入する流束Ｊ_Oの消費に依存して変化するため、次の式（２１）が成り立つ。なお、式（２１）において、ρ_p＝ｘ／（ΩＳｉＯｘ）と定義する。このρ_pの定義式において、Ω_SiOxはＳｉＯｘの分子容量（molecular volume）であり、その単位は分子量／密度である。またｘはＳｉＯｘのｘに対応し、Ωは定数である。
【数１６】

【００４９】
上記式（２１）を上記式（２０）に代入すると次の式（２２）が得られる。なお、式（２２）において、ηはΩ_SiOx／Ω_Siである。このηの定義式において、Ω_Siは１／ρであり、ρはシリコン単結晶１４の密度である。
【数１７】

【００５０】
上記式（２０）及び式（２２）から次のことが分かる。先ずｄがｄ₀になった時刻をｔ₃とすると、この時刻ｔ₃においてαがゼロとなるため、内壁酸化膜２２は非浸透となる。またボイド２１の外径がｒ_v(ｔ₃)で固定されている間は、内壁酸化膜２２の成長はシリコン単結晶１４のマトリックスから流入する空孔或いは酸素の消費量によって決まる。即ち、（Ｊ_v−γＪ_O）＞０である場合には式（２２）は次の式（２３）となり、（Ｊ_v−γＪ_O）＜０である場合には式（２２）は次の式（２４）となる。
【数１８】

【００５１】
上記式（２３）及び式（２４）において、Ω_SiOxはＳｉＯｘの分子容量（molecular volume）である。また式（２３）及び式（２４）において、内壁酸化膜の外周面での空孔平衡密度はＲ_pと無関係であり、Ｃ_veif＝Ｃ_veかつＣ_Oeif＝Ｃ_Oeであると想定している。式（２３）又は式（２４）から内壁酸化膜２２の外半径Ｒ_pを求め、内壁酸化膜２２の外半径Ｒ_pからボイド２１の半径ｒ_vを引くことにより、内壁酸化膜２２の厚さｄを求める。上記第８ステップから第１４ステップをシリコン単結晶１４の冷却が完了するまで繰返す（第１５ステップ）。なお、式（９）〜式（２４）は互いに関連させてコンピュータにより解く。
【００５２】
一方、シリコン単結晶１４内で空孔型点欠陥が優勢であって凝集した点欠陥を有する領域［Ｖ］と、格子間シリコン型点欠陥が優勢であって凝集した点欠陥を有する領域［Ｉ］との境界付近では、内壁酸化膜２２の厚さｄがボイド２１の半径と比較して極めて大きい領域がＯＳＦリングの位置にほぼ一致していることが分かっている。この内壁酸化膜２２の厚さｄがボイド２１の半径と比較して極めて大きい領域では、ボイド２１は殆ど成長しない（図６（ｃ）のＯＳＦ核のイメージ図）。これにより図７に示すように、上記ボイドの密度分布と、ボイドの半径ｒ_vと、内壁酸化膜の厚さｄとに基づいて、ＯＳＦリング２３のリング幅と、ＯＳＦリング２３の発生する位置とをコンピュータを用いて求めることができる。なお、図６（ｂ）はボイド２１の半径ｒ_vが内壁酸化膜２２の厚さより極めて大きいＣＯＰのイメージ図を示す。
【００５３】
上述のように、シリコン融液１２の対流を考慮してシリコン融液１２から成長するシリコン単結晶１４内の温度分布を求めた後に、シリコン融液１２から切離されたシリコン単結晶１４の冷却過程、即ちシリコン単結晶１４がシリコン融液１２から切離された後のシリコン単結晶１４の引上げ速度を考慮し、シリコン単結晶１４の徐冷及び急冷の効果を結果に反映して解析することにより、シリコン単結晶１４内におけるＯＳＦリング２３の分布を正確に予測できる。即ち、ボイド２１の密度分布をコンピュータを用いて求めた後に、ボイド２１の半径と、ボイド２１の周囲に成長する内壁酸化膜２２の厚さとを互いに関連させてコンピュータを用いて求め、更にボイド２１の密度分布と、ボイド２１の半径と、内壁酸化膜２２の厚さに基づいて、ＯＳＦリング２３のリング幅と、ＯＳＦリング２３のシリコン単結晶１４内における位置とをコンピュータを用いて求めることにより、シリコン単結晶１４内におけるＯＳＦリング２３の分布を正確に予測できる。この結果、シリコン融液１２から引上げられるシリコン単結晶１４内のＯＳＦリング２３を所望の分布にするために、引上げ機１１の設計段階で引上げ機１１の構造を検討することができる。
【００５４】
なお、上記実施の形態では、不純物（ドーパント）の添加を考慮しなかったが、シリコン融液に不純物を添加する場合には、シリコン融液に添加する不純物の密度に基づきシリコン単結晶内のシリコン融点における空孔又は格子間シリコンの密度を変更して、上記実施の形態の第１〜第１６ステップを行う。具体的には、ｒを不純物原子の半径とし、ｒ_Siをシリコン原子の半径とするとき、ｒ＞ｒ_Siの場合には、式（２５）に基づいて空孔の密度Ｃ_VDを求め、ｒ＜ｒ_Siの場合には、式（２６）に基づいて格子間シリコンの密度Ｃ_IDを求める。
【数１９】

【００５５】
上記式（２５）及び式（２６）において、（Ｃ_V）₀及び（Ｃ_I）₀は不純物を添加しない場合の空孔及び格子間シリコンの熱平衡密度であり、上記実施の形態の式（７）及び式（８）のＣ_ve及びＣ_ieに相当する。また式（２５）及び式（２６）において、ａは調整可能なパラメータであり、Ｎは不純物の密度であり、ｋ_Bはボルツマン定数であり、Ｔは絶対温度である。なお、ｒ＜ｒ_Siである不純物としては、ボロン（Ｂ）やリン（Ｐ）などが挙げられ、ｒ＞ｒ_Siである不純物としては、アンチモン（Ｓｂ）やヒ素（Ａｓ）などが挙げられる。
【００５６】
上記式（２５）及び式（２６）は、シリコン融液とシリコン単結晶の固液界面におけるシリコン単結晶内の空孔及び格子間シリコンの熱平衡密度をそれぞれ変更する式であり、初期値として上記固液界面での空孔の密度Ｃ_VD及び格子間シリコンの密度Ｃ_IDを上記式（２５）及び式（２６）から求め、その値をコンピュータに入力した後は、上記実施の形態と同一のプロセスでコンピュータを用いて計算する。これによりシリコン融液に不純物を添加した場合における不純物の添加密度に応じたＯＳＦリングの分布を正確に予測できる。
また、上記実施の形態では、単結晶としてシリコン単結晶を挙げたが、ＧａＡｓ単結晶，ＩｎＰ単結晶，ＺｎＳ単結晶若しくはＺｎＳｅ単結晶でもよい。
【００５７】
【実施例】
次に本発明の実施例を比較例とともに詳しく説明する。
＜実施例１＞
図５及び図６に示すように、石英るつぼ１３に貯留されたシリコン融液１２から直径６インチのシリコン単結晶１４を引上げる場合の、シリコン単結晶１４内のＯＳＦリングのリング幅及び位置を、図１〜図４のフローチャートに基づくシミュレーション方法により求めた。なお、上記シリコン融液１２には、不純物（ドーパント）としてボロン（Ｂ）を添加した。このボロンは、シリコン融液１２から引上げられるシリコン単結晶１４のトップ部の抵抗率及びボロン密度がそれぞれ０．０１５Ωｃｍ及び４．８９×１０¹⁸／ｃｍ³となるように添加した。
【００５８】
具体的には、先ず、シリコン単結晶引上げ機１１のホットゾーンをメッシュ構造でモデル化した。ここで、シリコン融液１２のシリコン単結晶１４直下のシリコン単結晶１４の径方向のメッシュを０．７５ｍｍに設定し、シリコン融液１２のシリコン単結晶１４直下以外のシリコン単結晶１４の径方向のメッシュを１〜５ｍｍに設定した。またシリコン融液１２のシリコン単結晶１４の長手方向のメッシュを０．２５〜５ｍｍに設定し、乱流モデル式の乱流パラメータＣとして０．４５を用いた。更にシリコン単結晶１４の引上げ開始時ｔ₀から冷却完了時ｔ₁までの引上げ長及び引上げ高さの段階的な変更を５０ｍｍずつとした。このような条件下で、シリコン融液１２の対流を考慮してシリコン単結晶１４内の温度分布を求めた。
【００５９】
次にシリコン単結晶１４の引上げ開始時ｔ₀から冷却完了時ｔ₁までの時間を所定の間隔Δｔ秒（微小な時間の間隔）毎に区切り、この時間間隔Δｔ秒毎に上記シリコン単結晶１４のメッシュの座標及び温度のデータから、シリコン単結晶１４の引上げ長及び引上げ高さとシリコン単結晶内１４の温度分布を求め、更にボイド２１の形成開始温度を求めてボイド２１の密度分布及びサイズ分布を求めた。即ち、シリコン単結晶１４をシリコン融液１２から切離した後のシリコン単結晶１４の徐冷及び急冷を考慮して、シリコン単結晶１４内のボイド２１の密度分布及びサイズ分布をコンピュータを用いてそれぞれ求めた。
【００６０】
一方、ボイド２１の周囲に成長する内壁酸化膜２２は、ボイド２１が発生するとすぐに成長するものとして扱ったので、内壁酸化膜２２の形成開始温度Ｔ_pをボイド２１の形成開始温度Ｔ_vと同一とした。ここで、「すぐに」とは、「ボイドが発生した次の繰返し計算から」という意味である。次に内壁酸化膜２２の形成開始温度より低くなったときの内壁酸化膜２２の外半径Ｒ_pをコンピュータを用いて求めた。更に内壁酸化膜２２の外半径Ｒ_pをボイド２１の半径から引くことにより内壁酸化膜２２の厚さｄを求めた。なお、内壁酸化膜の密度分布は上記ボイドの密度分布と同一とした。更に上記ボイド２１の密度分布及び半径と内壁酸化膜２２の厚さに基づいてＯＳＦリングのリング幅及び位置をコンピュータを用いて求めた。このシリコン単結晶を実施例１とした。
【００６１】
なお、上記計算において、ボロン原子の半径ｒはシリコン原子ｒ_Siの半径より小さい、即ちｒ＜ｒ_Siであるため、ボロンが添加されたシリコン融液とこのシリコン融液から引上げられるシリコン単結晶との固液界面における空孔の密度Ｃ_VDを式（２５）から求め、この値を初期値としてコンピュータに入力した。ここで、式（２５）中のａには、４．０×１０²ｅＶを代入した。
また上記計算において、シリコン単結晶１４中の酸素濃度を１．２×１０¹⁸〜１．３×１０¹⁸／ｃｍ³（旧ＡＳＴＭ）に徐々に増大させ、シリコン単結晶の引上げ速度を１．３〜１．４ｍｍ／分に徐々に増加させた。
【００６２】
＜実施例２＞
シリコン単結晶のトップ部の抵抗率及びボロン密度がそれぞれ０．０１６Ωｃｍ及び４．６０×１０¹⁸／ｃｍ³となるように、ボロンをシリコン単結晶に添加したことを除いて、実施例１と同様にして、ＯＳＦリングのリング幅及び位置をコンピュータを用いて求めた。このシリコン単結晶を実施例２とした。
＜実施例３＞
シリコン単結晶のトップ部の抵抗率及びボロン密度がそれぞれ０．０１７Ωｃｍ及び４．１０×１０¹⁸／ｃｍ³となるように、ボロンをシリコン単結晶に添加したことを除いて、実施例１と同様にして、ＯＳＦリングのリング幅及び位置をコンピュータを用いて求めた。このシリコン単結晶を実施例３とした。
【００６３】
＜比較例１＞
実施例１と同一形状のシリコン単結晶を実施例１と同一の条件で実際に引上げて、所定の熱処理を行ってＯＳＦリングを顕在化させた。このシリコン単結晶を比較例１とした。
＜比較例２＞
実施例２と同一形状のシリコン単結晶を実施例２と同一の条件で実際に引上げて、所定の熱処理を行ってＯＳＦリングを顕在化させた。このシリコン単結晶を比較例２とした。
＜比較例３＞
実施例３と同一形状のシリコン単結晶を実施例３と同一の条件で実際に引上げて、所定の熱処理を行ってＯＳＦリングを顕在化させた。このシリコン単結晶を比較例３とした。
【００６４】
＜比較試験１及び評価＞
実施例１のシリコン単結晶をトップ部から１２００ｍｍの位置で、引上げ方向に対して直角にスライスして得られたウェーハの半径方向におけるＣＯＰの半径、ボイドの半径及び内壁酸化膜の厚さの分布をコンピュータを用いて求めた。その結果を図９に示す。また比較例１のシリコン単結晶をトップ部から１２００ｍｍの位置で、引上げ方向に対して直角にスライスして得られたウェーハのＯＳＦリングのリング幅及び中心位置を実際に観察して測定した。その結果を図９に示す。
図９から明らかなように、実施例１のＣＯＰ発生領域からウェーハの半径方向外側に向うとき、内壁酸化膜の厚さがボイドの半径より大きくなる位置（以下、変節点という）が存在するけれども、この位置が比較例１のＯＳＦリングと一致することが分かった。
【００６５】
＜比較試験２及び評価＞
実施例１〜３のシリコン単結晶を引上げ方向に所定の間隔をあけかつ引上げ方向に対して直角にそれぞれスライスして得られた各ウェーハの上記変節点をコンピュータを用いて求めた。その結果を図１０（ａ）〜図１０（ｃ）に示す。また比較例１のシリコン単結晶を、上記と同様に、引上げ方向に所定の間隔をあけかつ引上げ方向に対して直角にそれぞれスライスして得られた各ウェーハのＯＳＦリングの中心位置を実際に観察して測定した。その結果を図１０（ａ）〜図１０（ｃ）に示す。
図１０（ａ）〜図１０（ｃ）から明らかなように、実施例１〜３の変節点が比較例１〜３のＯＳＦリングに略一致することが分かった。
【００６６】
【発明の効果】
以上述べたように、本発明によれば、融液の対流を考慮して融液から成長する単結晶内の温度分布をコンピュータを用いて求めるだけでなく、更に冷却過程における単結晶内の温度分布までも求めることによって、即ち融液から切離された単結晶の冷却過程における単結晶の徐冷及び急冷の効果を考慮することによって、ボイドの密度分布をコンピュータを用いて求めた後に、ボイドの半径と、ボイドの周囲に成長する内壁酸化膜の厚さとを互いに関連させてコンピュータを用いて求め、更にボイドの密度分布及び半径と内壁酸化膜の厚さに基づいてＯＳＦリングのリング幅及び位置をコンピュータを用いて求める。この結果、単結晶内でのＯＳＦリングの分布を正確に予測できる。
また融液に添加する不純物の密度に基づき融液と単結晶の固液界面における空孔又は格子間原子の密度を変更して、上記第１〜第１６ステップを行えば、融液に不純物を添加した場合、この不純物の添加密度に応じてＯＳＦリングの分布を正確に予測できる。
【図面の簡単な説明】
【図１】本発明実施形態シミュレーションによるＯＳＦリングの分布予測方法の第１段を示すフローチャート。
【図２】そのシミュレーションによるＯＳＦリングの分布予測方法の第２段を示すフローチャート。
【図３】そのシミュレーションによるＯＳＦリングの分布予測方法の第３段を示すフローチャート。
【図４】そのシミュレーションによるＯＳＦリングの分布予測方法の第４段を示すフローチャート。
【図５】本発明のシリコン融液をメッシュ構造としたシリコン単結晶の引上げ機の要部断面図。
【図６】そのシリコン単結晶内のボイド及び内壁酸化膜の模式図。
【図７】そのシミュレーションにより予測したＯＳＦリングの位置を示すシリコン単結晶の縦断面図。
【図８】従来例のシリコン融液をメッシュ構造としたシリコン単結晶の引上げ機の要部断面図。
【図９】実施例１のシリコン単結晶の半径方向の位置によって、ＣＯＰの半径の変化、ボイドの半径の変化及び内壁酸化膜の半径の変化をそれぞれ示す図。
【図１０】実施例１〜３のシリコン単結晶内のＯＳＦリングの位置をシミュレーションにより求めた値と、比較例１〜３のシリコン単結晶内のＯＳＦリングの位置を実測した値とを比較した図。
【符号の説明】
１１ シリコン単結晶引上げ機
１２ シリコン融液
１４ シリコン単結晶
２１ ボイド
２２ 内壁酸化膜
２３ ＯＳＦリング
Ｓ シリコンの三重点

Claims (2)

1. 引上げ機(11)による単結晶(14)の融液(12)からの引上げ開始時から前記単結晶(14)の冷却完了時までの前記引上げ機(11)のホットゾーンをメッシュ構造でモデル化する第１ステップと、
   前記ホットゾーンの各部材毎にメッシュをまとめかつこのまとめられたメッシュに対する前記各部材の物性値とともに前記単結晶(14)の引上げ長及びこの引上げ長に対応する前記単結晶(14)の引上げ速度をそれぞれコンピュータに入力する第２ステップと、
   前記各部材の表面温度分布をヒータの発熱量及び前記各部材の輻射率に基づいて求める第３ステップと、
   前記各部材の表面温度分布及び熱伝導率に基づいて熱伝導方程式を解くことにより前記各部材の内部温度分布を求めた後に融液(12)が乱流であると仮定して得られた乱流モデル式及びナビエ・ストークスの方程式を連結して解くことにより対流を考慮した前記融液(12)の内部温度分布を更に求める第４ステップと、
   前記単結晶(14)及び前記融液(12)の固液界面形状を前記単結晶の三重点(S)を含む等温線に合せて求める第５ステップと、
   前記第３ステップから前記第５ステップを前記三重点(S)が前記単結晶(14)の融点になるまで繰返し前記引上げ機(11)内の温度分布を計算して前記単結晶(14)のメッシュの座標及び温度を求めこれらのデータをそれぞれ前記コンピュータに入力する第６ステップと、
   前記単結晶(14)の引上げ長及び引上げ高さを段階的に変えて前記第１ステップから前記第６ステップまでを繰返し前記引上げ機(11)内の温度分布を計算して前記単結晶(14)のメッシュの座標及び温度を求めこれらのデータをそれぞれ前記コンピュータに入力する第７ステップと、
   前記単結晶(14)の前記融液(12)からの引上げ開始時から前記単結晶(14)の冷却完了時までの時間を所定の間隔毎に区切り前記区切られた時間間隔毎に第７ステップで求めた前記単結晶(14)のメッシュの座標及び温度のデータから前記単結晶(14)の引上げ長及び引上げ高さと前記単結晶(14)内の温度分布とを求める第８ステップと、
   前記単結晶(14)内の空孔及び格子間原子の拡散係数及び境界条件に基づいて拡散方程式を解くことにより前記所定の時間間隔の経過した後の空孔及び格子間原子の密度分布を求める第９ステップと、
   前記空孔の密度分布に基づいてボイド(21)の形成開始温度を求める第１０ステップと、
   前記単結晶(14)内のそれぞれのメッシュの格子点における温度が次第に低下して前記ボイド(21)の形成開始温度になったときの前記ボイド(21)の密度を求める第１１ステップと、
   前記単結晶(14)内のそれぞれのメッシュの格子点における温度が前記ボイド(21)の形成開始温度より低いときの前記ボイド(21)の半径を、このボイド(21)の成長速度が空孔の拡散速度に依存する拡散律速でありかつ前記ボイド(21)の形状が球状であるとして求める第１２ステップと、
   前記ボイド(21)の周囲に成長する内壁酸化膜(22)の形成開始温度を求める第１３ステップと、
   前記単結晶(14)内のそれぞれのメッシュの格子点における温度が次第に低下して前記内壁酸化膜(22)の形成開始温度より低くなったときの前記ボイド(21)の半径と前記内壁酸化膜(22)の厚さとを互いに関連させて求める第１４ステップと、
   第８ステップから第１４ステップを前記単結晶(14)の冷却が完了するまで繰返す第１５ステップと、
   前記ボイド(21)の密度分布及び半径と前記内壁酸化膜(22)の厚さに基づいてＯＳＦリング(23)のリング幅及び位置を求める第１６ステップと
   を含むコンピュータを用いてシミュレーションによりＯＳＦリングの分布を予測する方法。
2. 融液に添加する不純物の密度に基づき前記融液と単結晶の固液界面における空孔又は格子間原子の密度を変更して、請求項１に記載された第１〜第１６ステップを行う、コンピュータを用いてシミュレーションによりＯＳＦリングの分布を予測する方法。

Images