JP2009190926A - シリコン単結晶製造における数値解析方法 - Google Patents

Abstract

【課題】シリコン単結晶の引上げ条件を変更しても、シリコン単結晶及びシリコン融液の固液界面形状やこの固液界面近傍の温度分布の計算値が実測値と極めて良く一致する、シリコン単結晶製造における数値解析方法を提供する。
【解決手段】メッシュ構造でモデル化したホットゾーンの各部材の物性値をコンピュータに入力し、各部材の表面温度分布をヒータ１９の発熱量及び各部材の輻射率に基づいて求める。上記各部材の表面温度分布及び熱伝導率に基づいて各部材の内部温度分布を求めた後に対流を考慮したシリコン融液１３の内部温度分布を求め、固液界面形状をシリコン単結晶２３の三重点を含む等温線に合せて求め、更に上記ステップを三重点がシリコン単結晶の融点になるまで繰返す。シリコン融液のメッシュを所定の範囲に限定し、乱流モデル式がｋｌ−乱流モデル式であり、このｋｌ−乱流モデル式中の乱流プラントル数がチューニングパラメータである。
【選択図】図１

Description

本発明は、チョクラルスキー（以下、ＣＺという。）法によりシリコン単結晶を引上げて製造するために、固体であるシリコン単結晶と液体であるシリコン融液との固液界面形状や、この固液界面近傍の温度分布を、コンピュータによりシミュレーションを行って数値解析を行う方法に関するものである。

従来、この種の数値解析方法として、メッシュ構造でモデル化したホットゾーンの各部材毎にまとめられたメッシュに対する各部材の物性値をそれぞれコンピュータに入力し、各部材の表面温度分布をヒータの発熱量及び各部材の輻射率に基づいて求め、各部材の表面温度分布及び熱伝導率に基づいて各部材の内部温度分布を求めた後に対流を考慮した融液の内部温度分布を更に求め、単結晶及び融液の固液界面形状を単結晶の三重点を含む等温線に合せて求め、上記ステップを三重点が単結晶の融点になるまで繰返す、単結晶及び融液の固液界面形状のシミュレーション方法が開示されている（例えば、特許文献１参照）。このシミュレーション方法では、融液のメッシュを所定の範囲に限定し、乱流モデル式が式（１）で表されるｋｌ−モデル式であり、このモデル式の乱流パラメータＣとして０．４〜０．６の範囲内の任意の値が用いられる。このように構成されたシミュレーション方法では、計算により得られた単結晶及び融液の固液界面形状は実測値と極めて良く一致するので、このシミュレーション方法で求められた固液界面形状はシリコン単結晶の引上げ時の点欠陥の拡散を考慮した結晶内分布を予測する計算の基礎とすることができるようになっている。
特許第３８４６１５５号公報（請求項１、段落［００２４］）

しかし、上記従来のシミュレーション方法では、シリコン単結晶の引上げ条件を変更すると、シリコン単結晶及びシリコン融液の固液界面形状やこの固液界面近傍の温度分布の計算値が実測値と一致しない場合があった。本発明の目的は、シリコン単結晶の引上げ条件を変更しても、シリコン単結晶及びシリコン融液の固液界面形状やこの固液界面近傍の温度分布の計算値が実測値と極めて良く一致する、シリコン単結晶製造における数値解析方法を提供することにある。

本発明者らは、シリコン単結晶及びシリコン融液の固液界面形状やこの固液界面近傍のシリコン単結晶の温度データ等の種々の実験結果から、乱流プラントル数を適切に定めることにより、上記固液界面形状や固液界面近傍の結晶温度分布を再現し予測できることを見いだし、本発明をなすに至った。

請求項１に係る発明は、コンピュータを用いて計算するシリコン単結晶製造の引上げ機のホットゾーン内のチャンバ、石英るつぼ、カーボンるつぼ、るつぼ支持柱、シリコン融液、シリコン単結晶、ヒートキャップ、保温筒、ヒータ、湯漏れ受け皿及び水冷体の各部材をメッシュ構造でモデル化する第１ステップと、上記ホットゾーンの各部材毎にメッシュをまとめかつこのまとめられたメッシュに対する上記各部材の熱伝導率、輻射率、粘性率、体積膨張係数、密度及び比熱からなる物性値をそれぞれ上記コンピュータに入力する第２ステップと、各部材の表面温度分布をヒータの発熱量及び各部材の輻射率に基づいて求める第３ステップと、各部材の表面温度分布及び熱伝導率に基づいて熱伝導方程式を解くことにより各部材の内部温度分布を求めた後にシリコン融液が乱流であると仮定して得られた乱流モデル式及びナビエ・ストークスの方程式を連結して解くことにより対流を考慮したシリコン融液の内部温度分布を更に求める第４ステップと、シリコン単結晶及びシリコン融液の固液界面形状をシリコン単結晶の三重点を含む等温線に合せて求める第５ステップと、上記第３ステップから第５ステップを三重点がシリコン単結晶の融点になるまで繰返す第６ステップとを含むシリコン単結晶を製造するためにコンピュータを用いて数値解析する方法の改良である。その特徴ある構成は、シリコン融液のメッシュのうちシリコン単結晶の径方向のメッシュであってかつシリコン融液のシリコン単結晶直下の一部又は全部のメッシュを０．０１〜５．００ｍｍに設定し、シリコン融液のメッシュのうちシリコン単結晶の長手方向のメッシュであってかつシリコン融液の一部又は全部のメッシュを０．０１〜５．００ｍｍに設定し、上記乱流モデル式がｋｌ−乱流モデル式であり、この乱流モデル式中の乱流プラントル数がチューニングパラメータであるところにある。

請求項２に係る発明は、請求項１に係る発明であって、乱流モデル式が次の式（１）で表されるｋｌ−乱流モデル式であり、この乱流モデル式の乱流パラメータＣとして０．４〜０．６の範囲内の任意の値が用いられかつ乱流モデル式中の乱流プラントル数Ｐｒ_tがチューニングパラメータであることを特徴とする。

ここで、κ_tはシリコン融液の乱流熱伝導率であり、ｃはシリコン融液の比熱であり、Ｐｒ_tは乱流プラントル数であり、ρはシリコン融液の密度であり、ｄはシリコン融液を貯留するるつぼ壁からの混合距離であり、ｋはシリコン融液の平均流速に対する変動成分の二乗和である。

本発明に係るコンピュータを用いて数値解析する方法では、シリコン単結晶の引上げ条件を種々変更しても、この数値解析方法により得られた単結晶及び融液の固液界面形状やこの固液界面近傍の温度分布が実測値と極めて良く一致する。この結果、本発明の数値解析方法で求めた固液界面形状及びこの固液界面近傍の温度分布を、単結晶の引上げ時の点欠陥の拡散を考慮した結晶内分布を精密に予測する計算の基礎とすることができる。

次に本発明を実施するための最良の形態を図面に基づいて説明する。図１に示すように、シリコン単結晶引上げ機１１のチャンバ１２内には、シリコン融液１３を貯留する石英るつぼ１４が設けられる。この石英るつぼ１４はカーボンるつぼ１６及びるつぼ支持柱１７を介してるつぼ駆動手段１８に接続され、るつぼ駆動手段１８は石英るつぼ１４を回転させるとともに昇降させるように構成される。また石英るつぼ１４の外周面は石英るつぼ１４から所定の間隔をあけてヒータ１９により包囲され、このヒータ１９は保温筒２１により包囲される。ヒータ１９は石英るつぼ１４に投入された高純度のシリコン多結晶体を加熱・溶融してシリコン融液１３にする。またチャンバ１２は、保温筒２１等を収容するチャンバ本体１２ａと、このチャンバ本体１２ａの上端の開口部に接続された円筒状のケーシング１２ｂとを有する。このケーシング１２ｂには引上げ手段２２が設けられる。この引上げ手段２２は、ケーシング１２ｂの上端部に設けられた引上げヘッド（図示せず）から石英るつぼ１４の回転中心に向って垂下されたワイヤケーブル２２ａと、このワイヤケーブル２２ａの下端に取付けられシリコン単結晶２３をシリコン融液１３から引上げるためにシリコン融液１３に浸される種結晶２２ｂとを有する。この引上げ手段２２はシリコン単結晶２３を回転させながらシリコン融液１３から引上げるように構成される。なお、図１の符号２４はシリコン融液１３から引上げているシリコン単結晶２３の外周面を包囲するヒートキャップである。このヒートキャップ２４は、下方に向かうに従って細くなるコーン状に形成されヒータ１９からの輻射熱を遮る熱遮蔽部２４ａと、この熱遮蔽部２４ａの上縁に連設され外方に略水平方向に張り出すフランジ部２４ｂと、このフランジ部２４ｂを保温筒２１上に載置するリング板２４ｃとを有する。上記熱遮蔽部２４ａはコーン状ではなく円筒状に形成してもよい。また、必要により単結晶２３の熱履歴を制御するための水冷体（図示せず）や、チャンバ１２の底面に設けられシリコン融液１３が石英るつぼ１４から漏れたときに処理するための湯漏れ受け皿（図示せず）が設置される。

このように構成された引上げ機１１におけるコンピュータを用いた数値解析方法、即ちシリコン単結晶２３及びシリコン融液１３の固液界面形状とこの固液界面近傍のシリコン単結晶２３内の温度分布とをコンピュータを用いてシミュレーションする方法を図１〜図３に基づいて説明する。先ず第１ステップとして引上げ機１１のホットゾーンの各部材、即ちチャンバ１２，石英るつぼ１４，カーボンるつぼ１６，るつぼ支持柱１７，シリコン融液１３，シリコン単結晶２３，ヒートキャップ２４，保温筒２１，ヒータ１９，湯漏れ受け皿，水冷体等をメッシュ分割してモデル化する。具体的には上記ホットゾーンの各部材のメッシュ点の座標データをコンピュータに入力する。このときシリコン融液１３のメッシュのうちシリコン単結晶２３の径方向のメッシュであってかつシリコン融液１３のシリコン単結晶２３直下の一部又は全部のメッシュ（以下、径方向メッシュという。）を０．０１〜５．００ｍｍ、好ましくは０．２５〜１．００ｍｍに設定する。またシリコン融液１３のメッシュのうちシリコン単結晶２３の長手方向のメッシュであってかつシリコン融液１３の一部又は全部のメッシュ（以下、長手方向メッシュという。）を０．０１〜５．００ｍｍ、好ましくは０．１〜０．５ｍｍに設定する。

径方向メッシュを０．０１〜５．００ｍｍの範囲に限定したのは、０．０１ｍｍ未満では計算時間が極めて長くなり、５．００ｍｍを越えると計算が不安定になり、繰返し計算を行っても固液界面形状が一定に定まらなくなるからである。また長手方向メッシュを０．０１〜５．００ｍｍの範囲に限定したのは、０．０１ｍｍ未満では計算時間が極めて長くなり、５．００ｍｍを越えると固液界面形状の計算値が実測値と一致しなくなるからである。なお、径方向メッシュの一部を０．０１〜５．００の範囲に限定する場合には、シリコン単結晶２３直下のシリコン融液１３のうちシリコン単結晶２３外周縁近傍のシリコン融液１３を上記範囲に限定することが好ましく、長手方向メッシュの一部を０．０１〜５．００の範囲に限定する場合には、シリコン融液１３の液面近傍及び底近傍を上記範囲に限定することが好ましい。

次いで第２ステップとして上記ホットゾーンの各部材毎にメッシュをまとめ、かつこのまとめられたメッシュに対して各部材の物性値をそれぞれコンピュータに入力する。例えば、チャンバ１２がステンレス鋼にて形成されていれば、そのステンレス鋼の熱伝導率，輻射率，粘性率，体積膨張係数，密度及び比熱がコンピュータに入力される。また後述する乱流モデル式（１）の乱流パラメータＣをコンピュータに入力する。

第３ステップとして、ホットゾーンの各部材の表面温度分布をヒータ１９の発熱量及び各部材の輻射率に基づいてコンピュータを用いて求める。即ち、ヒータ１９の発熱量を任意に設定してコンピュータに入力するとともに、各部材の輻射率から各部材の表面温度分布をコンピュータを用いて求める。次に第４ステップとしてホットゾーンの各部材の表面温度分布及び熱伝導率に基づいて熱伝導方程式（２）をコンピュータを用いて解くことにより各部材の内部温度分布を求める。ここでは、記述を簡単にするためｘｙｚ直交座標系を用いたが、実際の計算では円筒座標系を用いる。

ここで、ρは各部材の密度であり、ｃは各部材の比熱であり、Ｔは各部材の各メッシュ点での絶対温度であり、ｔは時間であり、λ_x，λ_y及びλ_zは各部材の熱伝導率のｘ，ｙ及びｚ方向成分であり、ｑはヒータ１９の発熱量である。

一方、シリコン融液１３に関しては、上記熱伝導方程式（２）でシリコン融液１３の内部温度分布を求めた後に、このシリコン融液１３の内部温度分布に基づき、シリコン融液１３が乱流であると仮定して得られた乱流モデル式（１）及びナビエ・ストークスの方程式（３）〜（５）を連結して、シリコン融液１３の内部流速分布をコンピュータを用いて求める。

ここで、κ_tはシリコン融液１３の乱流熱伝導率であり、ｃはシリコン融液１３の比熱であり、Ｐｒ_tは乱流プラントル数であり、ρはシリコン融液１３の密度であり、Ｃは乱流パラメータであり、ｄはシリコン融液１３を貯留する石英るつぼ１４壁からの混合距離であり、ｋはシリコン融液１３の平均流速に対する変動成分の二乗和である。なお、混合距離とは、乱流・対流において、運動量、エネルギ流束や拡散などの輸送現象をモデル化するために、分子運動論の類推からプラントルが導入したパラメータであり、渦塊が混合によりその個性を失うまでに流れを横切って進む有効的な長さをいう。

ここで、ｕ，ｖ及びｗはシリコン融液１３の各メッシュ点での流速のｘ，ｙ及びｚ方向成分であり、ν_lはシリコン融液１３の分子動粘性係数（物性値）であり、ν_tはシリコン融液１３の乱流の効果による動粘性係数であり、Ｆ_x，Ｆ_y及びＦ_zはシリコン融液１３に作用する体積力のｘ，ｙ及びｚ方向成分である。

上記乱流モデル式（１）はｋｌ（ケイエル）−乱流モデル式と呼ばれ、一般的には、壁近傍の粘性を考慮する係数ｆμを掛けた次の式（１−１）で表される。この式（１−１）において、上記壁近傍の粘性を考慮する係数ｆμを１とした場合が上記式（１）となる。一方、乱流粘性係数ｖ_tは次の式（１−２）で表されるので、式（１−１）を乱流粘性係数ｖ_tを用いて表すと、次の式（１−３）で表される。

上記式（１）に戻って、式（１）中の乱流パラメータＣは０．４〜０．６の範囲内の任意の値が用いられる。また上記式（１）中の乱流プラントル数Ｐｒ_tはチューニングパラメータであり、１未満、好ましくは０．３〜０．８、更に好ましくは０．４５〜０．６５の範囲に設定される。更に上記ナビエ・ストークスの方程式（３）〜（５）はシリコン融液１３が非圧縮性であって粘度が一定である流体としたときの運動方程式である。ここで、乱流パラメータＣを０．４〜０．６の範囲に限定したのは、０．４未満又は０．６を越えると計算により求めた界面形状が実測値と一致しないという不具合があるからである。また、乱流プラントル数Ｐｒ_tを１未満に設定したのは、１以上ではシリコン単結晶２３の引上げ条件を変更したときにシリコン単結晶２３及びシリコン融液１３の固液界面形状やこの固液界面近傍の温度分布の計算値が実測値と一致しないからである。なお、少なくともシリコン単結晶２３の直径が同一であり、ホットゾーンの形態も同一である場合には、引上げ速度、シリコン単結晶２３の回転速度、石英るつぼ１４の回転速度等が変更されても、定めた所定の乱流プラントル数Ｐｒ_tを用いて予測できる。またシリコン単結晶２３の直径やホットゾーンの形態が変更されても、乱流プラントル数Ｐｒ_tを上記好ましい範囲の０．４５〜０．６５に収斂することが種々の実験結果から分かり、これにより未だ実現されていないシリコン単結晶２３のサイズやホットゾーンの形態に対しても乱流プラントル数Ｐｒ_tを０．４５〜０．６５という狭い範囲のいずれかの値を用いることにより、大枠予測が可能である。この場合の乱流プラントル数Ｐｒ_t推奨値は０．５５である。

上記求められたシリコン融液１３の内部流速分布に基づいて熱エネルギ方程式（６）を解くことにより、シリコン融液１３の対流を考慮したシリコン融液１３の内部温度分布をコンピュータを用いて更に求める。

ここで、ｕ，ｖ及びｗはシリコン融液１３の各メッシュ点での流速のｘ，ｙ及びｚ方向成分であり、Ｔはシリコン融液１３の各メッシュ点での絶対温度であり、ρはシリコン融液１３の密度であり、ｃはシリコン融液１３の比熱であり、κ_lは分子熱伝導率（物性値）であり、κ_tは式（１）を用いて計算される乱流熱伝導率である。

次に第５ステップとして、シリコン単結晶２３及びシリコン融液１３の固液界面形状を図２の点Ｓで示すシリコンの三重点Ｓ（固体と液体と気体の三重点(tri-junction)）を含む等温線に合せてコンピュータを用いて求める。更にコンピュータに入力するヒータ１９の発熱量を変更して（次第に増大して）、上記第３ステップから第５ステップを三重点がシリコン単結晶２３の融点になるまで繰返す。このようにして得られたシリコン単結晶２３及びシリコン融液１３の固液界面形状と、固液界面近傍のシリコン単結晶２３内の温度分布は、実測値とほぼ一致する。この結果、本発明で求められた固液界面形状及び結晶内温度分布はシリコン単結晶２３の引上げ時の点欠陥の拡散を考慮した結晶内分布を精密に予測する計算の基礎とすることができる。

次に本発明の実施例を比較例とともに詳しく説明する。
＜実施例１＞
図１に示すように、石英るつぼ１４に４５ｋｇのシリコン原料をチャージしてこの原料をヒータ１９により融解し、このシリコン融液１３から直径１５０ｍｍのシリコン単結晶２３を速度０．８ｍｍ／分で引上げる場合の、シリコン単結晶２３及びシリコン融液１３の固液界面形状とこの固液界面近傍のシリコン単結晶２３内の温度分布を、図３のフローチャートに基づくコンピュータを用いた数値解析方法により求めた。即ち、シリコン単結晶引上げ機１１のホットゾーンをメッシュ構造でモデル化した。ここで、シリコン融液１３のシリコン単結晶２３直下のシリコン単結晶２３の径方向のメッシュを０．７５ｍｍに設定し、シリコン融液１３のシリコン単結晶２３直下以外のシリコン単結晶２３の径方向のメッシュを１〜５ｍｍに設定した。またシリコン融液１３のシリコン単結晶２３の長手方向のメッシュを０．２５〜５ｍｍに設定した。またｋｌ−乱流モデル式の乱流パラメータＣとして０．５４を用い、乱流プラントル数Ｐｒ_tを０．３２とした。更に種結晶２２ｂの回転速度を２０ｒｐｍとし、石英るつぼ１４の回転速度を８ｒｐｍとした。結晶長方向の中央である４００Ｌを代表長さとして、本発明の効果を確認した。
＜比較例１＞
ｋｌ−乱流モデル式の乱流プラントル数Ｐｒ_tを１としたこと以外は、実施例１と同様にしてシリコン単結晶及びシリコン融液の固液界面形状をコンピュータを用いた数値解析方法により求めた。

＜比較試験１及び評価＞
実施例１及び比較例１の数値解析方法により求めたシリコン単結晶及びシリコン融液の固液界面形状の実測値に対する偏差を図４に示す。なお、図４において、シリコン融液の表面を０．０００ｍｍとした。また、図４の横軸はシリコン単結晶の中心（０．０００ｍｍ）からシリコン単結晶の半径方向への距離である。図４から明らかなように、比較例１の数値解析方法で得られた固液界面形状（破線で示す。）は実測値（黒丸で示す。）と大幅に相違しているのに対し、実施例１の数値解析方法で得られた固液界面形状（実線で示す。）は実測値とほぼ一致していることが分かった。
＜比較試験２及び評価＞
実施例１及び比較例１の数値解析方法により求めた固液界面近傍のシリコン単結晶内の温度分布の実測値に対する偏差を図５に示す。なお、図５の横軸はシリコン単結晶内（シリコン単結晶の回転中心）におけるシリコン融液表面からの高さ（シリコン融液表面を０．０００ｍｍとした）である。図５から明らかなように、比較例１の数値解析方法で得られた固液界面近傍におけるシリコン単結晶内の温度分布（破線で示す。）は実測値（二点鎖線で示す。）と大幅に相違しているのに対し、実施例１の数値解析方法で得られた固液界面近傍におけるシリコン単結晶内の温度分布（実線で示す。）は実測値とほぼ一致していることが分かった。

＜実施例２＞
種結晶の回転速度を１０ｒｐｍとし、石英るつぼの回転速度を１２ｒｐｍとしたこと以外は、実施例１と同様にしてシリコン単結晶及びシリコン融液の固液界面形状をコンピュータを用いた数値解析方法により求めた。
＜比較例２＞
ｋｌ−乱流モデル式の乱流プラントル数Ｐｒ_tを１としたこと以外は、実施例２と同様にしてシリコン単結晶及びシリコン融液の固液界面形状をコンピュータを用いた数値解析方法により求めた。
＜比較試験３及び評価＞
実施例２及び比較例２の数値解析方法により求めたシリコン単結晶及びシリコン融液の固液界面形状の実測値に対する偏差を図６に示す。なお、図６において、シリコン融液の表面を０．０００ｍｍとした。また、図６の横軸はシリコン単結晶の中心（０．０００ｍｍ）からシリコン単結晶の半径方向への距離である。図６から明らかなように、比較例２の数値解析方法で得られた固液界面形状（破線で示す。）は実測値（黒三角で示す。）と大幅に相違しているのに対し、実施例２の数値解析方法で得られた固液界面形状（実線で示す。）は実測値とほぼ一致していることが分かった。

＜実施例３＞
シリコン融液から直径１５０ｍｍのシリコン単結晶を速度０．４ｍｍ／分で引上げ、ｋｌ−乱流モデル式の乱流プラントル数Ｐｒ_tを０．７６とし、種結晶の回転速度を１０ｒｐｍとし、石英るつぼの回転速度を１０ｒｐｍとしたこと以外は、実施例１と同様にしてシリコン単結晶及びシリコン融液の固液界面形状をコンピュータを用いた数値解析方法により求めた。
＜比較例３＞
ｋｌ−乱流モデル式の乱流プラントル数Ｐｒ_tを１としたこと以外は、実施例３と同様にしてシリコン単結晶及びシリコン融液の固液界面形状をコンピュータを用いた数値解析方法により求めた。
＜比較試験４及び評価＞
実施例３及び比較例３の数値解析方法により求めたシリコン単結晶及びシリコン融液の固液界面形状の実測値に対する偏差を図７に示す。なお、図７において、シリコン融液の表面を０．０００ｍｍとした。また、図７の横軸はシリコン単結晶の中心（０．０００ｍｍ）からシリコン単結晶の半径方向への距離である。図７から明らかなように、比較例３の数値解析方法で得られた固液界面形状（破線で示す。）は実測値（黒丸で示す。）と僅かに相違しているのに対し、実施例３の数値解析方法で得られた固液界面形状（実線で示す。）は実測値とほぼ一致していることが分かった。

＜実施例４＞
種結晶の回転速度を５ｒｐｍとしたこと以外は、実施例３と同様にしてシリコン単結晶及びシリコン融液の固液界面形状をコンピュータを用いた数値解析方法により求めた。
＜比較例４＞
ｋｌ−乱流モデル式の乱流プラントル数Ｐｒ_tを１としたこと以外は、実施例４と同様にしてシリコン単結晶及びシリコン融液の固液界面形状をコンピュータを用いた数値解析方法により求めた。
＜比較試験５及び評価＞
実施例４及び比較例４の数値解析方法により求めたシリコン単結晶及びシリコン融液の固液界面形状の実測値に対する偏差を図８に示す。なお、図８において、シリコン融液の表面を０．０００ｍｍとした。また、図８の横軸はシリコン単結晶の中心（０．０００ｍｍ）からシリコン単結晶の半径方向への距離である。図８から明らかなように、比較例４の数値解析方法で得られた固液界面形状（破線で示す。）は実測値（黒三角で示す。）と僅かに相違しているのに対し、実施例４の数値解析方法で得られた固液界面形状（実線で示す。）は実測値とほぼ一致していることが分かった。

＜実施例５＞
石英るつぼに９０ｋｇのシリコン原料をチャージしてこの原料をヒータにより融解し、このシリコン融液から直径２００ｍｍのシリコン単結晶を速度０．６ｍｍ／分で引上げ、ｋｌ−乱流モデル式の乱流プラントル数Ｐｒ_tを０．５５とし、種結晶の回転速度を１５ｒｐｍとし、石英るつぼの回転速度を１０ｒｐｍとしたこと以外は、実施例１と同様にしてシリコン単結晶及びシリコン融液の固液界面形状とこの固液界面近傍のシリコン単結晶内の温度分布を、コンピュータを用いた数値解析方法により求めた。
＜比較例５＞
ｋｌ−乱流モデル式の乱流プラントル数Ｐｒ_tを１としたこと以外は、実施例５と同様にしてシリコン単結晶及びシリコン融液の固液界面形状とこの固液界面近傍のシリコン単結晶２３内の温度分布を、コンピュータを用いた数値解析方法により求めた。

＜比較試験６及び評価＞
実施例５及び比較例５の数値解析方法により求めたシリコン単結晶及びシリコン融液の固液界面形状の実測値に対する偏差を図９に示す。なお、図９において、シリコン融液の表面を０．０００ｍｍとした。また、図９の横軸はシリコン単結晶の中心（０．０００ｍｍ）からシリコン単結晶の半径方向への距離である。図９から明らかなように、比較例５の数値解析方法で得られた固液界面形状（破線で示す。）は実測値（黒丸で示す。）と大幅に相違しているのに対し、実施例５の数値解析方法で得られた固液界面形状（実線で示す。）は実測値とほぼ一致していることが分かった。
＜比較試験７及び評価＞
実施例５及び比較例５の数値解析方法により求めた固液界面近傍のシリコン単結晶内の温度分布の実測値に対する偏差を図１０に示す。なお、図１０の横軸はシリコン単結晶内（シリコン単結晶の回転中心）におけるシリコン融液表面からの高さ（シリコン融液表面を０．０００ｍｍとした）である。図１０から明らかなように、比較例７の数値解析方法で得られた固液界面近傍におけるシリコン単結晶内の温度分布（破線で示す。）は実測値（二点鎖線で示す。）と大幅に相違しているのに対し、実施例７の数値解析方法で得られた固液界面近傍におけるシリコン単結晶内の温度分布（実線で示す。）は実測値とほぼ一致していることが分かった。

＜実施例６＞
種結晶の回転速度を７ｒｐｍとし、石英るつぼの回転速度を１２ｒｐｍとしたこと以外は、実施例５と同様にしてシリコン単結晶及びシリコン融液の固液界面形状をコンピュータを用いた数値解析方法により求めた。
＜比較例６＞
ｋｌ−乱流モデル式の乱流プラントル数Ｐｒ_tを１としたこと以外は、実施例６と同様にしてシリコン単結晶及びシリコン融液の固液界面形状をコンピュータを用いた数値解析方法により求めた。
＜比較試験８及び評価＞
実施例６及び比較例６の数値解析方法により求めたシリコン単結晶及びシリコン融液の固液界面形状の実測値に対する偏差を図１１に示す。なお、図１１において、シリコン融液の表面を０．０００ｍｍとした。また、図１１の横軸はシリコン単結晶の中心（０．０００ｍｍ）からシリコン単結晶の半径方向への距離である。図１１から明らかなように、比較例６の数値解析方法で得られた固液界面形状（破線で示す。）は実測値（黒三角で示す。）と大幅に相違しているのに対し、実施例６の数値解析方法で得られた固液界面形状（実線で示す。）は実測値とほぼ一致していることが分かった。

上述のことから、比較例では、シリコン単結晶の引上げ条件を変更すると、シリコン単結晶及びシリコン融液の固液界面形状やこの固液界面近傍の温度分布の計算値が実測値と大幅に相違してしまう場合があるのに対し、実施例では、シリコン単結晶の引上げ条件を変更しても、シリコン単結晶及びシリコン融液の固液界面形状やこの固液界面近傍の温度分布の計算値が実測値と極めて良く一致することが分かった。またシリコン単結晶２３の直径が同一であり、ホットゾーンの形態も同一である場合には、引上げ速度、シリコン単結晶の回転速度、石英るつぼの回転速度等が変更されても、定めた所定の乱流プラントル数Ｐｒ_tを用いて予測できることが分かった。

本発明実施形態のコンピュータを用いた数値解析方法に用いられる引上げ機の断面構成図である。 石英るつぼに貯留されるシリコン融液をメッシュ構造とした引上げ機の要部断面図である。 そのコンピュータを用いた数値解析方法を示すフローチャートである。 実施例１の方法で求めた固液界面形状と比較例１の方法で求めた固液界面形状とが、実際の固液界面形状から鉛直方向にどの程度偏っているかを示す図である。 実施例１の方法で求めた固液界面近傍における結晶内温度分布と比較例１の方法で求めた固液界面近傍における結晶内温度分布とが、実際の固液界面近傍における温度分布とどの程度相違しているかを示す図である。 実施例２の方法で求めた固液界面形状と比較例２の方法で求めた固液界面形状とが、実際の固液界面形状から鉛直方向にどの程度偏っているかを示す図である。 実施例３の方法で求めた固液界面形状と比較例３の方法で求めた固液界面形状とが、実際の固液界面形状から鉛直方向にどの程度偏っているかを示す図である。 実施例４の方法で求めた固液界面形状と比較例４の方法で求めた固液界面形状とが、実際の固液界面形状から鉛直方向にどの程度偏っているかを示す図である。 実施例５の方法で求めた固液界面形状と比較例５の方法で求めた固液界面形状とが、実際の固液界面形状から鉛直方向にどの程度偏っているかを示す図である。 実施例５の方法で求めた固液界面近傍における結晶内温度分布と比較例５の方法で求めた固液界面近傍における結晶内温度分布とが、実際の固液界面近傍における結晶内温度分布とどの程度相違しているかを示す図である。 実施例６の方法で求めた固液界面形状と比較例６の方法で求めた固液界面形状とが、実際の固液界面形状から鉛直方向にどの程度偏っているかを示す図である。

符号の説明

１１ シリコン単結晶引上げ機
１２ チャンバ
１３ シリコン融液
１４ 石英るつぼ
１６ カーボンるつぼ
１９ ヒータ
２１ 保温筒
２３ シリコン単結晶
２４ ヒートキャップ
Ｓ シリコンの三重点

Claims (5)

1. コンピュータを用いて計算するシリコン単結晶製造の引上げ機のホットゾーン内のチャンバ、石英るつぼ、カーボンるつぼ、るつぼ支持柱、シリコン融液、シリコン単結晶、ヒートキャップ、保温筒、ヒータ、湯漏れ受け皿及び水冷体の各部材をメッシュ構造でモデル化する第１ステップと、
   前記ホットゾーンの各部材毎にメッシュをまとめかつこのまとめられたメッシュに対する前記各部材の熱伝導率、輻射率、粘性率、体積膨張係数、密度及び比熱からなる物性値をそれぞれ前記コンピュータに入力する第２ステップと、
   前記各部材の表面温度分布をヒータの発熱量及び前記各部材の輻射率に基づいて求める第３ステップと、
   前記各部材の表面温度分布及び熱伝導率に基づいて熱伝導方程式を解くことにより前記各部材の内部温度分布を求めた後にシリコン融液が乱流であると仮定して得られた乱流モデル式及びナビエ・ストークスの方程式を連結して解くことにより対流を考慮した前記シリコン融液の内部温度分布を更に求める第４ステップと、
   前記シリコン単結晶及び前記シリコン融液の固液界面形状を前記シリコン単結晶の三重点を含む等温線に合せて求める第５ステップと、
   前記第３ステップから前記第５ステップを前記三重点が前記シリコン単結晶の融点になるまで繰返す第６ステップと
   を含むシリコン単結晶製造するためにコンピュータを用いて数値解析する方法において、
   前記シリコン融液のメッシュのうち前記シリコン単結晶の径方向のメッシュであってかつ前記シリコン融液の前記シリコン単結晶直下の一部又は全部のメッシュを０．０１〜５．００ｍｍに設定し、
   前記シリコン融液のメッシュのうち前記シリコン単結晶の長手方向のメッシュであってかつ前記シリコン融液の一部又は全部のメッシュを０．０１〜５．００ｍｍに設定し、
   前記乱流モデル式がｋｌ−乱流モデル式であり、この乱流モデル式中の乱流プラントル数がチューニングパラメータである
   ことを特徴とするシリコン単結晶製造における数値解析方法。
2. 乱流モデル式が次の式（１）で表されるｋｌ−乱流モデル式であり、この乱流モデル式の乱流パラメータＣとして０．４〜０．６の範囲内の任意の値が用いられかつ前記乱流モデル式中の乱流プラントル数Ｐｒ_tがチューニングパラメータである請求項１記載のシリコン単結晶製造における数値解析方法。
   

   

   ここで、κ_tはシリコン融液の乱流熱伝導率であり、ｃはシリコン融液の比熱であり、Ｐｒ_tは乱流プラントル数であり、ρはシリコン融液の密度であり、ｄはシリコン融液を貯留するるつぼ壁からの混合距離であり、ｋはシリコン融液の平均流速に対する変動成分の二乗和である。
3. 乱流プラントル数Ｐｒ_tが１．０未満である請求項２記載のシリコン単結晶製造における数値解析方法。
4. 乱流プラントル数Ｐｒ_tが０．３〜０．８である請求項２記載のシリコン単結晶製造における数値解析方法。
5. 乱流プラントル数Ｐｒ_tが０．４５〜０．６５である請求項２記載のシリコン単結晶製造における数値解析方法。

Images