(57)【特許請求の範囲】
【請求項１】
熱硬化性樹脂の物性値に基づく熱硬化性樹脂の３次元流れの有限要素解析装置であって、３次元解析領域を立体的な有限要素に分割した形状データと前記熱硬化性樹脂の密度、粘度、比熱、熱伝導率、拡散係数および反応速度と、境界条件および初期条件とを入力する手段と、入力した前記形状データと前記熱硬化性樹脂の密度、粘度、比熱、熱伝導率、拡散係数および反応速度と境界条件および初期条件とを記憶する補助記憶装置と、立体要素を用いるコントロールボリューム法に基づいて離散化した運動量の式
【数１】
0003721765.tif 000041
【数２】
0003721765.tif 000042
【数３】
0003721765.tif 000043
（ここで、ｕ、ｖ、ｗは、それぞれｘ、ｙ、ｚ方向の樹脂の流速であり、ｎおよびｎ＋１がそれぞれｎ番目およびｎ＋１番目の時間ステップを表わし、ａiiu、ａiiv、ａiiwおよびａiju、ａijv、ａijwが運動量の式を離散化する過程で定められる流速の係数である。
ｂiu,(n)は運動量の式を離散化する過程で定められるｕに関する係数で、ｎステップの時間項、圧力項、重力項が含まれる。ｂiv,(n) 、ｂiw,(n)も同様にそれぞれｖおよびｗに関する係数である）
と、圧力項を除いてコントロールボリューム法に基づいて半離散化した運動量の式
【数４】
0003721765.tif 000044
【数５】
0003721765.tif 000045
【数６】
0003721765.tif 000046
（ここで、
【数７】
0003721765.tif 000047
）をコントロールボリューム法に基づいて離散化した連続の式
【数８】
0003721765.tif 000048
に代入したのち圧力項を離散化してえられる圧力の式
【数９】
0003721765.tif 000049
（ここで、ｐは圧力）
と、コントロールボリューム法に基づいて離散化した運動量の式から導出する流速補正式
【数１０】
0003721765.tif 000050
【数１１】
0003721765.tif 000051
【数１２】
0003721765.tif 000052
とを用いて、
前記補助記憶装置に記憶された形状データと
前記熱硬化性樹脂の密度、粘度、比熱、熱伝導率、拡散係数および反応速度と
前記境界条件および初期条件とから、コンピュータにより流速および圧力を求める反復計算を行い、前記流速と圧力を同一節点上で解を収束させて前記流速と圧力とを求める３次元流れの有限要素解析手段と、前記同一節点上で求められた流速と圧力を記憶する補助記憶装置とを備え、前記流速と圧力により熱硬化性樹脂の３次元流れを求めることを特徴とする３次元の有限要素解析装置。
【請求項２】
コントロールボリュームに流入流出する流束を計算することによって自由表面の移動を計算する解析手段を含む請求項１記載の３次元流れの有限要素解析装置。
【請求項３】
コントロールボリューム法に基づいて離散化したエネルギーの式
【数１３】
0003721765.tif 000053
と物質の拡散の式
【数１４】
0003721765.tif 000054
を反復計算に組み込み、熱硬化性樹脂の流速、圧力、温度および熱硬化反応率を同一節点上で求め、自由表面の移動を計算する解析手段を含む請求項２記載の３次元流れの有限要素解析装置。
【請求項４】
コンピュータを、請求項１記載の３次元流れの有限要素解析装置として動作させるプログラムを記録した記録媒体。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
金属の鋳造、プラスチック成形などの分野の、流動過程の有限要素解析方法および解析装置に関する。
【０００２】
【従来の技術】
従来、非圧縮性流体の流れの有限要素解析を行うばあい、安定な解をえるためには圧力と流速を求める格子点を食い違わせる必要があったこれは、Ｎａｖｉｅｒ−Ｓｔｏｋｅｓ式などの流れの式において圧力が１次微分、流速が２次微分であることに起因している。このため、有限要素法では流速の補間次数より圧力の補間次数を下げることが行われる。２次元での解析においても通常、たとえば４角形の有限要素を用いるばあい、４つの頂点で流速を求め、４角形の中心で圧力を求める必要がある。このとき、流速の補間次数は１次で圧力は０次となる。２次元解析ではＧａｌｅｒｋｉｎ法を用いて同一格子点で流速と圧力を求める方法(Computer Methods in Applied Mechanics and Engineering、 58, (1986)、 J. G. Rice and R. J. Schnipke、 An Equal-Order Velocity-Pressure Formulation that does not Exhibit Spurious Pressure Modes、 P.135-149)があるが、補間関数で重み付けを行うＧａｌｅｒｋｉｎ法は質量、運動量やエネルギーの保存を明示的に考慮した方法ではないため、安定な計算には細かな有限要素分割が必要であった。
【０００３】
３次元流れの解析法として流動コンダクタンスという量を仮定し、流速の変数を消去して未知変数を圧力のみの１つに減らして計算する方法が特開平８−９９４３１号公報、特開平９−１５０４４３号公報などに開示されているが、多孔質中の浸透流を表す近似的なダルシー流れの式にしか用いることができない。
【０００４】
これまで、３次元流れの有限要素解析で自由表面の移動を計算するばあい、流れに乗って移動する質量０のマーカー粒子を仮定して計算する方法が主流であった。また、薄板流れ（Ｈｅｌｅ−Ｓｈａｗ流れ）のように近似解析では、ＶＯＦ（Volume of Fluid）法のように格子点において流体が占める割合を計算し、この値から自由表面の位置を求めることができる。自由表面の変化が簡単なばあいは、変化した流体領域の形状に合わせて、有限要素分割をやり直したり、有限要素を追加または消去していく方法がある。
【０００５】
【発明が解決しようとする課題】
流速と圧力が求まる格子点が食い違っていると、結果を表示するばあい、このままでは、流速と圧力は異なったデータ系列となり、統一的に処理できないという問題があった。また、Ｇａｌｅｒｋｉｎ法を用いた場合は、前述のように定式化においては質量やエネルギーなどの保存則を考慮していないため、安定な解を求めるためには、有限要素分割を細かくする必要があり、結果的に計算時間の増大や解の精度の低下を招いていた。
【０００６】
また、自由表面の計算でマーカー粒子を用いるばあいは、複雑な流れのばあい、時間が経過すると粒子の分布に偏りが発生する。このため、最初から莫大な数のマーカー粒子を用いるか、適当な時間ごとにマーカー粒子の再配置を行う必要があった。これも、結果的に計算時間の増大を招く。解析領域の再有限要素分割や有限要素の追加、消去する方法は、流れが分岐したり、合流するばあいを計算するには、有限要素の形状がひずんだり、重なったりするのを避けるための多くの条件を考慮しておかねばらない。ＶＯＦ法に関しては薄板流れへの適用のみであり、３次元流れへの適用は見当たらない。
【０００７】
本発明は、このような問題点を解決するもので、３次元流れにおける流速と圧力を同一の格子点上で求め、自由表面の移動を容易に計算することができるようにすることにある。
【０００８】
【課題を解決するための手段】
本発明は、３次元流れ解析の有限要素解析法であって、立体型要素を用い、コントロールボリューム法に基づいて離散化した運動量の式と、圧力項を除いて半離散化した運動量の式を離散化した連続の式に代入したあと圧力項を離散化してえられる圧力の式と、離散化した運動量の式から導出する速度補正式とを用いて流速および圧力を求める反復計算を行い、流速と圧力を同一節点上で求めることを特徴とする解析方法である。
【０００９】
また、本発明により、３次元流れ解析の有限要素解析法に関し、立体要素を用い、コントロールボリューム法に基づいて離散化した運動量の式と、圧力項を除いて半離散化した運動量の式を離散化した連続の式に代入したあと圧力項を離散化してえられる圧力の式と、離散化した運動量の式から導出する速度補正式とを用いて流速および圧力を求める反復計算を行い、流速と圧力を同一節点上で求めることを特徴とする解析装置を提供する。
【００１０】
また、本発明は、コントロールボリュームでは流入流出する流束を見積もることにより保存則を考慮した解析となるため、本質的に計算が安定で、移流項を計算するのと同じ方法で自由表面の移動も容易に計算できる解析方法である。
【００１１】
また、本発明により、コントロールボリュームでは流入流出する流束を見積もることにより保存則を考慮した解析となるため、本質的に計算が安定で、移流項を計算するのと同じ方法で自由表面の移動も容易に計算できる解析装置を提供する。
【００１２】
また、本発明は、温度を計算するエネルギーの式と反応を計算する拡散の式を容易に組み込み、コントロールボリュームでは流入流出する流束を見積もることにより保存則を考慮した解析となるため、本質的に計算が安定で、移流項を計算するのと同じ方法で自由表面の移動も容易に計算できる解析方法である。
【００１３】
また、本発明により、温度を計算するエネルギーの式と反応を計算する拡散の式を容易に組み込み、コントロールボリュームでは流入流出する流束を見積もることにより保存則を考慮した解析となるため、本質的に計算が安定で、移流項を計算するのと同じ方法で自由表面の移動も容易に計算できる解析装置を提供する。
【００１４】
また、本発明により、えられた自由表面の移動状況、流速分布、圧力分布、温度分布および物質の濃度分布の結果をもとに成形条件を決定し、この決定した成形条件で成形品を製造する製造方法を提供する。
【００１５】
また、本発明により、３次元流れの有限要素解析法のプログラムを記録した媒体を提供する。
【００１６】
【発明の実施の形態】
以下に、添付図を参照しつつ、本発明による３次元流れの有限要素解析法とその解析装置の一例である、熱硬化性樹脂成形における樹脂流動過程の解析法および解析装置の好ましい実施の形態の例を詳細に説明する。
【００１７】
熱硬化性樹脂の自由表面が移動する３次元流れの式は以下で表される。
【００１８】
【数１】
0003721765.tif 000002
【００１９】
【数２】
0003721765.tif 000003
【００２０】
【数３】
0003721765.tif 000004
【００２１】
【数４】
0003721765.tif 000005
【００２２】
【数５】
0003721765.tif 000006
【００２３】
【数６】
0003721765.tif 000007
【００２４】
【数７】
0003721765.tif 000008
【００２５】
ここでｘ、ｙ、ｚは３次元の空間座標であり、ｕ、ｖ、ｗはそれぞれ座標軸方向の樹脂の流速、ｐ、Ｔ、ｃはそれぞれ圧力、温度および熱硬化性樹脂の既反応基濃度（反応率）である。ｆは流体占有率で、ｆ＝０の位置には流体が存在しないこと、０＜ｆ＜１の位置は自由表面、ｆ＝１の位置は流体が満ちていることを示す。また、ρは密度、ηは粘度、Ｃvは比熱、ｋは熱伝導率、Ｄは拡散係数である。また、ｇx、ｇy、ｇzはそれぞれ重力加速度のｘ方向成分、ｙ方向成分、ｚ方向成分であり、Ｈは反応熱、Ｒは反応速度である。
【００２６】
解析領域を有限要素で分割し、有限要素の各節点に対して、それぞれが重ならないようにコントロールボリュームを定義する。式（１）を、コントロールボリュームで積分し離散化する。積分記号の添え字ＣＶはコントロールボリュームで積分することを表す。左辺第１項の時間微分は以下のように後退差分で表す。
【００２７】
【数８】
0003721765.tif 000009
【００２８】
ここで、Ｎｉは節点ｉに対する補間関数、ｕiは節点ｉでのｘ方向の流速である。また、（ｎ）はｎ回目の時間ステップであることを示し、ｎ＋１ステップとｎステップの時間間隔がΔｔである。
【００２９】
式（１）の左辺第２、３、４項の移流項は非線型であり、単純に通常の補間関数を用いて積分しても解の空間的振動などの不具合が発生する。このため、本発明の一つである、流れの方向を考慮し、コントロールボリューム法にマッチした計算法を適用する。なお、この計算法は式（５）のＴ、式（６）のｃおよび式（７）の自由表面の移動に対しても同様に用いることができる。ｕ、ｖ、ｗ、Ｔ、ｃ、ｆをまとめてφで表せば、移流項のコントロールボリュームに対する積分は、ガウスの定理を適用して、面積分に変換し、近似的にφを積分の外に出してφとすると
【００３０】
【数９】
0003721765.tif 000010
【００３１】
である。ここで、φはコントロールボリューム界面Ｓでの流束ｑsの正負に応じて、上流側の節点の値を取る。また、コントロールボリュームｉに対する流束ｑsは
【００３２】
【数１０】
0003721765.tif 000011
【００３３】
と離散化して計算することができる。
【００３４】
右辺第１、２、３項の拡散項のコントロールボリュームに対する積分はガウスの定理により面積分に変換して補間関数で積分して離散化すれば、以下のとおりである。
【００３５】
【数１１】
0003721765.tif 000012
【００３６】
右辺第４項の圧力項は
【００３７】
【数１２】
0003721765.tif 000013
【００３８】
で計算する。右辺第５項の重力項の積分は以下となる。
【００３９】
【数１３】
0003721765.tif 000014
【００４０】
以上の式（８）〜式（１３）から、ｕiに関して整理すると
【００４１】
【数１４】
0003721765.tif 000015
【００４２】
と表すことができる。ここでｂiu,(n)はｕに関する係数で、ｎステップの時間項、圧力項、重力項が含まれる。同様にｖ、ｗに関して
【００４３】
【数１５】
0003721765.tif 000016
【００４４】
【数１６】
0003721765.tif 000017
【００４５】
である。
【００４６】
またＴ、ｃ、ｆに関しても同様に
【００４７】
【数１７】
0003721765.tif 000018
【００４８】
【数１８】
0003721765.tif 000019
【００４９】
【数１９】
0003721765.tif 000020
【００５０】
ここで、上付き添え字のＴはＴに関する係数、上付き添え字のｃはｃに関する係数、上付き添え字のｆはｆに関する係数であることを示す。
【００５１】
式（４）の連続の式を用いて本発明の一つである圧力の式を導く。式（１）の圧力項以外の補間関数を用いて半離散化すると
【００５２】
【数２０】
0003721765.tif 000021
【００５３】
と表される。ここでｄiu,(n)はｂix,(n)から圧力項に関するものを除いたものとなる。圧力ｐがコントロールボリューム内で一定と仮定し、積分の外に出すと
【００５４】
【数２１】
0003721765.tif 000022
【００５５】
ここで
【００５６】
【数２２】
0003721765.tif 000023
【００５７】
であり
【００５８】
【数２３】
0003721765.tif 000024
【００５９】
である。同様にｖ、ｗに関しても
【００６０】
【数２４】
0003721765.tif 000025
【００６１】
【数２５】
0003721765.tif 000026
【００６２】
である。節点ｉに対するコントロールボリュームで（４）式を積分して離散化すると
【００６３】
【数２６】
0003721765.tif 000027
【００６４】
ここでＮｉは節点ｉに関する補間関数である。これに式（２１）、式（２４）、式（２５）を代入すると次式の離散化した圧力の式をえる。
【００６５】
【数２７】
0003721765.tif 000028
【００６６】
この式をすべてのコントロールボリュームに対して重ね合せ、節点での圧力に関する全体連立方程式をえる。
【００６７】
【数２８】
0003721765.tif 000029
【００６８】
この式は、反復法で簡単に計算でき、求めた圧力の値が空間的振動を示すことはない。
【００６９】
図１は以上の計算式を用いて行う反復計算手順の例を示したフローチャートである。流れの解析では、はじめに３次元解析領域を立体的な有限要素に分割して形状データを入力する。つぎに、流体の物性値を入力し、境界条件を入力し、初期条件を入力し、計算終了条件や反復計算の緩和係数などの解析条件を入力する。解析領域の形状のみに依存する積分を実施しておく。運動量の式（１４）〜（１６）の係数を計算し仮の圧力値で流速を求める。求めた流速で圧力の式（２８）の係数を計算し、圧力を求める。求めた圧力で次式を用いて流速を補正する。
【００７０】
【数２９】
0003721765.tif 000030
【００７１】
【数３０】
0003721765.tif 000031
【００７２】
【数３１】
0003721765.tif 000032
【００７３】
必要であれば、温度の式（１７）、反応の式（１８）、自由表面の移動の式（１９）を解く。解が収束しなければ、流速、圧力を次式で更新し、必要であれば温度、反応率も次式で更新し
【００７４】
【数３２】
0003721765.tif 000033
【００７５】
運動量の式の係数の計算に戻る。式（３２）において、φは流速または圧力または温度または反応率を表し、ｒは緩和係数で０＜ｒ＜２の範囲の値を取り、上付き添え字のＮＥＷは新しく計算された値、ＯＬＤは古い計算値であることを示す。解が収束すれば、流体占有率の式（１９）を用いて流体占有率分布を更新して自由表面を移動させ、つぎの時間のステップの計算に移る。最後の時間ステップであれば、計算を終了し、えられた解析結果をグラフィック処理して、等高線図、ベクトル分布図、グラフなどの形式で表示したプリンタなどに出力する。
【００７６】
図２は前述した３次元流れの有限要素解析に用いる解析装置の例を示す説明図である。補助記憶装置を内蔵するコンピュータ２０１に入力装置２０２、表示装置２０３が接続されている。また、プリンタ２０４、外付けの補助記憶装置２０５が接続されていても良い。入力装置２０２により３次元解析領域形状データと、流体の物性値、境界条件、初期条件を入力し補助記憶装置に保存される。このデータをもとに、コンピュータで解析を行い、えられた解析結果は通常補助記憶装置に保存され、これをたとえばグラフィック処理して表示装置３により表示する。結果の表示はプリンタ装置で行っても良い。
【００７７】
実施の形態１
図３は本発明の一実施の形態にかかわる解析対象である立方キャビティを示す斜視説明図である。図３において、３１は立方キャビティを示しており、ｗはＺ軸方向の流体の流速を示している。
【００７８】
ここでは図３に示すような一辺の長さが１で、上（ｙ＝１）のｘｙ面がｚ方向に速度１で動くばあいの立体キャビディ定常流れに適用した実施例を説明する。
【００７９】
境界条件は上（ｙ＝１）のｘｙ面でｕ＝０、ｖ＝０、ｗ＝１を与え、それ以外の壁面ではｕ＝０、ｖ＝０、ｗ＝０の条件を与えた。初期条件は、境界条件が与えられている節点以外はすべてｕ＝０、ｖ＝０、ｗ＝０、ｐ＝０とした。流体の物性値は、ρ＝１、η＝０．０１とした。また、流速および圧力の緩和係数は両者とも０．３とし、収束の判定は、相対誤差が０．００１以下になることとした。
【００８０】
解析領域を分割するために図４に示した６面体要素を用いる。図４は本発明の実施の形態１にかかわる６面体要素を示す説明図である。図４において、１、２、３、４、５、６、７および８は節点を示しており、４１は６面体要素を示している。また、９は座標軸α、βおよびγの座標軸原点である。コントロールボリュームは、たとえば節点１に対しては図４に示したように、６面体要素の中心と、節点１を共有する３つの面の中心と、節点１を共有する３つの辺の中心で定義する３つの面の内側をこの要素のサブコントロールボリュームとし、節点１を共有する要素に対するすべてのサブコントロールボリュームを加えたものとなる。解析領域は２０×２０×２０＝８０００要素に均等に分割した。なお、この要素における節点ｉに対する補間関数Ｎiとしては一般的な次式を用いた。
【００８１】
【数３３】
0003721765.tif 000034
【００８２】
ここで、α、β、γは局所座標系で０から１までの値を取り、±は節点ｉでのＮiの値が１となるように＋か−となることを示す。全体座標系ｘ、ｙ、ｚで示されている積分は局所座標系にヤコビアンを用いて変換し、ガウスの求積法を用いて数値積分した。
【００８３】
この解析のばあい、定常流れであるため式（１）〜式（３）の時間微分は不要となり、時間ステップは１回のみとなる。また、温度の式（５）、反応の式（６）を計算する必要はない。
【００８４】
解析結果をグラフィック処理して図５の、立方キャビティでの圧力分布を示す説明図に示すような立方キャビティ５１についての圧力分布をえた。流速と圧力が同じ節点上で求まっていないときは圧力を流速が求まっている節点に補間する必要があるばあいが多いが、本発明のように流速と圧力が同じ節点上で求まっているばあいは圧力の補間などを行う必要がなく、グラフィック処理を簡単に行うことができた。この解析条件のばあいレイノルズ数は１００となるが、圧力分布はスムーズであり空間的な振動は見られない。図６は同時に計算された、立方キャビティでの流速ベクトル分布を示す説明図であり、６１は立方キャビティを示している。これにも空間的振動は見られず、良好な計算結果である。
【００８５】
なお、以上の解析は図２に示したハードウエア構成で、コンピュータにエンジニアリングワークステーションを用いて行った。収束までの反復回数は５２回で、計算時間は約３００秒とかなり短い時間で計算を完了した。
【００８６】
実施の形態２
ここでは図７に示すようなインサート物があるキャビティ（１／２モデル）の中に流体が流入する流れに適用した実施の形態を説明する。図７は、本発明にかかわる解析対象のインサート物をもつキャビティを示す説明図である。図７において、７１はキャビティであり、７２はインサート物であり、７３は対称面であり、７４はゲートである。各部の寸法は、Ｌ1が１０ｍｍ、Ｌ2が６ｍｍ、Ｌ3が１ｍｍである。
【００８７】
境界条件は入り口でｕ＝２（ｃｍ／ｓ）、ｖ＝０、ｗ＝０を与え、対称面では対称境界条件、それ以外の壁面ではｕ＝０、ｖ＝０、ｗ＝０の条件を与え、自由表面ではｐ＝０（Ｐａ）を与えた。初期条件は、境界条件が与えられている節点以外はすべてｕ＝０、ｖ＝０、ｗ＝０、ｐ＝０とした。流体の物性値は、ρ＝１（ｇ／ｃｍ3）、η＝１０（Ｐａ・ｓ）とした。また、流速および圧力の緩和係数はそれぞれ０．５と０．３とし、収束の判定は、相対誤差が０．０１以下になることとした。解析領域の分割には実施の形態１と同様な６面体要素を用い、約１００００要素に分割した。この解析のばあい、温度の式（５）、反応の式（６）は計算しない。
【００８８】
自由表面の移動の計算は以下のように移流項を計算するのと同じように式（９）を用いて行った。初期条件として、入り口節点は最初流体で満たされているとして、ｆ＝１とし、それ以外の節点ではｆ＝０とした。つぎに、ｆ＝１の節点に隣接する節点でｆ＜１の節点がｆ＝１になる時間を式（９）を用いて対応するコントロールボリュームに流入する流束を求めることによって行い、これを時間間隔Δｔとしてつぎの時間ステップの計算を行う。隣接するｆ＜０となる節点が多数あるばあいは最小のΔｔを選択する。各時間ステップの計算はｆ＞０の領域に対して行う。すべての節点でｆ＞０となるか、時間ステップが設定した最大数を超えれば計算を完了する。
【００８９】
図８は、解析結果をグラフィック処理して本発明の実施の形態２におけるインサート物をもつキャビティでの流動パターンを時間の順に（ａ）〜（ｄ）に示した説明図である。図８において、８１はキャビティであり、８３は対称面であり、８４はゲートである。キャビティ８１は図７に示したインサート物７２と同様のインサート物を含んでいる（図示せず）。これは、ｆ＝０．５の点が自由表面であるとして、時間とともに自由表面がどのように移動するかを示したものである。また、図９および図１０は、それぞれ本発明の実施の形態２にかかわるインサート物をもつキャビティでの圧力分布、流速分布を示す説明図である。９１はキャビティでの圧力分布、１０１はキャビティでの流速分布をそれぞれ示している。通常の３次元流れ解析に比較して要素分割は荒い方であるが、流動パターン、図９に示した流動途中の圧力分布はスムーズであり空間的な振動は見られない。図１０に示した流動途中の流速においても空間的振動は見られず、良好な計算結果である。これらの結果は、コントロールボリューム法を用いていることによって荒い要素分割でも安定な計算が可能となり、結果として計算時間が短くなることを示している。
【００９０】
なお、以上の解析は図２に示したハードウエア構成で、コンピュータにエンジニアリングワークステーションを用いて行った。非定常で、自由表面の移動も計算しているにもかかわらず計算時間は約１時間とかなり短い時間で計算を完了した。
【００９１】
実施の形態３
ここでは熱硬化性樹脂成形において、実施の形態２と同様に図７に示すキャビティ（１／２モデル）に樹脂が流入する過程の流れ解析の一例として行った実施の形態を説明する。
【００９２】
境界条件は入り口でｕ＝２（ｃｍ／ｓ）、ｖ＝０、ｗ＝０、Ｔ＝１２０℃、ｃ＝１を与え、対称面では対称境界条件、それ以外の壁面ではｕ＝０、ｖ＝０、ｗ＝０、Ｔ＝１８０（℃）の条件を与え、自由表面ではＰ＝０を与えた。初期条件は、境界条件が与えられている節点以外はすべてｕ＝０、ｖ＝０、ｗ＝０、ｐ＝０とした。熱硬化性樹脂にはエポキシ樹脂を用い物性値を測定した。測定した物性値は、ρ＝１．７（ｇ／ｃｍ3）、ηは非ニュートニアンであり、せん断速度γに依存し、温度と反応率にも依存するため以下の式で表され
【００９３】
【数３４】
0003721765.tif 000035
【００９４】
ここで
【００９５】
【数３５】
0003721765.tif 000036
【００９６】
【数３６】
0003721765.tif 000037
【００９７】
ｎ＝０．７７、Ｂ＝４．０９×１０-6（Ｐａ・ｓ）、τ*＝２．５２×１０-9、Ｔb＝１１７９４（Ｋ）、ｃg＝０．５７、δ＝６．９２、またｃv＝０．９７（Ｊ／ｇ）、ｋ＝０．６７（Ｊ／ｍ・ｓ・Ｋ）、Ｈ＝２６０、ｇx＝０、ｇy＝０、ｇz＝９８０（ｃｍ／ｓ2）、Ｄ＝０である。反応速度Ｒには次式を用いた。
【００９８】
【数３７】
0003721765.tif 000038
【００９９】
ここで
【０１００】
【数３８】
0003721765.tif 000039
【０１０１】
【数３９】
0003721765.tif 000040
【０１０２】
であり、Ａ1＝２．４４×１０9（１／ｓ）、Ｅ1＝１２１７６（Ｋ）、Ａ2＝１．４５×１０2（１／ｓ）、Ｅ2＝４１６７（Ｋ）、μ＝０．６７、ν＝０．９５とした。
【０１０３】
また、流速、圧力、温度および反応率の緩和係数はそれぞれ０．５、０．３、０．５および０．５とし、収束の判定は、相対誤差が０．０１以下になることとした。解析領域の分割には実施の形態１と同様な６面体要素を用い、約３０００要素に分割した。自由表面の移動の計算は実施の形態２のばあいと同様に行った。
【０１０４】
図１１は解析結果をグラフィック処理して本発明の実施の形態３におけるインサート物をもつキャビティでの流動パターンを時間の順に（ａ）〜（ｄ）に示した説明図である。図１１において、１１１はキャビティであり、１１３は対称面であり、１１４はゲートである。キャビティ１１１は図７に示したインサート物７２と同様のインサート物を含んでいる（図示せず）。これは、ｆ＝０．５の点が自由表面であるとして、時間とともに自由表面がどのように移動するかを示したものである。また、図１２および図１３は、それぞれ本発明の実施の形態３にかかわるインサート物をもつキャビティでの圧力分布、流速分布を示す説明図である。１２１はキャビティでの圧力分布、１３１はキャビティでの温度分布をそれぞれ示している。さらに、図１４は流動途中の反応率分布を示す説明図である。１４１はキャビティでの反応率分布を示している。図１３および図１４においても、キャビティは、図７に示したインサート物と同様のインサート物を含んでいる（図示せず）。図１２において、Ｐ1は０．９ＭＰａの等圧線を示しており、Ｐ2は０ＭＰａの等圧線を示している。また、図１３において、Ｔ1は１４０℃の等温線を示しており、Ｔ2は１８０℃の等温線を示している。さらに、図１４において、Ｒ1は０．０５の等反応率線を示している。通常の３次元流れ解析に比較して要素分割は荒い方であるが、流動パターン、図１２に示した流動途中の圧力分布を示す等圧線の形状は図１３においても、キャビティは、図７にスムーズであり空間的な振動は見られない。図１３に示した流動途中の温度においても等温線の形状に空間的振動が見られず、良好な計算結果である。図１４に示した流動途中の反応率においても等反応率線の形状に空間的振動は見られず、良好な計算結果である。これらの結果は、コントロールボリューム法を用いていることによって荒い要素分割でも安定な計算が可能となり、結果として計算時間が短くなることを示している。
【０１０５】
なお、以上の解析は図２に示したハードウエア構成で、コンピュータにエンジニアリングワークステーションを用いて行った。非定常で、自由表面の移動と温度と反応率を計算しているにもかかわらず計算時間は約３時間とかなり短い時間で計算を完了した。
【０１０６】
実施の形態４
実施の形態３のばあいと同様だが、樹脂の種類を変えて解析を行った。測定した樹脂の物性値は、粘度に関しては、ｎ＝０．８５、Ｂ＝４．９２×１０-4（Ｐａ・ｓ）、τ*＝２．８０×１０-7、Ｔｂ＝１４７９（Ｋ）、ｃg＝０．４４、δ＝３．３６であり、反応速度式に関しては、Ａ1＝６．７５×１０-1、Ｅ1＝３８２９、Ａ2＝３．８６×１０4（１／ｓ）、Ｅ2＝６４９６（Ｋ）、μ＝０．５２、ν＝０．１０１であり、その他の物性値は実施の形態３のばあいと同じである。
【０１０７】
図１５は本発明の実施の形態４におけるインサート物をもつキャビティでの圧力分布の解析結果を示す説明図である。図１４において、Ｐ3は０．５ＭＰａ等圧線を示しており、Ｐ4は０ＭＰａの等圧線を示している。実施の形態３の圧力分布の解析結果と比較して、インサート物上下の圧力差が小さく、流動圧力も小さい。
【０１０８】
えられたこれらの解析結果は、ここで用いた樹脂の方が成形性に優れ、流動過程においてインサート物を変形させたり移動させたりすることが少ないと考えられる。この樹脂を用いて成形品を製造した結果、インサート物が変形したり成形品表面から見える位置に移動することなどがない良好な成形品をえることができた。
【０１０９】
実施の形態５
本発明の３次元流れの有限要素解析法を図２の装置を用い、エディタープログラムを用いてＦＯＲＴＲＡＮ言語でコーディングして補助記憶装置にソースプログラムとして保存した。もちろんＦＯＲＴＲＡＮ以外の言語、たとえばＰＡＳＣＡＬ言語やＣ言語を用いても良い。保存したソースプログラムをＦＯＲＴＲＡＮ言語のコンパイラープログラムを用いて実行可能なロードモジュールプログラムに変換した。ソースプログラムまたはロードモジュールプログラムもしくはソースプログラムとロードモジュールプログラムをもとにＣＤ−ＲＯＭ成形用のマスター金型を作製し、このマスター金型によってソースプログラムまたはロードモジュールプログラムもしくはソースプログラムとロードモジュールプログラムを保存したＣＤ−ＲＯＭを製造した。同様に、プログラムを保存する媒体としてフロッピーディスクを選び、ソースプログラムまたはロードモジュールプログラムもしくはソースプログラムとロードモジュールプログラムを図２のコンピュータにフレキシブルディスクドライブを接続しフレキシブルディスクに保存した。プログラムを保存する媒体としてはこれらの他にもちろん、ＤＡＴテープ、カセットテープ、ＭＤ、ＤＶＤ−ＲＯＭなどを用いることができる。
【０１１０】
【発明の効果】
本発明にかかわる３次元流れの有限要素解析装置によれば、コントロールボリューム法に基づいた離散化を行い、運動量の式と連続の式から導出した解が空間的振動をおこさない圧力の式を用いることによって計算を安定に行うことができ、比較的荒い要素分割でも精度良く計算できるために、実用的な計算時間で３次元モデルを用いた解析領域形状または成形品形状に忠実な解析を行うことができる。
【０１１１】
また、コントロールボリュームに適した移流項および自由表面の移動の計算方法を用いるため、粘度が小さいばあいや流速が大きいばあいおよび複雑形状のキャビティのなかの自由表面の移動が、実用的な計算時間で精度良い解析結果をえることができる。
【０１１２】
また、圧力と流速、ならびに温度と熱硬化反応率を同じ節点で求めることができるため、圧力を求める節点が流速とずれている一般的な流れ解析法と異なり、結果をグラフィック表示するばあい、圧力を流速が求まる節点まで補間することなく容易に表示などを行うことができる。
【０１１３】
また、成形品形状の設計、成形条件の設定、材料の選定などの成形条件を決定でき、品質に優れ信頼性の高い成形品を効率的に製造することができる。
【０１１４】
また、３次元流れの有限要素解析法をプログラム化してＣＤ−ＲＯＭやフロッピー（登録商標）ディスクなどの媒体に保存すれば、容易に他の装置で利用することができる。
【図面の簡単な説明】
【図１】  本発明の３次元流れの有限要素解析法の手順の例を示すフローチャートである。
【図２】  本発明の３次元流れの有限要素解析装置のハードウエア構成の例を示す説明図である。
【図３】  本発明の一実施の形態における解析対象の立方キャビティを示す説明図である。
【図４】  本発明の一実施の形態における６面体要素を示す説明図である。
【図５】  本発明の一実施の形態における立方キャビティでの圧力分布を示す説明図である。
【図６】  本発明の一実施の形態における立方キャビティでの流速分布を示す説明図である。
【図７】  本発明の一実施の形態における解析対象のインサート物をもつキャビティの１／２モデルを示す説明図である。
【図８】  本発明の一実施の形態におけるインサート物をもつキャビティでの流動パターンを示す説明図である。
【図９】  本発明の一実施の形態におけるインサート物をもつキャビティでの圧力分布を示す説明図である。
【図１０】  本発明の一実施の形態におけるインサート物をもつキャビティでの流速分布を示す説明図である。
【図１１】  本発明の一実施の形態におけるインサート物をもつキャビティでの熱硬化性樹脂の流動パターンを示す説明図である。
【図１２】  本発明の一実施の形態におけるインサート物をもつキャビティでの熱硬化性樹脂の圧力分布を示す説明図である。
【図１３】  本発明の一実施の形態におけるインサート物をもつキャビティでの熱硬化性樹脂の温度分布を示す説明図である。
【図１４】  本発明の一実施の形態におけるインサート物をもつキャビティでの熱硬化性樹脂の反応率分布を示す説明図である。
【図１５】  本発明の一実施の形態におけるインサート物をもつキャビティでの熱硬化性樹脂の圧力分布を示す説明図である。
【符号の説明】
２０１  コンピュータ、２０２  入力装置、２０３  表示装置、２０４  内蔵補助記憶装置、２０５  外付け補助記憶装置。
【図１】
0003721765.tif 000055
【図２】
0003721765.tif 000056
【図３】
0003721765.tif 000057
【図４】
0003721765.tif 000058
【図５】
0003721765.tif 000059
【図６】
0003721765.tif 000060
【図７】
0003721765.tif 000061
【図８】
0003721765.tif 000062
【図９】
0003721765.tif 000063
【図１０】
0003721765.tif 000064
【図１１】
0003721765.tif 000065
【図１２】
0003721765.tif 000066
【図１３】
0003721765.tif 000067
【図１４】
0003721765.tif 000068
【図１５】
0003721765.tif 000069