(57)【要約】【課題】本発明は、物体に加える熱負荷の位置を時間とともに移動することにより、物体に発生する変位分布を有限要素法を用いて解析する有限要素解析方法および装置等に関し、解析時間の短縮化を損なうことなく精度を向上することを目的とする。【解決手段】物体に加える熱負荷の位置を時間とともに移動することにより、前記物体に発生する変位分布を有限要素法を用いて解析する有限要素解析方法において、前記物体をシェル要素で分割するとともに、前記熱負荷による温度分布を求め、前記温度分布から前記シェル要素の中立面変位分布を求め、前記物体の厚み方向任意位置での変位分布を、前記中立面変位と前記シェル要素の厚み方向の位置との関数で近似し、前記物体に発生する変位分布を求めることを特徴とする。【選択図】      図４
2005025588.tif 000002
【特許請求の範囲】
【請求項１】
物体に加える熱負荷の位置を時間とともに移動することにより、前記物体に発生する変位分布を有限要素法を用いて解析する有限要素解析方法において、
前記物体をシェル要素で分割するとともに、前記熱負荷による温度分布を求め、前記温度分布から前記シェル要素の中立面変位分布を求め、前記物体の厚み方向任意位置での変位分布を、前記中立面変位と前記シェル要素の厚み方向の位置との関数で近似し、前記物体に発生する変位分布を求めることを特徴とする有限要素解析方法。
【請求項２】
請求項１記載の有限要素解析方法において、
前記熱負荷の少なくとも１つの位置について、前記シェル要素の面積を、前記シェル要素の位置と前記熱負荷の位置との距離の関数として変化するようにする処理を有することを特徴とする有限要素解析方法。
【請求項３】
請求項１または請求項２記載の有限要素解析方法において、前記物体に対向して第２の物体が位置され、前記物体に属する節点と前記第２の物体に属する節点との間の剛性方程式を、相対変位増分の変化に対する力の変化を表す折れ線状のグラフによって与える処理を有することを特徴とする有限要素解析方法。
【請求項４】
物体に加える熱負荷の位置を時間とともに移動することにより、前記物体に発生する変位分布を有限要素法を用いて解析する有限要素法解析装置において、
前記物体をシェル要素で分割する要素分割手段と、
前記熱負荷による前記物体の温度分布を求める温度分布算出手段と、
前記温度分布から前記シェル要素の中立面変位分布を求め、前記物体の厚み方向任意位置での変位分布を、前記中立面変位と前記シェル要素の厚み方向の位置との関数で近似し、前記物体に発生する変位分布を求める変位分布算出手段と、を有することを特徴とする有限要素法解析装置。
【請求項５】
請求項４記載の有限要素法解析装置において、
前記熱負荷の少なくとも１つの位置について、前記シェル要素の面積を、前記シェル要素の位置と前記熱負荷の位置との距離の関数として変化するようにする面積変更手段を有することを特徴とする有限要素法解析装置。
【請求項６】
請求項４または請求項５記載の有限要素法解析装置において、
前記物体に対向して第２の物体が位置され、前記物体に属する節点と前記第２の物体に属する節点との間の剛性方程式を、相対変位増分の変化に対する力の変化を表す折れ線状のグラフによって与える剛性・内力算出手段を有することを特徴とする有限要素法解析装置。
【請求項７】
請求項１ないし請求項３のいずれか１項記載の有限要素解析方法による解析結果を用いて製造されたことを特徴とする露光装置。
【請求項８】
請求項４に記載された変位分布算出手段として機能するプログラムが格納され、かつ、コンピュータ読み取り可能な記録媒体。
【請求項９】
請求項５に記載された面積変更手段として機能するプログラムが格納され、かつ、コンピュータ読み取り可能な記録媒体。
【請求項１０】
請求項６に記載された剛性・内力算出手段として機能するプログラムが格納され、かつ、コンピュータ読み取り可能な記録媒体。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、物体に加える熱負荷の位置を時間とともに移動することにより、物体に発生する変位分布を有限要素法を用いて解析する有限要素解析方法および装置に関する。また、この方法による解析結果を用いて製造された露光装置およびこの方法を実行させるためのプログラムを記録したコンピュータ読み取り可能な記録媒体に関する。
【０００２】
【従来の技術】
従来、物体に加える熱負荷の位置を時間とともに移動することにより、物体に発生する変位分布を有限要素法を用いて解析する有限要素解析方法および装置として、例えば、特開２００１−３３７９４０号公報に開示されるものが知られている。
【０００３】
この公報に開示される有限要素解析方法では、互いに対向する第１物体および第２物体に加える負荷の位置が、時間とともに移動され、この移動にともなって第１物体と第２物体とに発生する各時間における第１物体と第２物体との変形が有限要素法を用いて解析される。
そして、第１物体および第２物体を、２次元の要素および節点に要素分割するようにしたので、物体を３次元の要素として要素分割する方法に比較して、非常に短時間で有限要素解析を行うことが可能になる。
【０００４】
【特許文献１】
特開２００１−３３７９４０号公報
【０００５】
【発明が解決しようとする課題】
しかしながら、上述した有限要素解析方法では、第１物体および第２物体を、２次元の要素および節点に要素分割、すなわち、物体を２次元のプレーン要素として要素分割しているため、物体を３次元のソリッド要素として要素分割する方法に比較して、有限要素解析の精度が低下するという問題があった。
【０００６】
本発明は、かかる従来の問題を解決するためになされたもので、解析時間の短縮化を損なうことなく精度を向上することができる有限要素解析方法および装置、ならびに、この方法による解析結果を用いて製造された露光装置およびこの方法をコンピュータに実行させるプログラムを記録したコンピュータ読み取り可能な記録媒体を提供することを目的とする。
【０００７】
【課題を解決するための手段】
請求項１の有限要素解析方法は、物体に加える熱負荷の位置を時間とともに移動することにより、前記物体に発生する変位分布を有限要素法を用いて解析する有限要素解析方法において、前記物体をシェル要素で分割するとともに、前記熱負荷による温度分布を求め、前記温度分布から前記シェル要素の中立面変位分布を求め、前記物体の厚み方向任意位置での変位分布を、前記中立面変位と前記シェル要素の厚み方向の位置との関数で近似し、前記物体に発生する変位分布を求めることを特徴とする。
【０００８】
請求項２の有限要素解析方法は、請求項１記載の有限要素解析方法において、前記熱負荷の少なくとも１つの位置について、前記シェル要素の面積を、前記シェル要素の位置と前記熱負荷の位置との距離の関数として変化するようにする処理を有することを特徴とする。
請求項３の有限要素解析方法は、請求項１または請求項２記載の有限要素解析方法において、前記物体に対向して第２の物体が位置され、前記物体に属する節点と前記第２の物体に属する節点との間の剛性方程式を、相対変位増分の変化に対する力の変化を表す折れ線状のグラフによって与える処理を有することを特徴とする。
【０００９】
請求項４の有限要素解析装置は、物体に加える熱負荷の位置を時間とともに移動することにより、前記物体に発生する変位分布を有限要素法を用いて解析する有限要素法解析装置において、前記物体をシェル要素で分割する要素分割手段と、前記熱負荷による前記物体の温度分布を求める温度分布算出手段と、前記温度分布から前記シェル要素の中立面変位分布を求め、前記物体の厚み方向任意位置での変位分布を、前記中立面変位と前記シェル要素の厚み方向の位置との関数で近似し、前記物体に発生する変位分布を求める変位分布算出手段とを有することを特徴とする。
【００１０】
請求項５の有限要素法解析装置は、請求項４記載の有限要素法解析装置において、前記熱負荷の少なくとも１つの位置について、前記シェル要素の面積を、前記シェル要素の位置と前記熱負荷の位置との距離の関数として変化するようにする面積変更手段を有することを特徴とする。
【００１１】
請求項６の有限要素法解析装置は、請求項４または請求項５記載の有限要素法解析装置において、前記物体に対向して第２の物体が位置され、前記物体に属する節点と前記第２の物体に属する節点との間の剛性方程式を、相対変位増分の変化に対する力の変化を表す折れ線状のグラフによって与える剛性・内力算出手段を有することを特徴とする。
【００１２】
請求項７の露光装置は、請求項１ないし請求項３のいずれか１項記載の有限要素解析方法による解析結果を用いて製造されたことを特徴とする。
請求項８のコンピュータ読み取り可能な記録媒体は、請求項４に記載された変位分布算出手段として機能するプログラムが格納されていることを特徴とする。
請求項９のコンピュータ読み取り可能な記録媒体は、請求項５に記載された面積変更手段として機能するプログラムが格納されていることを特徴とする。
【００１３】
請求項１０のコンピュータ読み取り可能な記録媒体は、請求項６に記載された剛性・内力算出手段として機能するプログラムが格納されていることを特徴とする。
【００１４】
【発明の実施の形態】
以下、本発明の実施形態を、図面を用いて詳細に説明する。
【００１５】
（第１の実施形態）
図１は、本実施形態における有限要素解析装置のブロック図である。
この実施形態の解析装置１００は、フレキシブルディスク装置（以下、「ＦＤ装置」と略記する。）１１１、入力装置１１２、出力装置１１３、フレキシブルディスク制御部（以下、「ＦＤ制御部」と略記する。）１１４、入力制御部１１５、出力制御部１１６、バス１１７、中央処理装置１１８および記憶装置１１９を備えて構成される。
【００１６】
ＦＤ装置１１１は、ＦＤ制御部１１４によって制御され、記録媒体の１つであるフレキシブルディスク（以下、「ＦＤ」と略記する。）１０１に記録された本発明にかかるプログラムを解析装置１００にインストールしたり、解析結果をこのＦＤ１０１または別のＦＤに記録したりする。
入力装置１１２は、キーボード、マウスおよびデジタイザなどの入力機器であり、入力制御部１１５によって制御される。ユーザは、入力装置１１２を操作してコマンドや解析対象の形状などのデータを入力する。
【００１７】
出力装置１１３は、ＣＲＴなどのディスプレイやプリンタなどの出力機器であり、出力制御部１１６によって制御される。出力装置１１３は、入力装置１１２に入力されたデータやコマンド、プログラム実行のためのメニューおよび解析結果などが出力される。
記憶装置１１９は、半導体メモリやハードディスなどであり、本発明にかかるプログラムおよび各プログラム実行中の一時的なデータや解析結果などが格納される。
【００１８】
中央処理装置１１８は、バス１１７を介して、ＦＤ制御部１１４、入力制御部１１５、出力制御部１１６および記憶装置１１９と接続され、これらに指示を与えたり、記憶装置１１９に格納された各プログラムに従い処理を行う。
図２は、本実施形態の評価モデルを示す図である。
図２（ａ）に示すように、本実施形態の評価モデルである物体１１Ａは、矩形状の半導体ウェハである。
【００１９】
この物体１１Ａの中央の領域には、図２（ｂ）に示すように、熱負荷領域が６×２０のマトリックスで構成されている。各マトリックスは、２５０μｍ角の大きさを有している。
このマトリックスに対し、熱負荷を矢印の経路に沿って順次与える。この評価モデルの場合は、１２０回である。なお、この評価モデルでは、物体１１Ａに対して無限境界を設定している。
【００２０】
図３は、本実施形態の解析装置の動作を示すフローチャートである。
先ず、ユーザは、ステップＳ１において、解析装置１００を入力装置１１２から起動する。
次に、ステップＳ２において、中央処理装置１１８は、解析対象の形状・寸法などの初期データを入力するようにユーザに出力装置１１３を介して指示する。
【００２１】
ユーザは、入力装置１１２から、解析対象の形状・寸法のデータ、解析対象の物性のデータ（熱伝導率、比熱、密度、熱膨張率、ヤング率、ポアソン比）、熱負荷のデータ（熱エネルギー、位置、位置の移動軌跡、移動時間、負荷面積）、解析対象の初期温度、解析対象の外界の温度などを入力する。
なお、予めこれらのデータを記憶装置１１９またはＦＤ１０１に格納しておき、解析装置１００は、起動に従って出力装置１１３に格納内容を表示して、ユーザは、これから選択して各データを入力するようにしてもよい。
【００２２】
例えば、解析対象の形状のデータとして、円盤や方形盤などのテーブルを用意する。物性のデータとして、シリコン、３−５族化合物半導体  、および、鋼やアルミ合金などの金属材料などの各物性値のテーブルを用意する。熱負荷のデータとして、値別の熱エネルギや異なる移動順に設けられた複数個のパターンなどのデータを用意する。
【００２３】
次に、ステップＳ３において、中央処理装置１１８は、入力された解析対象の形状・寸法のデータを基に、解析対象をモデル化する。
すなわち、図４の（ａ）に示すように、中央処理装置１１８は、物体１１Ａを解析対象部分について正方形の要素に分割する。そして、各要素に節点Ｎｉを割り当て、座標と初期温度と初期変位φ（ゼロ）を割り振る。
【００２４】
本実施形態では、物体１１Ａを、厚肉シェル要素で分割する。
一般に、構造解析に用いられる有限要素には、大きく連続体要素と構造体要素とがある。連続体要素は、ソリッド要素ともよばれ、要素を構成する節点は、最低でも要素の頂点の数だけあり、その形は、立体形状となる。また、各節点が、並進自由度３成分を持つことで変形を表す。
【００２５】
これに対し、構造体要素は、その構造上の特徴から変形しないと仮定できる個所を定式化で考慮し、要素形状を省略した要素で、要素形状が板または棒となる。この際、省略された個所に相当する形状因子を定式化で考慮するため、要素を構成する各節点は、並進自由度の他、回転自由度を有する。
シェル要素とは、上記の構造要素にあたり比較的に厚さの薄い構造物に対し、その厚さが変化しないと仮定された要素で、その要素の定式化で厚さが考慮され、形状は、厚さ方向を省略した板形状となる。このシェル要素の定式化には、多種あるが、比較的厚みが厚く、その剪断的挙動を考慮した要素を厚肉シェル要素という。この実施形態で言えば、図２に示す様に、Ｚ方向を厚さ方向とし、ＸＹ面内を要素分割する。
【００２６】
次に、ステップＳ４において、中央処理装置１１８は、記憶装置１１９に格納されている汎用有限要素法プログラム、例えば、エムエスシーソフトウェア株式会社の「ＭＳＣ．Ｍａｒｃ」（登録商標）を呼び出し、この解析用モデルを汎用有限要素法プログラムに適用して、要素１に熱負荷が加えられた初期の場合における温度分布および変位分布を解析する。
【００２７】
すなわち、図４の（ａ）に示すように、加えられた熱負荷Ｑによる厚肉シェル要素の厚み方向の温度は、節点Ｎｉの持つ３つの自由度（θｉ１，θｉ２，θｉ３）として求められる。自由度を持つ位置は、厚み方向に対して、θｉ１，θｉ２は、上面および下面となる。また、θｉ３は、節点Ｎｉと同じ位置、すなわち、厚肉シェル要素の中立面Ｔとなり、これらの自由度に対して計算される温度から厚み方向の温度分布が求められる。
【００２８】
また、本実施形態では、図４の（ｂ）に示すように、厚み方向の節点自由度で計算された温度分布から厚肉シェル要素の厚み方向の中立面Ｔに対する変位分布が求められる。
更に、本実施形態の特徴として、従来求めることができなかったシェル要素の厚み方向の上面または下面変位を、以下に示す定式化によって求める。
【００２９】
図４の（ｃ）に示すように、厚肉シェル要素の厚みをｗ、節点の並進成分の変位をΔＸ、節点の回転成分の変位をΔφとすると、上面または下面の変位ΔＸ’を次式１で表される関数で求めることができる。
ΔＸ’＝ΔＸ±（ｗ／２）・ｔａｎΔφ  ・・・（１）
そして、このようにして求められた上面または下面の変位ΔＸ’を、上面または下面の評価点に外挿する。
【００３０】
上記定式化によりシェル要素の厚み方向の上面、下面変位分布を得るためには、汎用有限要素法プログラム、例えば、上記の「ＭＳＣ．Ｍａｒｃ」（登録商標）のサブルーチンに設定すればよい。
なお、式１のｗ／２を中立面Ｔからの厚み方向距離にすることで、厚み方向任意位置での変位分布を求めることができ、解析時間の短縮化を損なうことなく、精度を向上することができる。
【００３１】
次に、中央処理装置１１８は、要素１に熱負荷が加えられた初期の場合における計算結果、すなわち、各節点の温度および変位を記憶装置１１９に格納する。
次に、ステップＳ５において、中央処理装置１１８は、要素２に熱負荷が加えられるｔ時間後の場合について、解析を行うため、熱負荷位置を更新する。
次に、ステップＳ６において、中央処理装置１１８は、熱負荷が最後の要素１２０に加えられる時間の場合であるか否かの終了を判定し、終了していない場合には、Ｓ３に戻り、終了している場合には、Ｓ７を実行する。
【００３２】
Ｓ３に戻った場合には、中央処理装置１１８は、記憶装置１１９に格納されている汎用有限要素法プログラムを呼び出し、要素１の初期の場合と同様に、要素２に熱負荷が加えられたｔ時間後の場合における温度および変形を解析する。
次に、中央処理装置１１８は、要素２に熱負荷が加えられた場合における計算結果、すなわち、各節点の温度および変位を記憶装置１１９に格納する。
【００３３】
このように、中央処理装置１１８は、要素２に熱負荷が加えられるｔ時間後の場合、要素３に熱負荷が加えられる２ｔ時間後の場合、要素３に熱負荷が加えられる３ｔ時間後の場合、・・・、と順次に最後の要素１２０まで、節点の温度、変位の更新（Ｓ３）、求解（Ｓ４）および結果格納を繰り返し、解析を行う。
ここで、節点の温度、変位の更新において、或る要素に熱負荷が加えられた場合における解析計算を行う場合では、当該要素の１つ前に熱負荷が加えられた場合における解析計算結果を用いる。
【００３４】
一方、ステップＳ６の判定結果が終了の場合には、ステップＳ７において、中央処理装置１１８は、格納されている変位から、各時間の変形を計算し、ステップＳ８において、最終の位置ズレ量を出力装置１１３に表示し、ステップＳ９でプログラムを終了する。また、ユーザの要求に応じて、各時間の変形を表示することもできる。
【００３５】
図５は、解析結果である温度分布を示す図である。
この解析結果は、マトリックスの最終の１２０番目に熱負荷をかけた直後の温度分布である。
図５（ａ）は、物体１１Ａを２次元のプレーン要素として要素分割した場合の解析結果を、図５（ｂ）は、本実施形態のシェル要素で要素分割した場合の解析結果を、図５（ｃ）は、物体１１Ａを３次元のソリッド要素として要素分割した場合の解析結果を示している。
【００３６】
この図から、プレーン要素の温度分布に比較して、シェル要素の温度分布が、ソリッド要素の温度分布に非常に近いことがわかる。
また、解析に要した計算時間は、プレーン要素の１に対し、厚肉シェル要素が５、ソリッド要素が７９であり、解析時間の短縮化を損なうことなく、精度が大きく向上しているのがわかる。
【００３７】
図６は、解析結果である位置ずれ量評価を示す図である。
すなわち、図７に示すように、ステッパーにおいてウエハ上に、例えば、電子ビームを照射して熱負荷を与えることで所定のパターンが焼き付けられるが、その際に、本来の焼き付け位置が、前の焼き付けによるウエハの熱膨張やウェハの膨張によるウエハとチャック間の滑りに起因してずれ、これを位置ずれとよんでいる。
【００３８】
この位置ずれ量を、今まさに焼き付けている位置でプロットしたのが図６である。
図６（ａ）は、図２（ａ）のＸ方向の位置ずれ量を、図６（ｂ）は、図２（ａ）のＹ方向の位置ずれ量を示している。
また、曲線ａは、プレーン要素の場合を、曲線ｂは、本実施形態の厚肉シェル要素の場合を、曲線ｃは、ソリッド要素の場合の位置ずれ量を示している。
【００３９】
この図から、プレーン要素の位置ずれ量の履歴に比較して、厚肉シェル要素の位置ずれ量の履歴が、ソリッド要素の位置ずれ量の履歴に非常に近いことがわかる。
上述した実施形態では、物体１１Ａを、厚み方向に要素分割することなく厚肉シェル要素で要素分割するとともに、熱負荷による温度分布を求め、温度分布からシェル要素の中立面変位分布を求め、物体の厚み方向任意位置での変位分布を、中立面変位とシェル要素の厚み方向の位置との関数（式１）で近似し、物体に発生する変位分布を求めるようにしたので、解析時間の短縮化を損なうことなく、温度分布の解析精度を大きく向上することができる。
【００４０】
なお、上述した実施形態では、半導体ウェハからなる物体１１Ａに本発明を適用した例について説明したが、本発明はかかる実施形態に限定されるものではなく、例えば、半導体ウェハへの焼き付けパターンの原形であるレチクル等の解析に広く適用することができる。
そして、上述した解析結果を用いて露光装置を設計および製造することにより、物体１１Ａの変位分布を考慮した位置決め精度の高い露光装置を提供することができる。
【００４１】
（第２の実施形態）
図８は、本発明の第２の実施形態の解析対象を示している。
この実施形態の解析対象は、円盤状の半導体ウェハからなる物体１１Ｂとされている。
そして、回路パターンを転写する際には、露光エネルギが要素１から要素ｋ＋ｊまで順次照射され、この露光エネルギによって物体１１Ｂに熱負荷が加えられる。
【００４２】
なお、この実施形態のブロック図の構成は、第１の実施形態のブロック図の構成と同様であるため、図１のブロック図を用いて重複する説明を省略する。
図９は、本実施形態の解析装置１００の動作を示すフローチャートである。
先ず、ユーザは、ステップＳ１において、解析装置１００を入力装置１１２から起動する。
【００４３】
次に、ステップＳ２において、中央処理装置１１８は、解析対象の形状・寸法などの初期データを入力するようにユーザに出力装置１１３を介して指示する。
ユーザは、入力装置１１２から、解析対象の形状・寸法のデータ、解析対象の物性のデータ（熱伝導率、比熱、密度、熱膨張率、ヤング率、ポアソン比）、熱負荷のデータ（熱エネルギー、位置、位置の移動軌跡、移動時間、負荷面積）、解析対象の初期温度、解析対象の外界の温度、熱負荷の位置に対応する要素からいくつ離れた要素までの面積を変更するかを示す値（以下、「面積変更値」とする。）、および、要素の面積の基準値などを入力する。
【００４４】
この実施形態では、要素を正方形とする。この場合には、要素の面積の基準値は、正方形の１辺の値（以下、「辺基準値」とする。）とすることができる。熱負荷の位置の移動は、図８に示すように、各要素に番号を割り当て、この番号によって指定することができる。
面積変更値は、要素距離という量を、辺基準値の要素においてある要素が別のある要素を１番目として何番目の要素であるかを示す値と定義し、この要素距離によって指定することができる。例えば、図１０（ａ）で要素３は、要素２から要素距離が２であり、要素５は、要素２から要素距離が４である。
【００４５】
なお、予めこれらのデータを記憶装置１１９またはＦＤ１０１に格納しておき、解析装置１００は、起動に従って出力装置１１３に格納内容を表示して、ユーザは、これから選択して各データを入力するようにしてもよい。
次に、ステップＳ３において、中央処理装置１１８は、入力された解析対象の形状・寸法のデータを基に、解析対象をモデル化する。
【００４６】
すなわち、図８に示すように、中央処理装置１１８は、物体１１Ｂを解析対象部分について辺基準値を持つ正方形の要素に分割する。そして、各要素に節点を割り当て、座標と初期温度とを割り振る。
本実施形態では、物体１１Ｂを、厚肉シェル要素で分割する。
そして、図４の（ａ）に示したように、加えられた熱負荷Ｑによる厚肉シェル要素の厚み方向の温度は、節点Ｎｉの持つ３つの自由度（θｉ１，θｉ２，θｉ３）として求められる。自由度を持つ位置は、厚み方向に対して、θｉ１，θｉ２は、上面および下面となる。また、θｉ３は、節点Ｎｉと同じ位置、すなわち、厚肉シェル要素の中立面Ｔとなり、これらの自由度に対して計算される温度から厚み方向の温度分布が求められる。
【００４７】
また、本実施形態では、図４の（ｂ）に示したように、厚み方向の節点自由度で計算された温度分布から厚肉シェル要素の厚み方向の中立面Ｔに対する変位分布が求められる。
次に、面積変更値と分割された要素との関係を説明する。例えば、面積変更値が３である場合には、面積変更の対象となる要素は、熱負荷の位置に対応する要素（例えば、図１０（ｃ）の要素３）を１番目とするから、熱負荷の移動方向上にある３番目までの要素（図４（ｃ）の要素１ないし要素５）である。
【００４８】
次に、熱負荷が加えられる物体１１Ｂについて、更に、要素分割を行う。
次に、面積変更値変数と要素距離変数を用意し、面積変更値変数に面積変更値を代入し、要素距離変数に１を代入する。
次に、面積変更値変数で２をべき乗した値で辺基準値を割る。この割った辺基準値を１辺に持つ正方形の要素で熱負荷の位置に対応する要素から要素距離変数にある要素を更に分割する。
【００４９】
次に、面積変更値変数から１を引き、要素距離変数に１を足し、面積変更値変数および要素距離変数を更新し、再び、面積変更値変数で２をべき乗した値で辺基準値を割る。この割った辺基準値を１辺に持つ正方形の要素で熱負荷の位置に対応する要素から要素距離変数にある要素を更に分割する。この操作を面積変更値変数が１となるまで行う。
【００５０】
このように処理することにより、例えば、図８に示したように、熱負荷が、物体１１Ｂの１番目の要素１から順にｋ＋ｊ番目の要素までつづら折りに移動する場合には、要素を更に分割すると初期では図１０（ａ）に示すようになる。
すなわち、面積変更値が３である場合には、熱負荷の位置に対応する要素１が６４個の要素に更に分割され、要素１から要素距離が２である要素２は、１６個の要素に更に分割され、要素１から要素距離が３である要素３は、４個の要素に更に分割される。
【００５１】
次に、このように更に分割された各要素に対して、節点を割り当て、座標および初期温度を割り振る。
このように要素を更により小さい面積を持つ要素に分割することを本明細書では、要素のミクロ化と呼ぶことにする。これに対し、要素をより大きい面積を持つ要素になるように複数の要素を併合することを要素のマクロ化と呼ぶことにする。
【００５２】
なお、上述した説明では、最大の面積を持つ要素の大きさを基に、熱負荷が加えられる要素および面積変更値内の要素を更に分割するようにしたが、熱負荷が加えられる要素の面積の大きさを基に、残余の要素を決定するようにしてもよい。
次に、ステップＳ４において、中央処理装置１１８は、記憶装置１１９に格納されている汎用有限要素法プログラム、例えば、エムエスシーソフトウェア株式会社の「ＭＳＣ．Ｍａｒｃ」（登録商標）を呼び出し、この解析用モデルを汎用有限要素法プログラムに適用して、要素１に熱負荷が加えられた初期の場合における温度分布および変形を解析する。
【００５３】
この解析では、図４の（ａ）に示したように、厚み方向の温度分布が求められる。
また、図４の（ｂ）に示したように、厚み方向の温度分布から、厚肉シェル要素の厚み方向の中立面Ｔに対する変位分布が求められる。
次に、中央処理装置１１８は、要素１に熱負荷が加えられた初期の場合における計算結果、すなわち、各節点の温度および変位を記憶装置１１９に格納する。
【００５４】
次に、ステップＳ５において、中央処理装置１１８は、要素２に熱負荷が加えられるｔ時間後の場合について、解析を行うため、熱負荷位置を更新する。
次に、ステップＳ６において、要素分割の変更が行われる。この要素分割の変更は、熱負荷の移動方向上にある要素に対して要素のミクロ化を行い、熱負荷の移動方向と反対方向上にある要素（既に熱負荷が加えられた要素）に対して、マクロ化を行うことにより行われる。
【００５５】
ここで、要素のミクロ化について説明する。
まず、要素のミクロ化は、要素距離変数に１を代入する。
次に、熱負荷の位置に対応する要素から要素距離変数にある要素において、１辺の長さが１／２になるようにこの要素を更に分割する。
次に、要素距離変数に１を足し、要素距離変数を更新し、再び、熱負荷の位置に対応する要素から要素距離変数にある要素において、１辺の長さが１／２になるようにこの要素を更に分割する。この操作を要素距離変数が面積変更値になるまで行う。
【００５６】
このように熱負荷の移動方向上にある要素に対して要素のミクロ化が行われる。例えば、上述の例では面積変更値が３であるから、図１０（ｂ）に示すように、熱負荷の位置にある要素２が６４個の要素に更に分割され、要素２から要素距離が２である要素３は、１６個の要素に更に分割され、要素２から要素距離が３である要素４は、４個の要素に更に分割される。
【００５７】
次に、このようにミクロ化された各要素に対して、節点を割り当てる。そして、節点に座標を割り振り、初期の解析結果である温度および変位も割り振る。
この場合において、要素のミクロ化によって、新たに節点が生じるが、この新節点の温度および変位は、当該新節点に最も近い２個の節点における初期の温度および変位の値を適当な補間方法で分配する。補間方法には、直線補間、面積補間、距離の逆数補間、３角形補間、３次曲線補間および多項式補間などがある。例えば、当該新節点に最も近い２個の節点における温度および変位の差分を当該新節点から各節点までの距離で重み付けを行って配分する。
【００５８】
図１０（ｅ）に示すように、節点ａ１ａ２ａ３ａ４からなる要素を節点ａ１ｂ１ｂ５ｂ４からなる要素、節点ｂ１ａ２ｂ２ｂ５からなる要素、節点ｂ２ａ３ｂ３ｂ５からなる要素および節点ｂ３ａ４ｂ４ｂ５からなる要素に要素のミクロ化が行われた場合に、新節点ｂ１の温度は、「節点ａ１の温度＋長さａ１ｂ１／長さａ１ａ２×節点ａ１および節点ａ２における温度の差分」となる。
【００５９】
また、新節点ｂ１の変位は、「節点ａ１の変位＋長さａ１ｂ１／長さａ１ａ２×節点ａ１および節点ａ２における変位の差分」となる。他の新節点ｂ２、ｂ３、ｂ４、ｂ５も同様に求められる。なお、新節点ｂ５については、２節点ａ１、ａ３あるいは２節点ａ２、ａ４のいずれか一方によって求める。
次に、要素のマクロ化について説明する。
【００６０】
まず、要素のマクロ化は、要素距離変数に２を代入する。
次に、熱負荷の位置に対応する要素から要素距離変数にある要素において、互いに隣接する４個の要素を併合することによって、新たな要素を作成する。
次に、要素距離変数に１を足し、要素距離変数を更新し、再び、熱負荷の位置に対応する要素から要素距離変数にある要素において、互いに隣接する４個の要素を併合することによって、新たな要素を作成する。この操作を要素距離変数が面積変更値になるまで行う。
【００６１】
このように熱負荷の移動方向と反対方向上にある要素に対して要素のマクロ化が行われる。例えば、上述の例では面積変更値が３であるから、図１０（ｂ）に示すように、要素２から要素距離が２である要素１は、１６個の要素になる。要素２から要素距離が３である要素はないので、要素距離が３である場合についてはマクロ化は行われない。
【００６２】
次に、このようにマクロ化された各要素の節点、座標、温度および変位は、対応するマクロ化前の節点の各値を用いる。
次に、ステップＳ７において、中央処理装置１１８は、熱負荷が最後の要素Ｋ＋ｊに加えられる時間の場合であるか否かの終了を判定し、終了していない場合には、Ｓ３に戻り、終了している場合には、Ｓ８を実行する。
【００６３】
Ｓ３に戻った場合には、中央処理装置１１８は、記憶装置１１９に格納されている汎用有限要素法プログラムを呼び出し、要素１の初期の場合と同様に、要素２に熱負荷が加えられたｔ時間の場合における温度および変形を解析する。
【００６４】
次に、中央処理装置１１８は、要素２に熱負荷が加えられた場合における計算結果、すなわち、各節点の温度および変位を記憶装置１１９に格納する。
このように、中央処理装置１１８は、要素２に熱負荷が加えられるｔ時間後の場合、要素３に熱負荷が加えられる２ｔ時間後の場合、要素３に熱負荷が加えられる３ｔ時間後の場合、・・・、と順次に最後の要素ｋ＋ｊまで、要素分割変更（Ｓ６）、節点の温度、変位の更新（Ｓ３）、求解（Ｓ４）および結果格納を繰り返し、解析を行う。
【００６５】
ここで、節点の温度、変位の更新において、或る要素に熱負荷が加えられた場合における解析計算を行う場合では、当該要素の１つ前に熱負荷が加えられた場合における解析計算結果を用いる。
なお、図１０（ｃ）に２ｔ時間後の場合の要素分割を示し、図１０（ｄ）に３ｔ時間後の場合の要素分割を示す。
【００６６】
一方、ステップＳ７の判定結果が終了の場合には、ステップＳ８において、中央処理装置１１８は、格納されている変位から、各時間の変形を計算し、ステップＳ９において、最終の位置ずれ量を出力装置１１３に表示し、ステップＳ１０でプログラムを終了する。また、ユーザの要求に応じて、各時間の変形を表示することもできる。
【００６７】
上述した実施形態では、物体１１Ｂを、厚肉シェル要素で分割するとともに、熱負荷による温度分布を求め、温度分布からシェル要素の中立面変位分布を求め、物体の厚み方向任意位置での変位分布を、中立面変位とシェル要素の厚み方向の位置との関数（式１）で近似し、物体に発生する変位分布を求めるようにしたので、プレーン要素１とする解析時間に対し、厚肉シェル要素が５、ソリッド要素が１１４と、解析時間の短縮化を損なうことなく、温度分布の解析精度を大きく向上することができる。
【００６８】
また、熱負荷の少なくとも１つの位置について、シェル要素の面積を、シェル要素の位置と熱負荷の位置との距離の関数として変化するようにしたので、全面を等しい面積の要素とする場合に較べて、節点数を少なくすることができる。そして、負荷の移動にともない要素の面積が変更されるので、各時間において、指定された負荷位置について負荷位置の周辺の温度分布および変形が詳細に解析される。よって、負荷の位置の周辺を詳細に解析したいという要求を満たすことができる。
【００６９】
（第３の実施形態）
図１１は、本発明の第３の実施形態の解析対象を示している。
図１１（ａ）に示すように、この実施形態の解析対象は、略円盤状の第１物体１１および第２物体１２が接触面１３を介して互いの盤面を対向させるように所定の力によって接している。
【００７０】
半導体ウェハに回路パターンを露光装置で転写する場合には、図１１（ｂ）に示すように、第１物体１１は、ウェハ２１に対応し、第２物体は、チャック２２に対応する。
図１１（ｂ）は、露光装置のウェハ２１を固定する部分を示す断面図である。図１１（ｂ）において、ウェハ２１は、チャック２２の面上に載せられ、静電引力や真空吸着などにより、ホルダー２４に固定される。ホルダー２４は、ステージ２５に支持・固定される。これらウェハ２１、チャック２２、ホルダー２４およびステージ２５は、必要に応じて露光装置の減圧された真空チャンバ内に納められている。
【００７１】
そして、回路パターンを転写する際には、露光エネルギがウェハ２１に順次に照射され、この露光エネルギによってウェハ２１に熱負荷が加えられる。
なお、この実施形態のブロック図の構成は、第１の実施形態のブロック図の構成と同様であるため、図１のブロック図を用いて重複する説明を省略する。
【００７２】
図１２は、この実施形態の解析装置１００の動作を示すフローチャートである。
先ず、ユーザは、ステップＳ１において、解析装置１００を入力装置１１２から起動する。
中央処理装置１１８は、解析対象の形状・寸法などの初期データを入力するようにユーザに出力装置１１３を介して指示する。
【００７３】
ユーザは、ステップＳ２において、入力装置１１２から、解析対象の形状・寸法のデータ、解析対象の物性のデータ（熱伝導率、比熱、密度、熱膨張率、ヤング率、ポアソン比）、物体間圧、摩擦係数、熱負荷のデータ（熱エネルギー、位置、位置の移動軌跡、移動時間、負荷面積）、解析対象の初期温度、解析対象の外界の温度、熱負荷の位置に対応する要素からいくつ離れた要素までの面積を変更するかを示す値（以下、「面積変更値」とする。）、および、要素の面積の基準値などを入力する。
【００７４】
この実施形態では、要素を正方形とする。この場合には、要素の面積の基準値は、正方形の１辺の値（以下、「辺基準値」とする。）とすることができる。熱負荷の位置の移動は、後述するように、第１物体１１を辺基準値をもつ正方形の要素で分割するので、図１３（ａ）に示すように、各要素に番号を割り当て、この番号によって指定することができる。
【００７５】
ステップＳ３において、中央処理装置１１８は、入力された解析対象の形状・寸法のデータを基に、解析対象をモデル化する。
すなわち、図１３（ｂ）に示すように、中央処理装置１１８は、第１物体１１および第２物体１２を解析対象部分について辺基準値を持つ正方形の要素に分割する。そして、各要素に節点を割り当て、座標と初期温度とを割り振る。
【００７６】
この実施形態では、第１物体１１を、厚肉シェル要素で分割する。また、第２物体１２を二次元の要素に要素分割する。
そして、図４の（ａ）に示したように、加えられた熱負荷Ｑによる厚肉シェル要素の厚み方向の温度は、節点Ｎｉの持つ３つの自由度（θｉ１，θｉ２，θｉ３）として求められる。自由度を持つ位置は、厚み方向に対して、θｉ１，θｉ２は、上面および下面となる。また、θｉ３は、節点Ｎｉと同じ位置、すなわち、厚肉シェル要素の中立面Ｔとなり、これらの自由度に対して計算される温度から厚み方向の温度分布が求められる。
【００７７】
また、この実施形態では、図４の（ｂ）に示したように、厚み方向の節点自由度で計算された温度分布から厚肉シェル要素の厚み方向の中立面Ｔに対する変位分布が求められる。
次に、面積変更値と分割された要素との関係を説明する。例えば、面積変更値が３である場合には、面積変更の対象となる要素は、熱負荷の位置に対応する要素（例えば、図１０（ｃ）の要素３）を１番目とするから、熱負荷の移動方向上にある３番目までの要素（図１０（ｃ）の要素１ないし要素５）である。
【００７８】
次に、熱負荷が加えられる第１物体１１について、更に、要素分割を行う。
次に、面積変更値変数と要素距離変数を用意し、面積変更値変数に面積変更値を代入し、要素距離変数に１を代入する。
次に、面積変更値変数で２をべき乗した値で辺基準値を割る。この割った辺基準値を１辺に持つ正方形の要素で熱負荷の位置に対応する要素から要素距離変数にある要素を更に分割する。
【００７９】
次に、面積変更値変数から１を引き、要素距離変数に１を足し、面積変更値変数および要素距離変数を更新し、再び、面積変更値変数で２をべき乗した値で辺基準値を割る。この割った辺基準値を１辺に持つ正方形の要素で熱負荷の位置に対応する要素から要素距離変数にある要素を更に分割する。この操作を面積変更値変数が１となるまで行う。
【００８０】
このように処理することにより、例えば、図１３（ａ）に示したように、熱負荷が、第１物体１１の１段目の一方の端にある要素１から順に他方の端にある要素ｇまで移動し、２段目の端の要素ｇに近い要素ｇ＋１から順に折り返してもう一方の端にある要素ｈまで移動し、３段目の端の要素ｈに近い要素ｈ＋１に移動し、以下、同様に、最終段の端にある要素ｋ＋ｊまで順につづら折りに移動する場合には、要素を更に分割すると初期では図１０（ａ）に示したようになる。
【００８１】
すなわち、面積変更値が３である場合には、熱負荷の位置に対応する要素１が６４個の要素に更に分割され、要素１から要素距離が２である要素２は、１６個の要素に更に分割され、要素１から要素距離が３である要素３は、４個の要素に更に分割される。
次に、このように更に分割された各要素に対して、節点を割り当て、座標および初期温度を割り振る。
【００８２】
次に、第１物体１１に属する節点とこれに対応する第２物体１２に属する節点との間に２節点要素を定義する。この２節点要素は、解析において２節点間の剛性方程式および熱伝導方程式で構成される支配方程式を決定するものである。以下、これを接触・摩擦要素と呼ぶ。
この接触・摩擦要素の支配方程式のうちの剛性方程式は、相対変位増分と摩擦力との関係とし、図１４（ａ）に示すような折れ線グラフで表現する。
【００８３】
図１４（ａ）において、縦軸は、摩擦力Ｆであり、横軸は、２節点間の相対変位増分ｄＵｖである。摩擦力は、−∞＜ｄＵｖ≦−Ｃｒｉでは一定値の−Ｆｃｏｎｔであり、−Ｃｒｉ≦ｄＵｖ≦Ｃｒｉでは原点を通り傾きＦｃｏｎｔ／Ｃｒｉの正比例関数であり、Ｃｒｉ≦ｄＵｖ＜＋∞では一定値のＦｃｏｎｔである。
Ｃｒｉは、すべりの判定値であり、必要とする変形の精度に依存する。この実施形態では、Ｃｒｉは、熱負荷による節点変位量の１／１００ないし１／１０００に設定する。精度を高める観点から好ましくは１／１０００に設定する。
【００８４】
そして、Ｆｃｏｎｔは、次式２によって求める。
Ｆｃｏｎｔ＝動摩擦係数×物体間圧×代表面積×接触状態関数  ・・・  （２）
ここで、動摩擦係数は、第１物体１１と第２物体１２との間の動摩擦係数である。物体間圧は、第１物体１１と第２物体とを接触させている物体間力であり、例えば、ウェハ２１をチャック２２で固定する場合には、静電引力などによって生じるチャック圧力に相当する。
【００８５】
代表面積は、全体面積に対し１節点当たりが担う面積であり、より具体的に説明すれば、図１０（ｅ）で示すように、節点ｂ５の代表面積は、当該節点ｂ５に隣接する各節点ａ１、ｂ１、ａ２、ｂ２、ａ３、ｂ３、ａ４、ｂ４との１／２の距離にある点ｃ１、ｃ２、ｃ３、ｃ４からなる四辺形ｃ１ｃ２ｃ３ｃ４の面積である。この実施形態では、要素が正方形なので四辺形となるが、要素の形状が異なれば、四辺形とはならない。
【００８６】
そして、接触状態関数は、図１４の（ｂ）に示すように、
▲１▼第１物体１１と第２物体１２とが非接触の範囲では、第１接触状態関数＝０、
▲２▼第１物体１１と第２物体１２との接触力が集中する範囲では、第１接触状態関数＝１＋Ｓａ  ／Ｓｂ  、ここで、Ｓａ  は非接触の面積、Ｓｂ  は接触力集中範囲の面積、
▲３▼上記▲１▼および▲２▼を除く範囲では、第１接触状態関数＝１、・・・  （３）
とする。
【００８７】
また、支配方程式のうちの熱伝導方程式は、
ｑ＝第２接触状態関数×（Ｔ１−Ｔ２）×代表面積÷接触熱抵抗［第２接触状態
関数］                                                    ・・・  （４）
ここで、Ｔ１  は、接触・摩擦要素の第１物体１１に属する節点の温度であり、Ｔ２  は、接触・摩擦要素の第２物体１２に属する節点の温度である。
【００８８】
そして、第２接触状態関数は、第１接触状態関数により、
第１接触状態関数＝０である場合には、第２接触状態関数＝０
第１接触状態関数＞０である場合には、第２接触状態関数＝１・・・  （５）
とする。
【００８９】
また、第２接触状態関数に依存する接触熱抵抗は、
接触熱抵抗［０］＝第１物体１１と第２物体１２との間の媒体の熱抵抗
接触熱抵抗［１］＝第１物体１１と第２物体１２との間の接触熱抵抗・・・  （６）
とする。
【００９０】
第１接触状態関数は、図１４（ｂ）に示すように、第１物体１１に熱負荷が加わることによって、第１物体１１の一定範囲が熱膨張によって第２物体１２から浮き上がり（非接触となり）、この非接触となった範囲の重さｆ１  が第２物体１２と接触するようになった範囲に分力ｆ１１、ｆ１２されると考えることによる。
ｆ１  ＝−（ｆ１１＋ｆ１２）、ここでは、ｆはベクトル            ・・・  （７）
そして、第２接触状態関数は、第１物体１１と第２物体１２とに対応する節点間の熱伝導を接触状態により変化させることを意味する。すなわち、熱伝導の状態は、２物体間が接触状態から浮き上がり状態に変化することに対し、接触熱抵抗から媒体の熱抵抗に変化させることを意味する。
【００９１】
ここで、非接触の範囲およびその面積Ｓａ  並びに接触力集中範囲およびその面積Ｓｂ  は、第１物体１１および第２物体１２を接触面１３の法線を含む或る断面で、汎用有限要素法プログラムによって接触解析を行って求める。なお、汎用有限要素法プログラムは、例えば、エムエスシーソフトウェア株式会社の「ＭＳＣ．Ｍａｒｃ」（登録商標）を使用すればよい。
【００９２】
接触解析の結果、式７の関係および非接触範囲の長さＬａ  が求まるので、面積Ｓａ  は長さＬａ  を直径として求められ、接触力集中範囲の長さＬｂ  が求まるので、面積Ｓｂ  は半径Ｌａ  ／２＋Ｌｂ  の円の面積から面積Ｓａ  を引くことによって求められる。
【００９３】
このような接触・摩擦要素の支配方程式、すなわち、剛性方程式および熱伝導方程式を求める接触・摩擦サブプログラムを記憶装置１１９に格納し、中央処理装置１１８は、ステップＳ２１において、この接触・摩擦サブプログラムを呼び出し、熱負荷の位置に従い各節点間について第１および第２接触状態関数を求め、ステップＳ２２において、２節点間に剛性方程式および熱伝導方程式として割り振る。
【００９４】
この場合において、この支配方程式が非線形となるため、求解は、ニュートン・ラプソン法（Ｎｅｗｔｏｎ−Ｒａｐｈｓｏｎ  ｍｅｔｈｏｄ  ）などの手法が用いられる。そして、この接触・摩擦サブプログラムは、全体剛性マトリックスを作る際には要素剛性（図１４（ａ）のグラフの傾き）を求解の処理（ステップＳ４）に渡し、一方、変位より求まる内力計算の際には要素内力として摩擦力を求解の処理（ステップＳ４）に渡す。
【００９５】
こうして中央処理装置１１８によって、初期における要素、節点、２節点間の支配方程式が決定される。
次に、ステップＳ４において、中央処理装置１１８は、記憶装置１１９に格納されている汎用有限要素法プログラム、例えば、上記の「ＭＳＣ．Ｍａｒｃ」（登録商標）を呼び出し、この解析用モデルを汎用有限要素法プログラムに適用して、要素１に熱負荷が加えられた初期の場合における温度分布および変形を解析する。
【００９６】
この解析では、図４の（ａ）に示したように、加えられた熱負荷Ｑによる厚肉シェル要素の厚み方向の温度は、節点Ｎｉの持つ３つの自由度（θｉ１，θｉ２，θｉ３）として求められる。自由度を持つ位置は、厚み方向に対して、θｉ１，θｉ２は、上面および下面となる。また、θｉ３は、節点Ｎｉと同じ位置、すなわち、厚肉シェル要素の中立面Ｔとなり、これらの自由度に対して計算される温度から厚み方向の温度分布が求められる。
【００９７】
また、この実施形態では、図４の（ｂ）に示したように、厚み方向の節点自由度で計算された温度分布から厚肉シェル要素の厚み方向の中立面Ｔに対する変位分布が求められる。
そして、接触・摩擦サブプログラムを汎用有限要素法プログラムのサブルーチンに設定することにより、接触判定、接触力及び摩擦力の計算に接触面法線方向の変位量を使用する必要がない。このため、接触面法線方向の変位量の求解を行う必要がないので、本発明では、計算時間の短縮化を損なうことなく、少ない計算機資源でプログラムを実行できる。
【００９８】
次に、中央処理装置１１８は、要素１に熱負荷が加えられた初期の場合における計算結果、すなわち、各節点の温度および変位を記憶装置１１９に格納する。
次に、ステップＳ５において、中央処理装置１１８は、要素２に熱負荷が加えられるｔ時間後の場合について、解析を行うため、熱負荷位置を更新する。
次に、ステップＳ６において、要素分割の変更が行われる。この要素分割の変更は、熱負荷の移動方向上にある要素に対して要素のミクロ化を行い、熱負荷の移動方向と反対方向上にある要素（既に熱負荷が加えられた要素）に対して、マクロ化を行うことにより行われる。
【００９９】
なお、要素のミクロ化およびマクロ化は、第１の実施形態と同様であるので詳細な説明を省略する。
次に、ステップＳ７において、中央処理装置１１８は、熱負荷が最後の要素Ｋ＋ｊに加えられる時間の場合であるか否かの終了を判定し、終了していない場合には、ステップＳ３に戻り、終了している場合には、ステップＳ８を実行する。
【０１００】
ステップＳ３に戻った場合には、中央処理装置１１８は、記憶装置１１９に格納されている汎用有限要素法プログラムおよび接触・摩擦サブプログラムを呼び出し、ｔ時間後における要素、節点、および、２節点間の支配方程式を汎用有限要素法プログラムに適用して、要素１の初期の場合と同様に、要素２に熱負荷が加えられたｔ時間後の場合における温度および変形を解析する。
【０１０１】
次に、中央処理装置１１８は、要素２に熱負荷が加えられた場合における計算結果、すなわち、各節点の温度および変位を記憶装置１１９に格納する。
このように、中央処理装置１１８は、要素２に熱負荷が加えられるｔ時間後の場合、要素３に熱負荷が加えられる２ｔ時間後の場合、要素３に熱負荷が加えられる３ｔ時間後の場合、・・・、と順次に最後の要素ｋ＋ｊまで、要素分割変更（ステップＳ６）、節点の温度、変位および支配方程式の更新（ステップＳ３、ステップＳ２１、ステップＳ２２）、求解（ステップＳ４、ステップＳ２１、ステップＳ２２）および結果格納を繰り返し、接触解析を行う。
【０１０２】
ここで、節点の温度、変位および支配方程式の更新において、或る要素に熱負荷が加えられた場合における解析計算を行う場合には、当該要素の１つ前に熱負荷が加えられた場合における解析計算結果を用いる。
一方、ステップＳ７の判定結果が終了の場合には、ステップＳ８において、中央処理装置１１８は、格納されている変位から、各時間の変形を計算し、ステップＳ９において、最終の位置ずれ量を出力装置１１３に表示し、ステップＳ１０においてプログラムを終了する。
【０１０３】
この実施形態においても第１および第２の実施形態と略同様な効果を得ることができるが、この実施形態では、第１物体１１を厚肉シェル要素として扱い、第２物体１２を２次元のプレーン要素として扱って、２節点間の剛性方程式を所定のグラフで与えるようにしたので、接触判定、接触力及び摩擦力の計算に接触面法線方向の変位量を使用する必要がない。
【０１０４】
なお、上述した説明では、要素１から順次にミクロ化・マクロ化を行ったが、１つまたは複数の任意の負荷位置について、ミクロ化・マクロ化を行うようにしてもよい。この場合において、ミクロ化を行わない負荷位置については、上述のミクロ化の処理を飛ばして接触解析を行い、次の負荷位置に移動するようにすればよい。
【０１０５】
また、上述した実施形態では、面積の変更する範囲は、要素距離が３である場合について説明したが、これに限定されるものではない。要素の面積の基準の大きさにもよるが、要素距離は、解析時間を実用的な範囲にさせる観点から、１ないし５にすることが好ましい。要素距離は、解析結果を実験とよく一致させる観点から、３ないし５、より好ましくは５にすることが好ましい。
【０１０６】
さらに、上述した実施形態では、記録媒体としてＦＤの場合について説明したが、これに限定されるものではない。メモリチップ（例えば、ＲＯＭチップ）、コンパクトディスクのメモリ（例えば、ＣＤ−ＲＯＭ、ＣＤ−Ｒ、ＣＤ−ＲＷ）、光磁気ディスク（ＭＯ）、光ディスク（ＰＤ）、ディジタル・ビデオ・ディスクのメモリ（例えば、ＤＶＤ−ＲＯＭ、ＤＶＤ−ＲＡＭ）など各種記録媒体を用いることができる。
【０１０７】
そして、上述した第１から第３の実施形態では、シェル要素に厚肉シェル要素を用いた例について説明したが、これに限定されるものではなく、単なるシェル要素を用いても良い。
【０１０８】
【発明の効果】
以上述べたように、本発明によれば、解析時間を大きく増大することなく、解析精度を大きく向上することができる。
【図面の簡単な説明】
【図１】本発明の実施形態の解析装置を示すブロック図である。
【図２】本発明の第１の実施形態の解析対象を示す説明図である。
【図３】本発明の第１の実施形態の動作を示す流れ図である。
【図４】本発明の第１の実施形態の解析用モデルを示す説明図である。
【図５】本発明の第１の実施形態による解析結果である温度分布を示す説明図である。
【図６】本発明の第１の実施形態による解析結果である位置ずれ量を示す説明図である。
【図７】位置ずれ量を説明するための説明図である。
【図８】本発明の第２の実施形態の解析対象を示す説明図である。
【図９】本発明の第２の実施形態の動作を示す流れ図である。
【図１０】本発明の実施形態のミクロ−マクロモデルを示す説明図である。
【図１１】本発明の第３の実施形態の解析対象を示す説明図である。
【図１２】本発明の第３の実施形態の動作を示す流れ図である。
【図１３】本発明の第３の実施形態の解析用モデルを示す説明図である。
【図１４】本発明の第３の実施形態の接触・摩擦モデルの説明図である。
【符号の説明】
１１Ａ，１１Ｂ  物体
１１  第１物体
１２  第２物体
１３  接触面
２１  ウェハ
２２  チャック
１００  解析装置
１０１  フレキシブルディスク
１１１  ＦＤ装置
１１８  中央処理装置
１１９  記憶装置
【図１】
2005025588.tif 000003
【図２】
2005025588.tif 000004
【図３】
2005025588.tif 000005
【図４】
2005025588.tif 000006
【図５】
2005025588.tif 000007
【図６】
2005025588.tif 000008
【図７】
2005025588.tif 000009
【図８】
2005025588.tif 000010
【図９】
2005025588.tif 000011
【図１０】
2005025588.tif 000012
【図１１】
2005025588.tif 000013
【図１２】
2005025588.tif 000014
【図１３】
2005025588.tif 000015
【図１４】
2005025588.tif 000016