(57)【要約】【課題】有限要素法を用いる数値解析装置において、Helmholtz分解を適用して精度よく解析する。【解決手段】方程式設定手段１で、解析対象物の物理量を表すベクトル場をHelmholtz分解した変数で解析対象物の変化を規定する微分方程式と制約条件とから最小二乗法により求めた極小条件式を設定する。境界条件設定手段２で、解析対象物の境界条件を設定する。初期条件設定手段３で、解析対象物の初期条件を設定する。数値解法手段４で、極小条件式を準ニュートン法により有限要素法で解く。すなわち、数値計算手段41と補正手段42で計算し、収束するまで繰返し手段43で繰り返す。このようにして、変位要素でHelmholtz分解による解析ができる。渦度レベルで解くことも可能となるので、精度高い解析が可能となる。【選択図】図１
2007286801.tif 000002
【特許請求の範囲】
【請求項１】
解析対象物の境界条件を設定する境界条件設定手段と、解析対象物の初期条件を設定する初期条件設定手段と、解析対象物の物理量を表すベクトル場をHelmholtz分解した変数で解析対象物の変化を規定する微分方程式と制約条件とから最小二乗法により求めた極小条件式を設定する方程式設定手段と、前記極小条件式を準ニュートン法により有限要素法で解く数値解法手段とを具備することを特徴とする数値解析装置。
【請求項２】
前記極小条件式は、不連続な関数が対称的に分布する循環的対称モデルに基づき２倍の等価集中力の項を含むことを特徴とする請求項１記載の数値解析装置。
【請求項３】
前記極小条件式は、１〜1/100の緩和係数を含むことを特徴とする請求項１記載の数値解析装置。
【請求項４】
前記数値解法手段は、前記極小条件式を補正する補正手段を含むことを特徴とする請求項１記載の数値解析装置。
【請求項５】
前記数値解法手段は、ポテンシャルを基底関数とするHelmholtz要素により計算する手段であることを特徴とする請求項１記載の数値解析装置。
【請求項６】
前記数値解法手段は、基本変数を基底関数とする基本変数要素により計算する手段であることを特徴とする請求項１記載の数値解析装置。
【請求項７】
前記数値解法手段は、偏導関数要素と変位要素および圧力要素を用いる偏導関数要素法により計算する手段であって、偏導関数要素値と変位要素値の分布を与えて、解析対象物の変化を規定する微分方程式と連続の式を満足するように圧力要素のノード値をガレルキン法で求める手段と、偏導関数要素値と変位要素値および圧力要素値の分布を与え、ならびに境界条件値に基づいて、解析対象物の変化を規定する微分方程式を満足するように変位要素のノード値をガレルキン法で更新する手段と、偏導関数要素値と変位要素値および圧力要素値の分布を与えて、解析対象物の変化を規定する微分方程式を満足するように偏導関数要素のノード値を更新する手段とにより、収束するまで反復計算する手段を含むことを特徴とする請求項１記載の数値解析装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、数値解析装置に関し、特に、応力を表す高階微分が連続な有限要素法（Ｃ1-連続な要素法またはＣ0-連続な応力を表す高階微分が連続な要素と変位を表す低階微分要素のハイブリッド法）の応力を表す高階微分部分の境界条件処理法としての循環的対称モデル(Cyclic Symmetry Model)を導入して、解析対象物の物理量を表すベクトル場をHelmholtz分解して有限要素法で数値解を求める数値解析装置に関する。
【背景技術】
【０００２】
連続体を離散解析する方法として、有限要素法が確立されている。有限要素法以前からある差分法は、線形要素を適用した有限要素法で説明できる。連続体理論は、無限回連続的微分可能な系を前提として確立されている。無限回連続的微分可能な系とは、基底となる関数が解析対象とする系内で連続であり、その導関数も連続であり、またその導関数も連続というように、無限階の導関数まで連続的に連続する系を言う。これに対して、Ｎ回連続的微分可能な系とは、Ｎ階微分まで連続的に連続する系を言い、Ｃn-連続的微分可能な系などと記す。Ｃ0-連続とは、関数のみ連続で、１階導関数が不連続な場合を言う。
【０００３】
有限要素法は、系を小さな領域に区分（区分領域と言う）し、領域ごとに有限級数を適用する方法である。その領域、すなわち当該関数の適用領域を、有限要素と呼ぶ。なお、領域あるいは区分領域の名称は、後述のように、検査領域や有限体積など他でも使うので、関数の適用範囲を規定する領域を特に有限要素または要素とよぶ。したがって、要素間の連続性が、有限要素で表した系の連続性を決めることになる。その系に適用する要素を、Ｃn-連続な要素（あるいはＣn-要素）と呼ぶ。連続性の異なる複数の要素を適用する系の連続性は、そのうちの最も連続性の低い要素で決定される。
【０００４】
無限回連続的微分可能を前提とする連続体理論に、このようなＣn-連続な有限要素を適用する方法を、弱解法（あるいは離散解法）と言い、その解を弱解と言う。また、その方程式の表し方を、弱形式と言う。これに対して、厳密解や代数解（存在する場合）を強解と呼ぶ。連続体理論は、微分方程式で表し、系のあらゆる位置で満足されることを要求するのに対して、弱形式は、定積分で表した形式を言い、その定積分の位置（領域）で積分値として釣り合っていればよしとする方法である。その定積分の領域には、要素間境界あるいは系の境界を含み、導関数の不連続性（の差）が積分値に加算される。これを等価集中力と呼び、その導入技法が有限要素法の真髄と言える。定積分の結果は、領域に出入りする物理量（質量、運動量、エネルギー等）の差（収支）に等しいので、収支の方に着目して検査する。この場合、図10（ａ）に示すように、積分領域を検査体積と呼び、出入り面（検査体積の全表面）を検査面と呼ぶ。
【０００５】
ゼロの微分はゼロとして微分可能と考えれば、有限要素の内部では無限回連続微分可能である。そこで、有限要素法では、要素内のみ通る検査面を設定して検査する。これには、グリーン・ガウスの定理を適用する方法が定式化されており、例えば２Ｄ（２次元解析）での流体の粘性項は、式(01)で表現される。ただし、Ａは検査領域、Ｃは検査面、Ｓは要素間境界、ｎjは方向ベクトル、νは動粘性係数、∇2はラプラス演算子である。粘性項のモデル化：
−∫Aν∇2ＵidＡ  ⇒  −∫CνΣj{(∂Ｕi/∂ｘj)ｎj}dＣ      （01）
式(01)は、有限要素法の技法のひとつとして知られる有限体積法の公式でもある。⇒はモデル化を表し、単なる数学的な積分では等しくならないが、不連続面の微分の性質を力学的に取り込んだ結果を表す。したがって、⇒の右辺が検査結果である。
【０００６】
式(01)右辺を元の微分方程式に戻るように、単に数学的に変形すれば、式(２)となる。モデル化項の陰的表現と陽的表現式：
−∫CνΣj{(∂Ｕi/∂ｘj)ｎj}dＣ
＝−∫Aν∇2ＵidＡ＋∫SνΣj{(∂Ｕi/∂ｘj)ｎj}dＳ      (02)
式(01)左辺と式(02)右辺を比べれば、モデル化項が抽出される。すなわち、式(02)右辺第２項がモデル化項として加えられていることが判る。従来の有限体積法では、式(01)右辺を用いた。式(01)右辺（式(02)左辺）には、モデル化項が陽には表れていない。そこで、この表式を、モデル化項の陰的表現式と呼び、式(02)右辺の表式を、陽的表現式と呼んで区別する。陰的表現式は、検査面での積分式であり、陽的表現式は、検査体積内の積分式と言える。両者同じである。数学的に同じになるようにした。
【０００７】
線形要素では、式(02)第１項はゼロである。はじめに述べた差分法の差分係数は、式(02)第２項を検査体積（２Ｄでは検査面積、１Ｄでは検査区間）で除したものに等しい。（差分点を未知数として。）陽的表現式のモデル化項からは、要素間境界のせん断力の不連続性を補い、力の伝達のつじつまが合うよう導入した項であることが判る。外力には、要素内せん断力変化率（式(02)第１項）と、このモデル化項の合計で対応する。したがって、このモデル化項を等価集中力と呼ぶ。
【０００８】
系の境界にも、式(01)または式(02)の形で適用する。自然界は無限な広がりで連続しているが、数値計算するために系に境界を設け、そこに力学的メカニズムを設けて、系の内と外を絶縁（連立方程式の上で）させる。このメカニズムを境界条件と言う。境界条件の下に解く問題を境界値問題と言う。境界値問題の境界条件は、仕事量を表す共役な変数のうちのいずれか一方を既知と指定することで、指定を行う。境界条件によって系外とは絶縁されるので、系の外は積分しても、積分しなくても影響しない。あるいはゼロを積分したと考えても同じである。なぜなら、系外の連立方程式は解く必要がないからである。したがって、図10（ｂ）の検査領域を考え、式(02)を式(03)として積分する。ただし、Γは系の基本境界を表す。
境界を含んだモデル化：
−∫Aν∇2ＵidＡ
⇒  −∫Aν∇2ＵidＡ＋∫SνΣj{(∂Ｕi/∂ｘj)ｎj}dＳ
＋∫ΓνΣj{(∂Ｕi/∂ｘj)ｎj}dΓ      （03）
【０００９】
複数の共役な変数の組がある問題では、組の数だけ境界条件がある。それより少ない条件では、境界条件とは言わない。それは単なる系内のメカニズムであり、メカニズムによって連立方程式の変数は絶縁されない。１Ｄのはり（梁）の問題では、図10（ｃ）が無限の大きさの系であり、連立方程式も無限に大きい。図10（ｄ）に示すように、これに境界条件を設けることで計算可能な系を得る。（系とは連立方程式と同義。）はりでは、共役な変数は（たわみ：せん断力）（傾角：モーメント）の２組あり、それぞれ共役な変数の積が仕事量である。境界条件は、固定端で（たわみは既知：せん断力は未知）（傾角は既知：モーメントは未知）、ピン端で（たわみは既知：せん断力は未知）、（傾角は未知：モーメントは既知）となる。簡単のため、はりの剛性を１とすれば、モーメントは２階微分で、せん断力は３階微分に等しい。したがって、これら未知数の共役な変数の関係は、変数を、｛ｗ(0) ｗ(1) ｗ(2) ｗ(3)｝とすれば、共役な変数の組みは、｛ｗ(0) ｗ(3)｝と｛ｗ(1) ｗ(2)｝である。重要なことは、これら共役な変数の積は、いずれも仕事量の単位になることである。
【００１０】
はりの釣合い方程式は４階の微分方程式であり、境界上のこれら変数はその積分常数としてのパラメータである。棒の釣合い方程式は２階の微分方程式であり、共役な変数は変位（ｕ(0)）と力（ｕ(1)）である。簡単のため棒の問題で説明する。棒の共役な変数の組みは、｛u(0) u(1)｝であり、上述のように、棒の釣合い式は２階の微分方程式である。したがって棒の有限体積法には、式(02)が適用できる。
【００１１】
有限要素法には、有限体積法のほかガレルキン法が知られている。有限体積法やガレルキン法およびその他の技法は、総称して重み付き残差法と呼ばれ、共通点が理論的に証明されている。本発明の説明は、ガレルキン法に拠る。釣合い式に変位を乗じれば、単位長さ当たりの仕事量釣合い式となり、これを系全体で積分すれば、系の仕事量釣合い式となる。有限要素法では、式(02)の等価集中力の仕事量も加算する。外力の仕事量を省略して、内部のひずみエネルギーの項のみ記せば、式(03)のＡをΩ（系）に変えて式(04)で表される。２Ｄ、３Ｄ（３次元解析）も説明できるよう、要素の辺積分項の形のまま示した。
棒の有限要素モデル化：
∫Ωｕ∇2ｕdΩ
⇒  ∫Ωｕ∇2ｕdΩ−∫Sｕ(∂ｕ/∂ｘ)dＳ
−∫Γｕ(∂ｕ/∂ｘ)dΓ      （04）
【００１２】
これをガレルキン法では（05）式として連立方程式を組む。
ガレルキン法残差式：
∫Ωδｕ∇2ｕdΩ
⇒  ∫Ωδｕ∇2ｕdΩ−∫Sδｕ(∂ｕ/∂ｘ)dＳ
−∫Γδｕ(∂ｕ/∂ｘ)dΓ      （05）
式(05)は、線形要素であれば、等価集中力の項が表現できる。また、変位は連続するので、適合要素である。要素内で変位が一定の要素は、変位が連続しないので、非適合要素と呼ぶ。適合要素は安定的に計算でき、非適合要素で計算するためには、さらに制約条件が必要である。非適合要素は、本発明の説明には不要である。適合要素が安定的に計算できる理由は、等価集中力の共役変数の一方（変位場）が連続だからである。
【００１３】
線形要素では、ノード未知数は材の両端に設定する。変位が２次式の要素では、材中点にもノードを設定する。３次式では３等分点に２つのノードを設定する。いずれもＣ0-連続である。このように、有限要素法は、変位の連続する適合要素（その微分の連続性は問わない）を対象として確立されている。２次要素の中点ノードや、３次要素の３等分点ノードは、材端の傾角として変数変換して適用しても、数学的には同じである。したがって、Ｃ0-連続であることに変りはない。それら変数が材端で連続する要素はＣ1-連続となり、通常必要ないので適用することはない。しかし、数学的に従来どおりの方法で適用すれば、要素間境界では変形を拘束することとなる。すなわち、応力が２材で等しい条件を与える。したがって、集中外力のあるノードでは不正確となる。
【００１４】
系の自然境界（応力が与えられる境界）では、外力と釣り合わせるよう決定できることは、線形要素に比べて利点となる。しかし、それはＣ1-連続要素の利点ではなく、高次要素の利点であって、Ｃ0-連続な要素も同じである。しかし、基本境界（変位が与えられる境界）では、未知数のまま連立方程式に組み込むこととなる。この場合、境界上の未知数ｕ(1)に対して、連立方程式上何の処置も施していないので、自然境界の場合と同じ処置である。しかし、結果として必要なのは、固定端反力と同じ応力である。ここに矛盾が生じるので、数値計算はできない。
【００１５】
有限要素法による解法では、制約条件としての境界条件と連続の条件（連続式）を満足させることが重要である。境界条件は、自然境界条件（Neumann Condition）と基本境界条件（Dirichlet Condition）に区別される。そのうち、基本境界条件が特に重要である。通常の有限要素法では、基本境界条件で指定される変数（変位、速度）を基底関数とし、変位が要素間で連続し、その導関数は不連続となる要素を用いる。これをゼロ回連続的微分可能な要素（Ｃ0-continuity Element）と呼ぶ。したがって、基本境界条件の方は満足させることができる。連続の式は、速度の発散＝０で表される。例えば、ｘ-ｙ面が境界面となる面上では、ｘ方向およびｙ方向には、勾配はゼロである。したがって、連続の式を満足させるには、ｚ方向勾配もゼロにしたい。連続の式は、３大保存則のうちの質量保存則を表すから重要である。
【００１６】
しかし、Ｃ0-要素では不可能であり、計算結果に従属的に誤差が表れる。勾配を従属的な値としてではなく、独立な変数としてコントロールできる要素は１回連続的微分可能な要素（Ｃ1-continuity Element）と呼ばれる。この要素では、ｚ方向勾配もゼロと置けるので、連続の式は満足される。ただし、ｚ方向の勾配変数は３つの変位に各１つあるから、残りの２つは変数として、自然境界上と同じように取り扱って解析しなければならない。自然境界では、変位を未知数とすれば、応力を表すその微分は、通常ゼロ値（外力がある場合はその値）で与えられる。しかし、基本境界では、応力を表す勾配変数も、またその微分も、ゼロであってはいけない。
【００１７】
したがって、通常の有限要素法では、高次要素であっても、境界上で勾配は未知数としないで済むＣ0-要素（例えば２次要素）が適用される。このように、Ｃ1-要素が適用できる技法はまだ存在しない。Ｃ1-要素が活躍する例としては、板曲げ（曲面でも）が挙げられる。しかし、板曲げは、変位と傾角が基本境界条件として指定されるので、上述のソリッドの理論とは異なる。ただし、板曲げでも、Ｃ2-要素を適用すると、応力まで変数とすることになるので、上述と同じ矛盾が展開される。
【００１８】
応力を未知数とする方法は、差分法で使われる。すなわち、ハイブリッド法である。例えば、板曲げの方程式は、
∇2∇2ｗ＝ｑ
であるから、∇2ｗを未知数としてポアソン方程式を解く。すなわち、応力レベルの微分を未知数とする。このとき、固定端の境界で∇2ｗに対称条件を与えて解く。対称条件とは、境界線を局所原点としたとき、垂直方向には偶関数のみとすることである。したがって、対称線上に鏡を立て、鏡面内の差分点は対応する実差分点と同じ値になるとして解く。このほか、対称条件を与えて解く方法に、フーリエ解析がある。しかし、フーリエ解析は、無限回連続的微分可能な条件に近づくよう解くので、等価集中力モデルとは理論的根拠は異なる。
【００１９】
有限要素法の真髄は、有限級数を適用する技法にあり、このような対象条件に対して、差分法やフーリエ解析のように境界上で連続と仮定するモデルとは異なる。また、基本境界や自然境界でもなく、対象条件でもない境界条件に、サイクリックモデル(周期境界条件モデル)がある。サイクリックモデルは半無限をモデル化する技法であり、離れた位置の２つの境界の値が同じ値（変数）になると仮定して解くモデルである。
【００２０】
Helmholtzの定理について説明する。Helmholtzの定理によれば、任意のベクトル場Ｕは、スカラーポテンシャルφとベクトルポテンシャルΨに分解できる。Helmholtz分解のメリットは次の通りである。Ｕが３成分であるのに対し、φとΨの導関数で、３成分の体積変化あり（非回転成分、Dilatation）と３成分の体積変化なし（回転成分、Rotation）の合計６成分に、物理現象の特性に応じた分解表示ができる。Ψ成分は、連続の式を代数的に満足する。したがって、非圧縮ではφは通常は不要である。また、速度を求めるのに、圧力要素は不要である。圧力は、計算結果の速度を独立変数として与えて、従属変数として反復なしで求める。ただし、Helmhltz分解にはCoulombゲージと呼ばれる制約条件が付いている。Coulombゲージは、ベクトルポテンシャルの発散＝０で表される。上述のように、Helmholtz分解は理論的な利点は大変多いが、この制約条件を満足させることが最大の難関である。その方法が確立されていないことが、Helmholtz分解が一般的に用いられていない原因である。Helmholtz分解のもうひとつの特徴は、ベクトルポテンシャルは２Ｄではベクトルではなくスカラー（流れ関数と呼ばれる）になることである。スカラーには発散＝０の条件は意味がないので、ゲージ条件はない。したがって、２Ｄではメリットのみ得ることができる。
【００２１】
Helmholtz表示による連続体理論は完成されている。Helmholtz分解は、３成分の変位を合計６成分で表すので、力学性状が正確に表現できる。すなわち、体積変化のある変形（圧縮性）とない変形（非圧縮）を表す２種類の変位または速度に分け、３×２の要素で表す。したがって、それぞれの変位に異なる境界条件を与えることができる。例えば、非回転成分を表す変位にはslip条件を与え、回転成分を表す変位にはnon-slip条件を与えるなど、異なる境界条件を別々に与えることができる。これは、Helmholtz分解する以外には方法はない。
【００２２】
有限要素法にはＣ1-要素を適用する方法のほか、ポテンシャル要素とその導関数を要素化するハイブリッド法がある。それぞれＣ0-要素（線形要素など）で済ませる方法である。ただし、この方法ではポテンシャル要素を微分した値は導関数要素の値に完全には一致しない。このため、ポテンシャル要素は導関数要素の弱解（その周りの積分値）として求める。ハイブリッド法には、１次導関数を要素化する方法のほか、２次導関数を要素化する方法がある。差分法で発達した方法であるが、有限要素法にも適用できる。また、この２次導関数要素を適用する方法の方が１次導関数要素を適用する方法より一般的である。２次導関数を用いるハイブリッド法（渦度-流れ関数法）では２回積分値の弱解としてポテンシャル（流れ関数）要素を求める。渦度-流れ関数法は、ゲージ条件の存在しない２Ｄでのみ確立されている。
【００２３】
３Ｄではゲージ条件を無視して計算する方法が試みられているが、ゲージ条件を無視すると一般的化して適用できない。Ｃ1-要素を適用する方法では未知数の数（未知数自由度）は十分ある（８ degrees of freedom / node×３要素＝24未知数）ので、ノード上でのみ条件式（Coulombゲージ式）を満足させる方法がある。すなわち、24ある自由度のうちの１つ（1/24）をCoulombゲージに充てる方法である。このように、特定の点で条件式を満足させる方法は選点法と呼ばれる。選点法の問題点は汎関数の資格がない（極大または極小問題として表せない）ことである。また、Ｃ1-要素は高次要素なので、任意の境界形状へ適合させる特別な技法が必要となるなどの問題点がある。以下に、これらに関連する従来技術の例をいくつかあげる。
【００２４】
特許文献１に開示された「六面体有限要素数値解析方法」は、Ｃ1-連続なHelmholtz要素を用いるボクセル解法において、有限パッチ要素で補綴して任意境界に高精度に近似適合させる方法である。３重３次の基底関数を用いるボクセルで格子を構成する。ベクトル場をHelmholtz分解して、スカラーポテンシャルφとベクトルポテンシャルΨで表す。８個のボクセルからなるパッチ系を構成し、実境界の点（実境界面と基本要素の稜との交点、実境界面と基本要素の表面との相貫線上の１点または複数点、基本要素内の実境界面上の１点または複数点）にパッチ系の中心を重ねて、その点の関数の１階微分値がすべて０になるように補綴する。ボクセル基本系の外周境界からのエネルギー流出がなく、高次ボクセル要素をオイラー型解法に適用できる。このようにＣ1-要素を任意の境界形状へ適合させる技法はあり、また、連続の式を高精度に満足させることができるが、簡単とは言えない。Ｃ0-要素によるハイブリッド法は、境界形状へ適合させることは比較的容易であるが、連続の式の満足度は劣る欠点がある。
【００２５】
非特許文献１には、「高レイノルズ数の非圧縮流体をNavier-Stokes方程式と多重格子法を用いて解く方法」が報告されている。２次元の非圧縮Navier-Stokes方程式の渦度-流れ関数方程式は、高レイノルズ数の細かいメッシュの流れの解析において、強結合多重格子法の有効性の研究に使われる。方形空洞（正方キャビティ）中の強制流は、モデル問題として使われる。レイノルズ数が10000でも、メッシュが257×257点の条件で、解が得られた。流れ場における２次的な渦度の出現があるため、１次元格子集団座標変換の使用のためには、一様なメッシュの細分化が好ましい。この文献は、Coulombゲージのない２Ｄ問題を、流れ関数と渦度を未知数として、差分法で解いたものである。したがって、Helmholtz分解のメリットのみ得られる問題を解いていることになる。しかし、その解は、他のあらゆる研究者が比較資料として参照し、ベンチマーク的存在として認められている。２Ｄの流れ関数-渦度法は有限要素法でも可能である。有限要素法では、Ｃ0-要素を用いるハイブリッド法でも、Ｃ1-要素のみでハイブリッドなしで解く方法でも可能である。
【００２６】
非特許文献２と非特許文献３には、Ｃ1-要素のみ用いて、ハイブリッドなしで解く有限要素法で、非特許文献１と同じ２Ｄの、レイノルズ数が10000の問題を、メッシュが41×41点の条件で解いている。その結果は、非特許文献１の結果をほぼ完全に再現している。その結果を図11（ａ）に示す。この方法では、境界条件と連続の式を完全に満足する。
【００２７】
非特許文献６には、Ｃ0-要素をハイブリッドして解く有限要素法で、非特許文献１と同じ２Ｄのレイノルズ数が10000の問題を、メッシュが121×121点の条件で解いている。その結果は、非特許文献１の結果をよく再現しているが完全ではない。その結果を図11（ｂ）に示す。メッシュが201×201点の条件で解くと著しく改良される。その結果を図11（ｃ）に示す。これは渦度の境界条件値を、流れ関数の数値微分値から求めて、設定していることにも原因がある。したがって、さらに細かなメッシュでは、境界条件値も改善されるので、さらに改良されることになる。
【００２８】
非特許文献４と非特許文献５には、Ｃ1-要素のみ用いて、ハイブリッドなしで解く有限要素法で、Coulombゲージはノード上で、選点法で満足させる方法を用い、プロポーションが２Ｄ問題と同じ１×１の断面で奥行き６とした１×１×６の３Ｄキャビティを、対称条件を入れて奥行き３として、レイノルズ数100の問題を、メッシュが７×７×７点の条件で解いている。このプロポーションの問題は非特許文献７に可視化実験があり、非特許文献１でのレイノルズ数100の計算結果は、実験結果の対称面の状態をよく表していることが示されている。Ｃ1-要素の計算結果を非特許文献１でのレイノルズ数100の計算結果と対称面で比較した結果はよく一致している。これを図11（ｄ）に示す。
【００２９】
【特許文献１】特開2005-202634号公報
【非特許文献１】U. Ghia, K. N. Ghia and C. T. Shin "High-Re Solutions for Incompressible Flow using the Navier-Stokes Equations and a Multigrid Method," Journal of Computational Physics 48, (1982), pp387-411.
【非特許文献２】今村, "Ｃ1-連続なFEMモデルの非圧縮流体への適用" 第11回数値流体力学シンポジウム論文集, pp.567-568, 1997.
【非特許文献３】今村・棚橋, "Helmholtz要素を用いる渦度保存スキームと非回転成分保存スキーム" 第17回風工学シンポジウム論文集, pp.225-230, 2002.
【非特許文献４】J. Imamura, T. Tanahashi "Numerical Scheme for Continuum Mechanics Using C1-Continuity Helmholtz-Element," World Conference Ψ2Ψ3 on Computational Mechanics (CCMVI, Beijing), RS-01-11-6, (2004)
【非特許文献５】J. Imamura, T. Tanahashi, Numerical scheme for incompressible fluid using Helmholtz-elements, 4th Int. Symposium on Adv. Fluid Information, AFI/TFI 2004, Sendai, Japan, MS3-11, pp.115-120, (2004).
【非特許文献６】今村・棚橋, "Ψ-Ｕ-ω法の検討" 第16回数値流体力学シンポジウム講演論文集、講演番号A12-1, 2002.
【非特許文献７】望月・山辺・山田, "二次元キャビティ内の流れ（可視化実験）・・・・・日本機械学会論文集(B) 52巻 476号,(1986), pp1589-1592.
【発明の開示】
【発明が解決しようとする課題】
【００３０】
しかし、従来の数値解析装置（３自由度の変位や速度を用いる方法）では、次のような問題がある。力学の問題は運動方程式を解く問題である。運動量は加速度に密度を乗じて表される。運動方程式は運動量が応力と釣り合ってゼロとなるべきことを示す方程式である。応力は加速度の時間積分で求められる速度や変位の、空間的な変化（問題により１次〜４次の導関数）で表される。したがって、密度を１とすれば、運動方程式は加速度項と空間微分項の釣り合い方程式で表される。通常加速度は時間ピッチΔtの間は一定として、時間ステップをｎ＝０,１,２,３,・と進めながら解析する。ｎ＝０ステップは初期状態であり、すべての値をゼロと置く。
【００３１】
したがって、運動方程式は完全に満足されている。時系列問題は初期値問題と呼び、境界値問題が境界上で仕事量を表す共役な変数（導関数）の半数が独立変数で、残り半数が従属変数であるのに比べ、初期値問題は初期値側で基底関数のすべての導関数が既知であるとする問題である。連続体理論は時間軸にも無限回連続的微分可能を前提とするが、離散解析（有限要素法）は加速度が時空要素の時間境界で不連続と仮定して解析するものである。したがって、運動方程式には誤差（残差）を伴う。その誤差分散を、境界条件の下に、極小化するように解く必要がある。その極小化を表す微分方程式が方程式設定手段である。
【００３２】
従来の数値解析装置（３自由度の変位や速度を用いる方法）では、空間的な自由度の制約（高次の導関数が表せないこと）から方程式設定手段に制約がある。また、物理特性（体積変化のある、なし）に応じた対応ができないなど、多くの制約がある。それらの制約によって生じる問題や課題を列挙すれば、運動方程式の微分レベルである渦度輸送方程式解法に対応できないので非特許文献１の解が再現できない。ロッキング問題（Volume Locking, Shear Locking, Mesh Locking）がある。ロッキング問題への対応策として、圧力（スカラーポテンシャルと同じ働きをする）要素を高次化する必要があるが、そのための離散化inf-sup条件は難解である。ロッキング問題はいすれも、３自由度の変位要素で生じる問題である。Helmholtz分解して要素化すれば、問題自体が存在しなくなる。しかし、Helmholtz分解する連続体理論は確立されているが、数値計算技法はまだ確立されていない。また、自由度不足によって、壁面上で連続の式の誤差をゼロとすることはできない。
【００３３】
通常の３自由度の変位要素では、ポアソン比νが0.5に近い非圧縮材料でVolume Locking問題が生じる。その原因は、(１−２ν)が分母に表れる体積弾性係数は本来直ひずみのみに関与すべきであるにも拘らず、せん断ひずみも同じ変位要素で表すと、その使い分けができなくなることに因る。Shear Lockingも曲げ変位要素とせん断変位要素を分けないことに因るもので、同根である。Mesh Lockingは変位が線形要素で、圧力が一定要素で生じ、圧力を高次化するためにはinf-sup条件を満足するよう高次化しなければならないなどが課題である。Mesh Lockingは過拘束問題とも呼ばれ、例えば系の角（隅）の要素で境界条件と連続の式を満足させるために線形要素の未知数を使いきってしまうと、未知数すべてがゼロとなり、かつその拘束状態が拡大するもので、過拘束問題と言うよりも、自由度不足問題と言える。Helmholtz要素ではこれらの問題は存在しない。
【００３４】
固体の非回転変位は直ひずみを表すので欠かせない。したがって、固体-流体の連成問題で、非圧縮を３要素で表現する流体側は、非圧縮で壁面slip条件とするか、圧縮性を考慮してnon-slip条件とする必要があり、連成条件は２極端のいずれかしか選べない。Helmholtz要素では非回転成分はslip条件として両者の縁を切り、回転成分はnon-slip条件とすることができる。
【００３５】
非特許文献１では渦度輸送方程式を解くので、運動方程式を解くより、よりミクロな空間変化を対象として解き、解精度もよい。しかし、２Ｄしか解析できない。３Ｄ解析するためには、Coulombゲージを満足させる方法、および渦度のソレノイダル条件（渦度の発散＝０）を満足させる方法の開発が必要である。また、流れ関数の数値微分で渦度の境界条件を与えているが、数値微分は桁落ちし易く、有効桁数の確保が難しいので、メッシュを非常に小さくしなければならない。これに対しては、数値積分をベースとする方法の開発が必要である。
【００３６】
本発明の目的は、上記従来の問題を解決するために、非特許文献１と同じHelmholtz分解に基づく方法により、体積変化のある変位要素と体積変化なしの変位要素に分けてロッキング問題を根本的に解決し、同時に壁面境界条件を２種類の変位それぞれにslip条件とnon-slip条件を独立に選べるものとして固体-流体の連成問題での矛盾を解決し、圧力を用いずに速度を求める方法によって離散化inf-sup条件の検討を無用とし、方程式設定手段として渦度輸送方程式を選び、または部分的に渦度輸送方程式を含む方程式設定手段により、ミクロな変動を対象として解析できるものとし、ゲージ条件（Coulombゲージ、速度ソレノイダル条件、渦度ソレノイダル条件）を制約条件として方程式を解く数値解法手段により、非特許文献１で回避した３Ｄ問題への適用を可能とすることである。渦度輸送方程式を高精度に解くためには、境界条件としての渦度を、ベクトルポテンシャル（２Ｄで流れ関数）の数値微分で与える（非特許文献１の方法）のではなく、渦度の境界条件値を未知数として、その数値積分値がベクトルポテンシャルの分布になるよう、数値積分をベースとする数値解法手段が必要である。
【００３７】
このために、有限要素法のモデル化技法を踏まえて、応力レベルの導関数である渦度を未知数として、応力の対称性とサイクリック性を同時に導入する循環的対称モデル（Cyclic Symmetry Model）をＣ1-連続な要素法で確立することも目的となる。さらに、循環的対称モデルを、応力レベルと速度レベルのＣ0-連続な要素によるハイブリッド法の応力レベルの要素に適用して、新たなハイブリッド法を確立することも目的となる。そのハイブリッド法確立のために、微分項ごとに独立した有限要素（Ｃ0-連続な偏導関数要素９個）を充てて連続体を解く方法（偏導関数要素法）を確立する。偏導関数要素法は、応力レベルの導関数が未知数とならないＣ0-要素によるハイブリッド法、すなわち速度要素とベクトルポテンシャル要素のハイブリッド法（Helmholtz要素のハイブリッド法）では、循環的対称モデルなしに計算できるものとなる。
【課題を解決するための手段】
【００３８】
上記の課題を解決するために、本発明では、数値解析装置を、解析対象物の境界条件を設定する境界条件設定手段と、解析対象物の初期条件を設定する初期条件設定手段と、解析対象物の物理量を表すベクトル場をHelmholtz分解した変数で解析対象物の変化を規定する微分方程式と制約条件とから最小二乗法により求めた極小条件式を設定する方程式設定手段と、極小条件式を準ニュートン法により有限要素法で解く数値解法手段とを具備する構成とした。
【発明の効果】
【００３９】
上記のように循環的対称モデルを入れて構成したことにより、Helmholtz分解による解析が、変位要素でできる。一般になじみ易い方法で、しかも渦度レベルで解くことが可能となるので、精度高い解析が可能となる。さらに、Ｃ1-連続な変位要素を、変位レベルとその導関数レベルに分けて、それぞれにＣ0-連続な要素を適用してハイブリッド（複合）し、導関数レベルの要素（偏導関数要素）に循環的対称モデルを適用することで、さらになじみ易い方法となる。
【発明を実施するための最良の形態】
【００４０】
以下、本発明を実施するための最良の形態について、図１〜図９を参照しながら詳細に説明する。
【実施例１】
【００４１】
本発明の実施例１は、解析対象物の物理量を表すベクトル場をHelmholtz分解した変数で表し、解析対象物の変化を規定する微分方程式に非圧縮の運動方程式（Navier-Stokes方程式）を選び、Coulombゲージを制約条件として、最小二乗法による極小条件式を設定し、解析対象物の境界条件と初期条件を設定し、極小条件式を準ニュートン法により有限要素法で解く数値解析装置である。
【００４２】
図１は、本発明の実施例１における数値解析装置の機能ブロック図と動作手順を示す流れ図である。図１において、方程式設定手段１は、解析対象物の物理量を表すベクトル場をHelmholtz分解した変数で解析対象物の変化を規定する微分方程式と制約条件とから最小二乗法により求めた極小条件式を設定する手段である。境界条件設定手段２は、解析対象物の境界条件を設定する手段である。初期条件設定手段３は、解析対象物の初期条件を設定する手段である。数値解法手段４は、極小条件式を準ニュートン法により有限要素法で解く手段である。数値計算手段41は、極小条件式を準ニュートン法で数値計算する手段である。補正手段42は、その計算結果を補正する手段である。繰返し手段43は、数値計算手段41と補正手段42の計算を、収束するまで繰り返す手段である。
【００４３】
図２と図３は、基本原理を説明する図である。図４は、修正なしの計算結果を示す図である。図５は、修正ありの計算結果を示す図である。図６は、計算結果を示す図である。図７は、計算結果を示す図である。
【００４４】
上記のように構成された本発明の実施例１における数値解析装置の機能と動作を説明する。最初に、図１を参照しながら、数値解析装置の機能の概要を説明する。方程式設定手段１で、解析対象物の物理量を表すベクトル場をHelmholtz分解した変数で解析対象物の変化を規定する微分方程式と制約条件とから最小二乗法により求めた極小条件式を設定する。境界条件設定手段２で、解析対象物の境界条件を設定する。初期条件設定手段３で、解析対象物の初期条件を設定する。数値解法手段４で、極小条件式を準ニュートン法により有限要素法で解く。すなわち、数値計算手段41と補正手段42で計算し、収束するまで繰返し手段43で繰り返す。
【００４５】
数値解析装置の動作手順を説明する。ステップ１で、解析対象物の物理量を表すベクトル場をHelmholtz分解した変数で解析対象物の変化を規定する微分方程式と制約条件とから最小二乗法により求めた極小条件式を設定する。ステップ２で、運動方程式の初期条件と境界条件を設定する。ステップ３で、目的方程式を、有限要素法により準ニュートン法により解く。すなわち、数値計算手段41と補正手段42で計算し、収束するまで繰返し手段43で繰り返す。ステップ４で、計算結果を出力する。
【００４６】
解析対象物の変化を規定する微分方程式について説明する。本発明の実施例１は、解析対象物の変化を規定する微分方程式に非圧縮の運動方程式（Navier-Stokes方程式）を選ぶ。非圧縮のNavier-Stokes方程式は密度を１として、加速度と応力変化の釣り合いで表される。すなわち、i方向の速度をＵiとし、応力変化をＦiとすれば、
(∂Ｕi/∂t)＋Ｆi＝0，(ｉ＝ｘ，ｙ，ｚ)
で表される。Ｆiは動粘性係数をνとし、圧力をＰとして
Ｆi≡Σj{∂(ＵiＵj−ν∂Ｕi/∂ｘj＋δijＰ)/∂ｘj}，
(ｊ＝ｘ，ｙ，ｚ)，（ｘx＝ｘ，ｘy＝ｙ，ｘz＝ｚ）
すなわち、
Ｆx≡{∂(ＵxＵx−ν∂Ｕx/∂ｘ＋Ｐ)/∂ｘ}
＋{∂(ＵxＵy−ν∂Ｕx/∂y)/∂ｙ}
＋{∂(ＵxＵz−ν∂Ｕx/∂z)/∂ｚ}
である。
【００４７】
次に、速度ＵiをHelmholtz分解して表す方法を説明する。Helmholtzの定理では、Ｕを任意のベクトル場（Ｕ≡{Ｕx,Ｕy,Ｕz}）とし、φをスカラーポテンシャルとし、Ψをベクトルポテンシャルとすると、
Ｕ＝gradφ＋curlΨ，（divΨ＝０）
となる。非圧縮では、
Ｕ＝curlΨ，（divΨ＝０）
となるのでΨのみで表される。（  ）内はCoulombゲージと呼ばれる制約条件式であり、常に付帯する。
【００４８】
次に、Helmholtz表示されたNavier-Stokes方程式と制約条件式を最小二乗法により求めた極小条件式を設定する方法について説明する。Navier-Stokes方程式の残差ベクトルをReqと記し、制約条件式の残差ベクトルを時間ピッチΔtで除してRccと記す。ReqとRccの単位は同じになるように選んでいるので、最小二乗法が適用できる。ベクトルの二乗は内積をとるとする。
最小二乗方程式：
(1/2)∫Ω[δ(Req)2＋δ(Rcc)2]dΩ＝０      (１)
(Req)＝(∂Ｕ/∂t)＋Ｆ
(Rcc)＝(1/Δt)(divΨ)
Ｕ≡{Ｕx,Ｕy,Ｕz}, Ｆ≡{Ｆx,Ｆy,Ｆz}, Ψ≡{Ψx,Ψy, Ψz}
実際には、式(１)を準ニュートン法的に解く。
【００４９】
次に、数値解法手段で、極小条件式を準ニュートン法で解く方法について説明する。準ニュートン法は定まった方法があるものではなく、ニュートン・ラプソン法に準ずる方法であればよい。ただし、反復計算の初期値が解の近傍に設定されている必要がある。本発明の実施例では、Navier-Stokes方程式の項にはガレルキン法を適用する。ガレルキン法は最小二乗法に準じた方程式をニュートン・ラプソン法で解く方法である。その式は、υを緩和係数（υ＝1〜1/100）として、次の方程式となる。
【００５０】
準ニュートン法方程式：
∫Ω[δＵキ(Req)＋(υ/2)δ(Rcc)2]dΩ
＋∫S[νΣjδＵ?{∂Ｕ/∂ｘj}]dＳ
＋∫Γ[νΣjδＵ?{∂Ｕ/∂ｘj}]dΓ＝０      (２)
Ｓは要素間境界、Γは系の境界である。それらの積分項は要素間等価集中力モデルと系の境界の等価集中力モデルである。Γの積分項はδＵが境界上では変数ではなく、境界値として与えられるので消える。
【００５１】
次に、補正ステップについて説明する。緩和係数υは経験的にυ＝１〜1/100の間で設定する。その設定には、選点法が1/24の自由度をCoulombゲージに使うことが参考になる。しかし、緩和係数を入れても、非圧縮性の制約条件は剛体変形を表すから、数値誤差も系全体に瞬時に伝播する。その誤差を取り除くのが補正ステップである。誤差はRcc（Ｕの非回転成分）がゼロになっていないことであるから、補正値をΔΨとして、最小二乗法で求めてΨを補正する。その際、ΔΨで表したＵの回転成分が最小二乗法の意味でゼロとなる制約条件を付帯させる。その式、および更新する式は、次の方程式となる。
【００５２】
補正ステップの方程式：
(1/2)∫Ω[δ(curlΔΨ/Δt) 2＋υδ(divΔΨ/Δt＋Rcc)2]dΩ＝０
Rcc＝(1/Δt)divΨ
更新ステップ：
Ψ ? Ψ＋ΔΨ
【００５３】
次に、初期値の設定手段について説明する。初期値は系のすべての要素のノード値がゼロとして設定する。これは条件式をすべて満足する解である。以下には、非特許文献４と同じ３×３×６のプロポーションの３Ｄキャビティで説明する。この問題では上蓋が移動壁になっているので、初期値は移動量もゼロとする。時間ステップをΔｔとして、Δｔ後には少しの移動量を与える。例えば1/10程度とし、１になるまで順次増加させる。Δｔは1/100〜1/1000程度とし、レイノルズ数によって経験的に与える。重要なことは、初期値がΔｔ後の解の近傍にあるようにいろいろな値を設定することである。
【００５４】
次に、境界条件の設定手段について説明する。これも、非特許文献４と同じ３×３×６のプロポーションの３Ｄキャビティで説明する。境界条件のテーブルは境界面の向きで種類分けする。すなわち、ボクセルではｉ軸を垂線とする面を(i-1, i+1)面として３種類ある。例えば、ｘ軸を垂線とする面は(ｚ,ｙ)面であり、同様にｙ軸、ｚ軸を垂線とする面の合計３種類である。系の稜線上のノードの境界条件は、稜線に接する２面のテーブルを指定することで与える。系の角点のノードは、同様にして、３面のテーブルを指定することで与える。
【００５５】
２Ｄキャビティをｘ-ｙ面とし、ｚ軸方向に奥行きをつけて３Ｄキャビティを考える。ノード未知数ベクトルを図２（ｄ）により、
{Ψi(000)，Ψi(100)，Ψi(010)，Ψi(110)，
Ψi(001)，Ψi(101)，Ψi(011)，Ψi(111)}
とし、境界値で固定するベクトル要素は１、未知数とする要素は０のバイナリーでテーブル値を与えるとする。そのうち、Ψi(000)および面内２方向の勾配を表す成分は壁面上ですべて１である。したがって、残りの垂直軸方向勾配のみ１または０の区別が必要である。
【００５６】
先ず、移動壁面の速度Ｕ1について考える。移動速度は面内で一定であるが、垂直壁との連続性より、端部でゼロとなる必要がある。かつ、
Ｕ1＝Ψ3(010)−Ψ2(001)
の成分のうち、Ψ2(001)はｚ方向の対称性を考えればゼロでなければならない。よって、移動速度はΨ3(010)値で与えられる。Ｕ2，Ｕ3についてもそれぞれ対応する軸（面）で同様に考えれば、結局、面内の勾配はすべてゼロでなければならない。したがって、Ψ1(000)値もゼロとなる。
【００５７】
ここで、(ｘ,ｚ)面内のゲージ条件について考えてみる。Ψ1(100)＋Ψ2(010)＋Ψ3(001)＝０のうち、ノード点でΨ1(100)およびΨ3(001)はゼロであるから、必然的にノード値Ψ2(010)もゼロでなければならない。かつ、面内変化もゼロでなければならないから、結局ノードベクトルの値はすべて境界条件として与えられることになる。かつ、壁面上でCoulombゲージは完全に満たされる。もし、手前(ｘ,ｙ)壁面を対称面とするなら、ｚ方向勾配に関するノード値はすべてゼロとなるから、それ以外のノード成分に留意してテーブル値を設定すればよい。先ずＵ3＝０およびその(ｘ,ｙ)面内変化を表すノード成分のテーブル値を１とする。
【００５８】
前同様、さらにゲージ条件について考えれば、Ψ3(001)の項は(ｘ,ｙ)面内ですべて値ゼロである。したがって、Ψ1(100)およびΨ2(010)の項についても(ｘ,ｙ)面内で分布ゼロにできればCoulombゲージは完全に満足される。これはΨ1(000)およびΨ2(000)の(ｘ,ｙ)面内での相対変化なしとすることで可能となる。結局、系の１点でΨi(000)値をゼロに固定すれば、すべての壁面でΨi(000)値はゼロとなる。よって、完成した対象面内テーブルは次のようになり、結局２Ｄキャビティと同じになる。
【００５９】
{Ψ1(000)＝１，Ψ1(100)＝１，Ψ1(010)＝１，Ψ1(110)＝１，
Ψ1(001)＝１，Ψ1(101)＝１，Ψ1(011)＝１，Ψ1(111)＝１}
{Ψ2(000)＝１，Ψ2(100)＝１，Ψ2(010)＝１，Ψ2(110)＝１，
Ψ2(001)＝１，Ψ2(101)＝１，Ψ2(011)＝１，Ψ2(111)＝１}
{Ψ3(000)＝０，Ψ3(100)＝０，Ψ3(010)＝０，Ψ3(110)＝０，
Ψ3(001)＝１，Ψ3(101)＝１，Ψ3(011)＝１，Ψ3(111)＝１}
【００６０】
境界値の大きさ自体は、２Ｄキャビティ断面では、移動壁端部でゼロとし、３次式で立ち上げて補間する。３Ｄでも同様とし、それらをｚ方向端部でも３次補間して、いずれの端部でも垂直壁の値（ゼロ値）との連続性を確保する。以上のように、Helmholtz要素では完璧な形で境界条件を与えることができる。
【００６１】
次に、繰り返し手段について説明する。解析対象物の物理量を表すベクトル場をHelmholtz分解した変数で解析対象物の変化を規定する微分方程式と制約条件とから最小二乗法により求めた極小条件式を設定する方程式設定手段１と、解析対象物の境界条件を設定する境界条件設定手段２と、解析対象物の初期条件を設定する初期条件設定手段３とで組立てられた連立方程式を、極小条件式を準ニュートン法により有限要素法で解く数値計算手段41と、その解を補正する補正手段42とで解くことで、第１近似解を得る。次いで、第１近似解を初期値として、第２近似解を計算する。これを繰り返す。第２近似解と第１近似解の差が小さくなれば、収束したと判断して、繰り返し計算を終了し、Δｔ後の速度の状態を、加速度から計算する。
【００６２】
Δｔを小さく設定しているので、Δｔ後の速度の状態は、次のΔｔ後の初期値の資格がある。よって、時間を進めて計算し、所望の時間まで、または所望の速度分布の状態になるまで計算する。計算が完了すれば結果を出力する。
【００６３】
図４は、補正ステップ42がない場合の計算結果のコンター図である。対称となるべき分布が対称でないことや、誤差によって解が振動している。これらから正しい解でないと判定される。図５は、補正ステップを入れた場合の計算結果のコンター図である。この計算では、υ＝1/100に設定した。図６(ａ)は、対称面の速度プロファイルを、非特許文献１の２Ｄ解と比較した図である。レイノルズ数（Re）が100程度ではよく合うが、1000程度になると若干の誤差が生じている。図６(ｂ)と(ｃ)は、２Ｄでレイノルズ数1000を計算した結果である。非特許文献１の結果とよく合っていることが分かる。２Ｄ計算と３Ｄ計算の違いは，Coulombゲージの制約がある/なしの違いである。
【００６４】
図７(ａ)は、３Ｄ計算のz方向のコンターである。これによれば、両端の壁の影響がz方向の1/4点まで及んでおり、定性的に誤差が含まれていることが判る。図７(ｂ)は、Ｃ1-連続なHelmholtz要素を、Ｃ1-連続な速度要素で計算した結果である。非特許文献１の結果とよく合っていることが分かる。したがって、Helmholtz要素を用いる方法でυ＝1/100では、Coulombゲージが十分に保証されていないことが判る。しかし、重要なことは、ゲージ条件が十分保証されていないにもかかわらず、対称性は補正ステップによって確保されることである。
【００６５】
上記のように、本発明の実施例１では、数値解析装置を、解析対象物の物理量を表すベクトル場をHelmholtz分解した変数で解析対象物の変化を規定する微分方程式と緩和係数を乗じた制約条件とから最小二乗法により求めた極小条件式を設定し、解析対象物の境界条件と初期条件を設定し、極小条件式を準ニュートン法によりＣ1-連続な有限要素法で解く構成とし、解いた結果に含まれる誤差を補正する方法とした。緩和係数を乗じることで安定して解けることと、結果に含まれる誤差を補正するところが、本発明の実施例１の特徴である。
【実施例２】
【００６６】
本発明の実施例２は、解析対象物の物理量を表すベクトル場を、速度（基本変数）を変数としてＣ1-連続な要素で表し、解析対象物の変化を規定する微分方程式に非圧縮の渦度輸送方程式（Navier-Stokes方程式の回転）を選び、速度のソレノイダル条件を制約条件として、最小二乗法による極小条件式を設定し、解析対象物の境界条件と初期条件を設定し、極小条件式を準ニュートン法により有限要素法で解く数値解析装置である。
【００６７】
Helmholtz分解の連鎖則を説明する。これは、Helmholtz分解の一般化である。Helmholtzの定理では、Ｕを任意のベクトル場とし、φをスカラーポテンシャルとし、Ψをベクトルポテンシャルとすると、
Ｕ＝gradφ＋curlΨ，（divΨ＝０）
となる。Helmholtzの定理は、任意のベクトル場に適用できるので、項の単位が同じなら、他のベクトル場も表せる。そこで、微分レベルを１段階ずらせば、
Ｕ(1)＝gradλ＋curlβ，（divβ＝０）
となる。Ｕ≡Ｕ(0)と表し、その１階微分のレベルをＵ(1)とした。λは、新たなスカラー変数として導入したものである。
【００６８】
αを非回転速度（irrotational）(dilatation)とし、βを、回転速度(solenoidal) (rotation)とする。
Ｕ(0)＝α＋β，α≡gradφ，β≡curlΨ，（divΨ＝０）
Ｕ(1)＝κ＋ω，κ≡gradλ，ω≡curlβ，（divβ＝０）
いずれもHelmholtz分解のひとつである。（  ）の中は制約条件式として常に付帯する。以上の操作は、連鎖的かつエンドレスに分解可能である。その規則性をHelmholtz分解の連鎖の法則と呼ぶ。したがって、Helmholtz分解の計算装置は微分方程式のアナロジーを利用して、速度レベルで分解表示する計算装置へ移行させることができる。非圧縮流体では次のようになる。
【００６９】
Ｕ(0)＝β，β≡curlΨ，（divΨ＝０）
Ｕ(1)＝ω，ω≡curlβ，（divβ＝０）
divβ＝(∂βx/∂ｘ)＋(∂βy/∂ｙ)＋(∂βz/∂ｚ)
(ωx，ωy，ωz)T≡curl(βx，βy，βz)T
≡{(∂βy/∂ｚ)−(∂βz/∂ｙ)，(∂βz/∂ｘ)−(∂βx/∂ｚ)，
(∂βx/∂ｙ)−(∂βy/∂ｘ)}T
【００７０】
（１）ポテンシャルφ、Ψを基底関数に選ぶ場合は、Helmholtz要素法という。
（２）基本変数α、βを基底関数に選ぶ場合は、基本変数要素法という。
【００７１】
次に、速度レベルもαとβにHelmholtz分解すると流体の運動方程式がうまく解ける理由を説明する。変位ベクトルＵは３成分である。
Ｕ≡{Ｕx, Ｕy, Ｕz}T，（ＵxはＵのｘ成分）（Ｔは転置行列）
変位の偏導関数は、全部で９個ある。これらを、Ｕxyなどで表す。
Ｕxx≡∂Ｕx/∂ｘ，Ｕxy≡∂Ｕx/∂ｙ，Ｕxz≡∂Ｕx/∂ｚ
Ｕyx≡∂Ｕy/∂ｘ，Ｕyy≡∂Ｕy/∂ｙ，Ｕyz≡∂Ｕx/∂ｚ
Ｕzx≡∂Ｕz/∂ｘ，Ｕzy≡∂Ｕz/∂ｙ，Ｕzz≡∂Ｕz/∂ｚ
【００７２】
変位の偏導関数を行列で表したものを、速度勾配行列という。
[Ｕij]≡[∂Ｕi/∂ｘj]
速度勾配行列は、対称行列と反対称行列（交代行列）の和で表すことができる。これを和分解と呼ぶ。
[Ｕij]＝(1/2)[SUij]＋(1/2)[AUij]
対称行列は、その転置行列と同じになる。反対称行列は、転置すると符号が反転する。したがって、反対称行列は、対角成分がゼロとなる。
【００７３】
連続体の力学的挙動は、Ｕxyなどで表すより、和分解して[Ｕij]の対称行列の成分や反対称行列の成分で表す方が、その観察結果を方程式で表し易い。すなわち、[Ｕij]の反対称行列の成分で、渦度を表すことができる。[Ｕij]の対称行列の成分で、速度ひずみ（固体ではせん断ひずみ）を表すことができる。対称行列の成分は、さらに対角成分と非対角成分に分けると、力学的挙動を説明し易い。速度勾配行列対角成分は、直ひずみと呼ばれる。対称行列の非対角成分は、せん断ひずみと呼ばれる。このように２つのグループに分ける理由は、ひずみに比例する応力を表すとき、比例常数（体積弾性係数とせん断弾性係数）が異なるからである。
【００７４】
速度勾配行列の対角成分は、さらに２つのグループに分けると都合がよい。対角成分はＵxxなどであるから、それらの和（Ｕxx＋Ｕyy＋Ｕzz）は、divＵであり、体積変化を表す。すなわち、各成分はゼロでなくても、divＵがゼロとなることが要求されるグループと、体積変化のあるグループに分ける。前者は、後者の一般解と言える。すなわち、後者は特殊解であり、圧縮性の解は両方の和で表される。非圧縮性の解は前者のみで表される。体積力（ボディーフォース）はゼロを仮定し、境界でのみ外力が与えられる場合を考える。
【００７５】
以上をツリー状に表し、ベクトルで示すと次のとおりとなる。
速度勾配行列
速度ひずみテンソル（対称行列）
対角成分
体積変化あり（ε）：ε≡(Ｕxx，Ｕyy，Ｕzz)T，(divＵ≠０)
体積変化なし（κ）：κ≡(Ｕxx，Ｕyy，Ｕzz)T，(divＵ＝０)
非対角成分
せん断歪ベクトル（γ）：γ≡imiＵ
渦度テンソル（反対称行列）
渦度ベクトル（ω）：ω≡curlＵ
【００７６】
ただし、演算子imiとcurlは、次のとおりである。
imiＵ≡{(∂Ｕz/∂ｙ)＋(∂Ｕy/∂ｚ), (∂Ｕx/∂ｚ)＋(∂Ｕz/∂ｘ),
(∂Ｕy/∂ｘ)＋(∂Ｕx/∂ｙ)}T
curlＵ≡{(∂Ｕz/∂ｙ)−(∂Ｕy/∂ｚ), (∂Ｕx/∂ｚ)−(∂Ｕz/∂ｘ),
(∂Ｕy/∂ｘ)−(∂Ｕx/∂ｙ)}T
imiＵ≡{SU23，SU31，SU12}T
curlＵ≡{AU23，AU31，AU12}T
【００７７】
物理現象を和分解行列で表せば、和分解行列の成分が物理現象ごとの卓越する成分を表す。精密に表す場合にはすべての成分が必要であるが、通常は卓越成分で表すだけで十分である。固体では、ばらばらになって飛散する現象以外は渦度を無視できる。流体では特別な場合以外はせん断ひずみを無視できる。非圧縮流ではさらに、対角成分は一般解としてのdivＵ＝０の成分だけ必要で、特殊解成分は必要ない。そこで、divＵ＝０の成分をα要素ではなく、β要素で表す。したがって、divβ＝０がβ要素の制約条件（ゲージ条件、ソレノイダル条件）である。
【００７８】
本発明の実施例２は、Ｃ1-連続な基本変数要素βを非圧縮流体に適用する例で説明する。
【００７９】
解析対象物の変化を規定する微分方程式について説明する。本発明の実施例２は、解析対象物の変化を規定する微分方程式に非圧縮の渦度輸送方程式を選ぶ。渦度輸送方程式は密度を１として、非圧縮のNavier-Stokes方程式の回転で表される。すなわち、ｉ方向の速度をＵiとし、応力変化をＦiとすれば、Navier-Stokes方程式は、
(∂Ｕi/∂t)＋Ｆi＝0，(ｉ＝ｘ，ｙ，ｚ)
で表される。Ｆiは動粘性係数をνとし、圧力をＰとして
Ｆi≡Σj{∂(ＵiＵj−ν∂Ｕi/∂ｘj＋δijＰ)/∂ｘj}，
(ｊ＝ｘ，ｙ，ｚ)，（ｘx＝ｘ，ｘy＝ｙ，ｘz＝ｚ）
【００８０】
すなわち、
Ｆx≡{∂(ＵxＵx−ν∂Ｕx/∂ｘ＋Ｐ)/∂ｘ}
＋{∂(ＵxＵy−ν∂Ｕx/∂y)/∂ｙ}
＋{∂(ＵxＵz−ν∂Ｕx/∂z)/∂ｚ}
であるから、Navier-Stokes方程式をベクトル表示すれば、
Ｕ≡{Ｕx,Ｕy,Ｕz}, Ｆ≡{Ｆx,Ｆy,Ｆz}として、
∂Ｕ/∂t＋Ｆ＝０
であり、渦度輸送方程式はその回転が釣り合うのであるから
curl(∂Ｕ/∂t＋Ｆ)＝０
である。
【００８１】
次に、非圧縮である渦度輸送方程式を基礎変数レベルの変数でHelmholtz表示する。速度は非圧縮では、Ｕ＝β, （divβ＝０）であり、回転をとると圧力は消えるから、
Ｇi≡Σj{∂(βiβj−ν∂βi/∂ｘj)/∂ｘj}，
(ｊ＝ｘ，ｙ，ｚ)，（ｘx＝ｘ，ｘy＝ｙ，ｘz＝ｚ）
Ｇ≡{Ｇx,Ｇy,Ｇz}
として、
curl(∂β/∂t＋Ｇ)＝0, （divβ＝０）
で表される。
【００８２】
次に、Helmholtz表示された渦度輸送方程式と制約条件式を最小二乗法により求めた極小条件式を設定する方法について説明する。渦度輸送方程式の残差ベクトルをReqと記し、制約条件式（ソレノイダル条件）の残差ベクトルを時間ピッチΔtで除してRccと記す。ReqとRccの単位は同じになるように選んでいるので、最小二乗法が適用できる。ベクトルの二乗は内積をとるとする。
最小二乗方程式：
(1/2)∫Ω[δ(Req)2＋δ(Rcc)2]dΩ＝０      (３)
Req＝curl(∂β/∂t＋Ｇ) 
Rcc＝(1/Δt)divβ
実際には、式(３)を準ニュートン法的に解く。
【００８３】
次に、数値解法手段で、極小条件式を準ニュートン法で解く方法について説明する。準ニュートン法は定まった方法があるものではなく、ニュートン・ラプソン法に準ずる方法であればよい。ただし、反復計算の初期値が解の近傍に設定されている必要がある。本発明の実施例では、渦度輸送方程式の項にはガレルキン法を適用する。ガレルキン法は最小二乗法に準じた方程式をニュートン・ラプソン法で解く方法である。その式は、υを緩和係数（υ＝１〜1/100）として、次の方程式となる。
【００８４】
準ニュートン法方程式：
∫Ω[δcurlβ?(Req)＋(υ/2)δ(Rcc)2]dΩ
＋∫S[νΣjδcurlβ?{∂(curlβ)/∂ｘj}]dＳ
＋∫Γ[νΣjδcurlβ?{∂(curlβ)/∂ｘj}]dΓ
＋∫Γ[νΣjδcurlβ?{∂(curlβ)/∂ｘj}]dΓ＝０      (４)
Ｓは要素間境界、Γは系の境界である。それらの積分項は要素間等価集中力モデルと系の境界の等価集中力モデルである。Γの積分項が２つあるのは、応力の対称性を入れるためである。これを循環的対称モデルと呼ぶ。
【００８５】
循環的対称モデルについて説明する。基本境界では対称条件を入れて解く。これを図３（ａ）で説明する。陽的モデル化項は部分積分項の符号を替えて外力として機械的に加えればよいモデルである。対称かつ不連続な境界では、隣の系の等価集中力もあるのでこれを加えるモデルである。これは図３(ｂ)の周期境界条件（離れたａ端とｃ端を連続させる条件）を入れた対象な系を解くことと同じである。基本境界のΓの積分項は、Ｃ1-連続な要素では、δcurlβのうちβ(000)および境界面に沿ったノード未知数は境界上では変数ではなく、境界値として与えられるので消える。残る境界面垂直方向勾配を表すノード未知数のうち、例えば、境界面がx-z平面ならばβxおよびβzの面内勾配は既知であるから、∂βy/∂yは連続の式より決定できる。また、その面内微分（２階微分）のノードも同様にして決定できる。したがって、ｘ-ｚ境界面ではβy要素のノード未知数はすべて決定されることになる。したがって、変数となるのは残り２要素の境界面垂直方向勾配（および境界に沿う２階微分）のノード未知数である。これらノード未知数は循環的対称モデルに組み込まれて、連立方程式で内的に決定されることとなる。
【００８６】
これを簡単な１Ｄ問題のＣ1-要素の場合で説明する。Ｃ1-要素のノード未知数をαとその勾配θとし、curlβの代わりに関数ｆ(α, θ)の場合で考えれば、（δｆ）は（∂ｆ/∂α）と（∂ｆ/∂θ）より成るが、境界でαが与えられる問題では変分（∂ｆ/∂α）はない。したがって、連立方程式の未知数αに対応する行はα＝０の境界条件で置き換える。未知数θは変数のままとすれば、（∂ｆ/∂θ）の行は連立方程式に残るので、内的に求めることになる。２Ｄ、３Ｄの場合でも全く同じである。
【００８７】
次に、補正ステップについて説明する。緩和係数υは経験的にυ＝1で設定する。しかし、非圧縮性の制約条件は剛体変形を表すから、数値誤差も系全体に瞬時に伝播する。その誤差を取り除くのが補正ステップである。誤差はRccがゼロになっていないことであるから、補正値をΔβとして、最小二乗法で求めてβを補正する。その際、Δβで表したβの回転成分が最小二乗法の意味でゼロとなる制約条件を付帯させる。その式、および更新する式は、次の方程式となる。
【００８８】
補正ステップの方程式：
(1/2)∫Ω[δ(curlΔβ/Δt) 2＋υδ(divΔβ/Δt＋Rcc)2]dΩ＝０
Rcc＝(1/Δt)divβ
更新ステップ：
β ?  β＋Δβ
【００８９】
次に、初期値の設定手段について説明する。初期値は系のすべての要素のノード値がゼロとして設定する。これは条件式をすべて満足する解である。以下には、非特許文献４と同じ３×３×６のプロポーションの３Ｄキャビティで説明する。この問題では上蓋が移動壁になっているので、初期値は移動量もゼロとする。時間ステップをΔｔとして、Δｔ後には少しの移動量を与える。例えば1/10程度とし、１になるまで順次増加させる。Δｔは1/100〜1/1000程度とし、レイノルズ数によって経験的に与える。重要なことは、初期値がΔｔ後の解の近傍にあるようにいろいろな値を設定することである。
【００９０】
次に、境界条件の設定手段について説明する。これも、非特許文献４と同じ３×３×６のプロポーションの３Ｄキャビティで説明する。基本変数要素の境界条件設定は簡単であり、基本境界では境界面に沿った勾配および変位が外的に与えられる。したがって、境界面垂直方向の残りの未知数を循環的対称モデルで内的に決定することとなる。ただし、ゲージ条件（速度のソレノイダル条件）に関する導関数は外的に与える。
【００９１】
ここで、(ｘ,ｚ)面内のゲージ条件について考えてみる。β1(100)＋β2(010)＋β3(001)＝０のうち、ノード点でβ1(100)およびβ3(001)はゼロであるから、必然的にノード値β2(010)もゼロでなければならない。かつ、面内変化もゼロでなければならないから、β2(110),β2(011) ,β2(011)もゼロでなければならない。β1(100)およびβ3(001)にゼロ以外の値が与えられる場合にも同様とする。
【００９２】
境界値の大きさ自体は、２Ｄキャビティ断面では、移動壁端部でゼロとし、３次式で立ち上げて補間する。３Ｄでも同様とし、それらをｚ方向端部でも３次補間して、いずれの端部でも垂直壁の値（ゼロ値）との連続性を確保する。
【００９３】
次に、繰り返し手段について説明する。解析対象物の物理量を表すベクトル場をHelmholtz分解した変数で解析対象物の変化を規定する微分方程式と制約条件とから最小二乗法により求めた極小条件式を設定する方程式設定手段１と、解析対象物の境界条件を設定する境界条件設定手段２と、解析対象物の初期条件を設定する初期条件設定手段３とで組立てられた連立方程式を、極小条件式を準ニュートン法により有限要素法で解く数値計算手段41と、その解を補正する補正手段42とで解くことで、第１近似解を得る。次いで、第１近似解を初期値として、第２近似解を計算する。これを繰り返す。第２近似解と第１近似解の差が小さくなれば、収束したと判断して、繰り返し計算を終了し、Δｔ後の速度の状態を、加速度から計算する。
【００９４】
Δｔを小さく設定しているので、Δｔ後の速度の状態は、次のΔｔ後の初期値の資格がある。よって、時間を進めて計算し、所望の時間まで、または所望の速度分布の状態になるまで計算する。計算が完了すれば結果を出力する。
【００９５】
図８は、Ｃ1-連続な速度要素で計算した結果である。非特許文献１の結果とよく合っていることが分かる。基本変数要素の数値計算結果について説明する。(10)式と(11)式により、基本変数要素を用いて計算した結果の速度プロファイルを、図８に示す。Helmholtz要素法で解析対象とした奥行き６倍の３Ｄキャビティを、同じくRe数100と1000のケースについて幾つかの格子数で数値計算した。いずれもυ＝1.0としている。Re数100の要素数（８×８×12：対称条件なし）の速度プロファイルを、図８（ａ)に示す。図８（ｂ）〜（ｄ）は、結果の対称性を見るための速度コンターである。視線方向-（横軸−縦軸）の区別を下段に示した。対称性は充分保たれていることが判る。図８（ａ）の速度プロファイルが２Ｄと多少ずれるのは、ｚ方向の格子数が粗いことも一因である。
【００９６】
Re数1000（要素数12×12×20：対称条件あり）の速度プロファイルを、図８（ｅ）に示す。Ghiaらとのずれは多少あるが、Helmholtz要素の結果と比べると、主渦の速度プロファイルが直線となって、２Ｄの場合と一致している。He1mholtz要素の場合と同じように、その速度のコンターを、各幾何学的対称面で作図した。これを図８（ｆ）に示す。速度分布図より、ｚ方向への壁面の影響は、対称断面にまで大きく及んではいない。Helmholtz要素との大きな違いは、υ＝1.0でも安定的に計算できることにある。
【００９７】
本発明の実施例２では、数値解析装置を、解析対象物の物理量を表すベクトル場を基本変数（変位、速度）レベルでHelmholtz分解して変数表示し、解析対象物の変化を規定する微分方程式に渦度レベルの高階の微分方程式を選び、分解した変数で循環的対称モデルを適用して表した微分方程式と緩和係数を乗じた制約条件とから最小二乗法により求めた極小条件式を設定し、解析対象物の境界条件と初期条件を設定し、極小条件式を準ニュートン法によりＣ1-連続な有限要素法で解く構成とし、解いた結果に含まれる誤差を補正する方法とした。緩和係数はυ＝１でも安定して解けることと、循環的対称モデルによって、適合要素として要求される連続性よりも高い連続性のＣn-連続な有限要素を適用できるように有限要素法を一般化したところが、本発明の実施例２の特徴である。高階の微分方程式を対象として解析できることで、よりミクロな変化を捉えることができるようになったことが、本発明の実施例２の最も大きな特徴である。
【実施例３】
【００９８】
本発明の実施例３は、解析対象物の物理量を表すベクトル場を、速度（基本変数）を変数としてＣ0-連続な要素で表し、速度の偏導関数９つも同じくＣ0-連続な要素で表して、両者をハイブリッドして解く偏導関数要素法と呼ぶ方法による。実施例３では、要素として、速度要素群と偏導関数要素のほか、Ｃ0-連続な圧力要素も適用する。解析対象物の変化を規定する微分方程式に非圧縮の運動方程式（Navier-Stokes方程式）、および運動方程式の勾配の方程式、並びに連続の式（速度のソレノイダル条件）を選ぶ。運動方程式の空間微分項（局所加速度項以外の項）を、そのうちの移流速度には速度要素の値を代入して、速度勾配には偏導関数要素の値を代入し、粘性項にも偏導関数要素の値を代入して、圧力勾配には圧力要素の値を代入して表す。実施例３は、４ステップより構成される。ステップ１では、圧力要素の値を求める。非圧縮流れの圧力要素の値は、圧力ポアソン方程式を解いて得られることが知られている。圧力ポアソン方程式は上述のように表された空間微分項の発散をとることで表される。ステップ２では、速度要素の値を求める。速度要素の加速度値は、上述のように表された空間微分項と釣り合わなくてはならないから、最小二乗法による極小条件式を設定し、解析対象物の境界条件と初期条件を設定し、極小条件式を準ニュートン法により有限要素法で解くことで求める。ステップ３では、ステップ２で求めた速度要素の値から、ソレノイダル誤差を除くための補正を行う。ステップ４では、偏導関数要素の値を求める。各偏導関数要素の加速度は、上述のように表された空間微分項の各勾配にそれぞれ釣り合わなくてはならない。また、各偏導関数要素の加速度は、速度要素の加速度の勾配とも等しくなくてはならないから、２つの残差分散を同時に、最適に最小化するために、最小二乗法による極小条件式を設定し、系の基本境界では循環的対称モデルによってモデル化し、解析対象物の初期条件を設定し、極小条件式を準ニュートン法により有限要素法で解いて求める。ステップ１〜ステップ４は反復計算するので、順序は必ずしもこのとおりでなくてもよい。
【００９９】
本発明の実施例３は、Helmholtz分解の連鎖則を適用して計算する実施例２の方法を、Ｃ0-連続な要素で実施する装置である。Helmholtzの定理では、Ｕ(0)を速度のベクトル場とし、φをスカラーポテンシャルとし、Ψをベクトルポテンシャルとすると、
Ｕ(0)＝gradφ＋curlΨ，（divΨ＝０）
となる。Helmholtz分解の連鎖則を適用して、微分レベルを１段階ずらせば、
Ｕ(1)＝gradλ＋curlβ，（divβ＝０）
となる。Ｕ≡Ｕ(0)と表し、その１階微分のレベルをＵ(1)とした。λは、新たなスカラー変数として導入したものである。
【０１００】
αを非回転速度（irrotational）(dilatation)とし、βを、回転速度(solenoidal) (rotation)とする。
Ｕ(0)＝α＋β，α≡gradφ，β≡curlΨ，（divΨ＝０）
Ｕ(1)＝κ＋ω，κ≡gradλ，ω≡curlβ，（divβ＝０）
いずれもHelmholtz分解のひとつである。（  ）の中は制約条件式として常に付帯する。非圧縮流体では次のようになる。
Ｕ(0)＝β，β≡curlΨ，（divΨ＝０）
Ｕ(1)＝ω，ω≡curlβ，（divβ＝０）
divβ＝(∂βx/∂ｘ)＋(∂βy/∂ｙ)＋(∂βz/∂ｚ)
(ωx，ωy，ωz)T≡curl(βx，βy，βz)T
≡{(∂βy/∂ｚ)−(∂βz/∂ｙ)，(∂βz/∂ｘ)−(∂βx/∂ｚ)，
(∂βx/∂ｙ)−(∂βy/∂ｘ)}T
【０１０１】
（１）ポテンシャルφ、Ψを基底関数に選ぶ場合は、Helmholtz要素法という。
（２）基本変数α、βを基底関数に選ぶ場合は、基本変数要素法という。
（３）微分レベルの異なる複数要素を（１）または（２）に適用する方法を、ハイブリッド法あるいは偏導関数要素法という。
本発明の実施例３は、（２）の基本変数要素とその導関数要素をハイブリッドする方法である。
【０１０２】
変位ベクトルＵは３成分である。
Ｕ≡{Ｕx, Ｕy, Ｕz}T，（ＵxはＵのｘ成分）（Ｔは転置行列）
変位の偏導関数は、全部で９個ある。これらを、Ｕxyなどで表す。
Ｕxx≡∂Ｕx/∂ｘ，Ｕxy≡∂Ｕx/∂ｙ，Ｕxz≡∂Ｕx/∂ｚ
Ｕyx≡∂Ｕy/∂ｘ，Ｕyy≡∂Ｕy/∂ｙ，Ｕyz≡∂Ｕx/∂ｚ
Ｕzx≡∂Ｕz/∂ｘ，Ｕzy≡∂Ｕz/∂ｙ，Ｕzz≡∂Ｕz/∂ｚ
【０１０３】
変位の偏導関数を行列で表したものを、速度勾配行列という。
[Ｕij]≡[∂Ｕi/∂ｘj]
速度勾配行列は、対称行列と反対称行列（交代行列）の和で表すことができる。これを和分解と呼ぶ。
[Ｕij]＝(1/2)[SUij]＋(1/2)[AUij]
対称行列は、その転置行列と同じになる。反対称行列は、転置すると符号が反転する。したがって、反対称行列は、対角成分がゼロとなる。
【０１０４】
連続体の力学的挙動は、Ｕxyなどで表すより、和分解して[Ｕij]の対称行列の成分や反対称行列の成分で表す方が、その観察結果を方程式で表し易い。すなわち、[Ｕij]の反対称行列の成分で、渦度を表すことができる。[Ｕij]の対称行列の成分で、速度ひずみ（固体ではせん断ひずみ）を表すことができる。対称行列の成分は、さらに対角成分と非対角成分に分けると、力学的挙動を説明し易い。速度勾配行列対角成分は、直ひずみと呼ばれる。対称行列の非対角成分は、せん断ひずみと呼ばれる。このように２つのグループに分ける理由は、ひずみに比例する応力を表すとき、比例常数（体積弾性係数とせん断弾性係数）が異なるからである。
【０１０５】
速度勾配行列の対角成分は、さらに２つのグループに分けると都合がよい。対角成分はＵxxなどであるから、それらの和（Ｕxx＋Ｕyy＋Ｕzz）は、divＵであり、体積変化を表す。すなわち、各成分はゼロでなくても、divＵがゼロとなることが要求されるグループと、体積変化のあるグループに分ける。非圧縮性の解は前者のみで表される。体積力（ボディーフォース）はゼロを仮定し、境界でのみ外力が与えられる場合を考える。
【０１０６】
以上をツリー状に表し、ベクトルで示すと次のとおりとなる。
速度勾配行列
速度ひずみテンソル（対称行列）
対角成分
体積変化あり（ε）：ε≡(Ｕxx，Ｕyy，Ｕzz)T，(divＵ≠０)
体積変化なし（κ）：κ≡(Ｕxx，Ｕyy，Ｕzz)T，(divＵ＝０)
非対角成分
せん断歪ベクトル（γ）：γ≡imiＵ
渦度テンソル（反対称行列）
渦度ベクトル（ω）：ω≡curlＵ
【０１０７】
ただし、演算子imiとcurlは、次のとおりである。
imiＵ≡{(∂Ｕz/∂ｙ)＋(∂Ｕy/∂ｚ), (∂Ｕx/∂ｚ)＋(∂Ｕz/∂ｘ),
(∂Ｕy/∂ｘ)＋(∂Ｕx/∂ｙ)}T
curlＵ≡{(∂Ｕz/∂ｙ)−(∂Ｕy/∂ｚ), (∂Ｕx/∂ｚ)−(∂Ｕz/∂ｘ),
(∂Ｕy/∂ｘ)−(∂Ｕx/∂ｙ)}T
imiＵ≡{SU23，SU31，SU12}T
curlＵ≡{AU23，AU31，AU12}T
【０１０８】
物理現象を和分解行列で表せば、行列の成分が物理現象ごとの卓越する成分を表す。精密に表す場合にはすべての成分が必要であるが、通常は卓越成分で表すだけで十分である。固体では、ばらばらになって飛散する現象以外は渦度を無視できる。流体では特別な場合以外はせん断ひずみを無視できる。非圧縮流ではさらに、対角成分は一般解としてのdivＵ＝０の成分だけ必要で、特殊解成分は必要ない。そこで、divＵ＝０の成分をα要素ではなく、β要素で表す。したがって、それがβ要素の制約条件（ゲージ条件、ソレノイダル条件）となっている。
【０１０９】
本発明の実施例３は、Ｃ0-連続な基本変数要素βと、同じくＣ0-連続な９つの偏導関数要素群をハイブリッドする方法であるが、非圧縮流体に適用する例で説明する。Ｃ0-連続な基本変数要素をβ≡{βx, βy, βz}≡{βi}と記号し、偏導関数要素を、
{∂βx/∂x, ∂βx/∂y, ∂βx/∂z}≡{β1,1, β1,2, β1,3},
{∂βy/∂x, ∂βy/∂y, ∂βy/∂z}≡{β2,1, β2,2, β2,3},
{∂βz/∂x, ∂βz/∂y, ∂βz/∂z}≡{β3,1, β3,2, β3,3}
と記号し、一般化して{βi,j}とも記号する。
【０１１０】
解析対象物の変化を規定する微分方程式について説明する。本発明の実施例３は、解析対象物の変化を規定する微分方程式に非圧縮の運動方程式の９つの勾配の方程式と連続の式を選ぶ。Helmholtz分解は、上述のように、問題の特性に応じて卓越項を選ぶのに適している。しかし、元の速度勾配行列の非対角項すべてを選べば、渦度とせん断変形いずれをも選ぶことに等しい。せん断変形も体積変化はないから、同じβ要素で表すことができる。（あるいは、β要素の勾配でせん断変形も表すことができる。）したがって、Helmholtz分解で重要な点は、体積変化のあるひずみは速度勾配行列の対角項にのみ表れることである。そのうち、非圧縮の成分（divＵ＝０の成分）はβ要素で表すとした。かつ、その制約条件がβ要素の制約条件（ゲージ条件、ソレノイダル条件）とした。よって、非圧縮流れはせん断変形も入れて、β要素で表すことができる。実施例３が実施例１や実施例２と異なる点は、せん断変形を入れたことで、圧力要素が必要となる点である。すなわち、i方向の速度をＵiとし、応力変化をＦiとすれば、Navier-Stokes方程式は
(∂Ｕi/∂t)＋Ｆi＝0，(ｉ＝ｘ，ｙ，ｚ)
で表される。空間微分項Ｆiは動粘性係数をνとし、圧力をＰとして、
Ｆi≡Σj{∂(ＵiＵj−ν∂Ｕi/∂ｘj＋δijＰ)/∂ｘj)}，
(ｊ＝ｘ，ｙ，ｚ)，（ｘx＝ｘ，ｘy＝ｙ，ｘz＝ｚ）
【０１１１】
すなわち、
Ｆx≡{∂(ＵxＵx−ν∂Ｕx/∂ｘ＋Ｐ)/∂ｘ}
＋{∂(ＵxＵy−ν∂Ｕx/∂y)/∂ｙ}
＋{∂(ＵxＵz−ν∂Ｕx/∂z)/∂ｚ}
であるから、Navier-Stokes方程式をベクトル表示すれば、
Ｕ≡{Ｕx,Ｕy,Ｕz}, Ｆ≡{Ｆx,Ｆy,Ｆz}として、
∂Ｕ/∂t＋Ｆ＝０
であり、運動方程式の勾配の輸送方程式は、
∂(∂Ｕ/∂t＋Ｆ)/∂xj＝０
であり、合計９つの方程式となる。連続の式は、
divＵ＝0
である。
【０１１２】
非圧縮では、Ｕ＝β, （divβ＝０）であるから、
Ｆi≡Σj{∂(βiβj−ν∂βi/∂ｘj＋Ｐ)/∂ｘj}，
(i＝ｘ，ｙ，ｚ)，(ｊ＝ｘ，ｙ，ｚ)，（ｘx＝ｘ，ｘy＝ｙ，ｘz＝ｚ）
Ｆ≡{Ｆx,Ｆy,Ｆz}
で表される。連続の式は、
divβ＝０
である。
【０１１３】
次に、初期値の設定手段について説明する。初期値は系のすべての要素（{βi}も{βi,j}も）のノード値がゼロとして設定する。これは条件式をすべて満足する解である。以下には、非特許文献４と同じ３×３×６のプロポーションの３Ｄキャビティで説明する。この問題では上蓋が移動壁になっているので、初期値は移動量もゼロとする。時間ステップをΔtとして、Δt後には少しの移動量を与える。例えば1/10程度とし、１になるまで順次増加させる。Δtは1/100〜1/1000程度とし、レイノルズ数によって経験的に与える。重要なことは、初期値がΔt後の解の近傍にあるように、いろいろな値を設定することである。
【０１１４】
次に、境界条件の設定手段について説明する。これも、非特許文献４と同じ３×３×６のプロポーションの３Ｄキャビティで説明する。Ｃ0-連続な基本変数要素{βi}の境界条件設定は簡単であり、基本境界では境界面に沿った{βi}のノード値が外的に与えられる。したがって、偏導関数要素のうち、境界面に沿った勾配を表す偏導関数要素のノード値も外的に決定される。偏導関数要素のうち、境界上で未定なノード値は境界面垂直方向の勾配を表す３つの要素である。例えば、境界面がx-z平面ならばβx,y,βy,y,βz,y要素である。このうち、βy,y要素のノード値は、連続の式βx,x＋βy,y＋βz,z＝０から決定される。したがって、βx,y,βz,y要素のノード値は変数として残る。この２つの値は、境界面の循環的対称モデルで内的に決定することとなる。初期値はいずれもゼロである。
【０１１５】
境界値の大きさ自体は、２Ｄキャビティ断面では、移動壁端部でゼロとし、３次式で立ち上げて補間する。３Ｄでも同様とし、それらをｚ方向端部でも３次補間して、いずれの端部でも垂直壁の値（ゼロ値）との連続性を確保する。
【０１１６】
次に、非圧縮の運動方程式の空間微分項（局所加速度項以外の項）を、その移流速度には速度要素の値を代入して、速度勾配には偏導関数要素の値を代入し、粘性項にも偏導関数要素の値を代入して、圧力勾配には圧力要素の値を代入して表した式を次に示す。
Ｈi≡Σj(βjβi,j−ν∂βi,j/∂ｘj）＋∂Ｐ/∂ｘi，   (４.１)
(i＝ｘ，ｙ，ｚ)，(ｊ＝ｘ，ｙ，ｚ)，（ｘx＝ｘ，ｘy＝ｙ，ｘz＝ｚ）
ＨiはＦiと同じであるが、式Ｆiはすべて速度要素のみでも表すことができるが、Ｆiの速度勾配に偏導関数要素を代入した場合に限り、特にＨiと記号する。
【０１１７】
次に、ステップ１の圧力ポアソン方程式のガレルキン法による重みつき残差式を示す。まず、(４.１)式の左辺で示せば、
∫ΩδＰ[∂Ｈi/∂ｘi]dΩ＝０   (４.２)
となる。(４.２)式に等価集中力を入れて陰的モデルで示せば、
∫Ωδ（∂Ｐ/∂ｘi ）[Ｈi]dΩ＝０   (４.３)
となる。(４.３)式を(４.１)式右辺で表せば、
∫Ωδ（∂Ｐ/∂ｘi ）[∂Ｐ/∂ｘi＋Σj(βjβi,j−ν∂βi,j/∂ｘj）] dΩ＝０   (４.４)
となる。(４.４)式は次の最小二乗法による極小条件式の二次形式となっている。
(1/2)∫Ω [δ(Ｈi)2]dΩ＝０   
よって、(４.４)式に、速度要素（初期値は境界値のみ）による移流速度と偏導関数要素（初期値はゼロ）による速度勾配を代入して、連立方程式を解いて圧力分布を求めることができる。
【０１１８】
次に、ステップ２の、速度要素（初期値は境界値のみ）による移流速度と偏導関数要素の値（初期値はゼロ）、並びにステップ１で求めた圧力分布を与えて、(４.１)式のＨiを求め、速度要素の加速度との和の極小条件式を設定し、速度要素のノード値を求める方法について説明する。釣り合い誤差の極小条件式は次に示す。
(1/2)∫Ω[δ(Req)2]dΩ＝０      (５)
(Req)＝∂βi/∂t＋Ｈi
(５)式の極小条件式は一般的な最小二乗法と同じ二次形式となるので、解法上の問題はない。よって、速度要素のノード値が求まる。
【０１１９】
次に、ステップ３の速度要素の補正について説明する。(５)式に代入したＨiの値は非圧縮条件を満足する値である。なぜなら、圧力ポアソン方程式は局所加速度の発散をゼロとして組立てた方程式だからである。したがって、(５)式を解いた結果の加速度も非圧縮となる。これを補正する意味は、数値誤差を取り除くことにある。圧縮性は速度ひずみテンソルの対角項のみで表される。よって、補正値は非対角項に影響を与えてはならない。そこで、非対角項の補正による影響を極小化する条件を制約条件式として、(５)式で求めた速度要素の発散の残差を極小化する。それら極小条件式を次に示す。
(1/2)∫Ω[δ(Req)2＋δ(Rcc)2]dΩ＝０      (５.1)
(Req)＝Σjδij ｛∂（Δβi＋βi)/∂xj)
(Rcc)＝（１−δij ）(∂Δβi/∂xj)
ただし、βiは(５)式で求めた値、Δβiは補正量（変数）、δijはクロネッカーのデルタである。(５.１)式を解いて求めた速度要素の補正量で速度要素をアップデートする。
βi  SYMBOL 223 \f "Wingdings" \s 12 βi＋Δβi
【０１２０】
次に、ステップ４の、各偏導関数要素の加速度と空間微分項Ｈiの各勾配の釣合い方程式から、速度要素の加速度の勾配とも等しくなくてはならないとする制約条件式の下で、最小二乗法により極小条件式を設定する手段について説明する。空間微分項Ｈiの勾配の釣合い方程式の残差ベクトルをReqと記し、制約条件式(速度要素の勾配との差) の残差ベクトルを時間ピッチΔtで除してRccと記す。ReqとRccの単位は同じになるように選んでいるので、最小二乗法が適用できる。ベクトルの二乗は内積をとるとする。
【０１２１】
最小二乗方程式：
(1/2)∫Ω[δ(Req)2＋δ(Rcc)2]dΩ＝０      (６)
(Req)＝∂βi,k/∂t＋∂Ｈi/∂xk
{βi,k} ≡{βx,x, βx,y, βx,z, βy,x, βy,y, βy,z, βz,x, βz,y, βz,z}
∂{Ｈi}/∂xk≡∂{Σj(βjβi,j−ν∂βi,j/∂ｘj）＋∂Ｐ/∂ｘi}/∂xk，
(i＝ｘ，ｙ，ｚ)，(ｊ＝ｘ，ｙ，ｚ)，（ｘx＝ｘ，ｘy＝ｙ，ｘz＝ｚ），
(k＝ｘ，ｙ，ｚ)
(Rcc)＝(1/Δt)(βi,j-∂βi/∂xj)
実際には、式(６)にさらに、要素間境界の等価集中力および基本境界での循環的対称モデルを加えて、準ニュートン法的に解く。
【０１２２】
次に、数値解法手段で、極小条件式を準ニュートン法で解く方法について説明する。準ニュートン法は定まった方法があるものではなく、ニュートン・ラプソン法に準ずる方法であればよい。ただし、反復計算の初期値が解の近傍に設定されている必要がある。本発明の実施例では、方程式のδ(Req)2の項にはガレルキン法を適用する。ガレルキン法は最小二乗法に準じた方程式をニュートン・ラプソン法で解く方法である。その式は、υを緩和係数（υ＝1）として、次の方程式となる。
【０１２３】
準ニュートン法方程式：
∫Ω[δ{βi,k} ?(Req) ＋(υ/2)δ(Rcc)2]dΩ
＋∫S[νΣjδβi,k?{∂βi,k/∂ｘj}]dＳ
＋∫Γ[νΣjδβi,k?{∂βi,k/∂ｘj}]dΓ
＋∫Γ[νΣjδβi,k?{∂βi,k/∂ｘj}]dΓ＝０      (７)
Ｓは要素間境界、Γは系の境界である。それらの積分項は要素間等価集中力モデルと系の境界の等価集中力モデルである。Γの積分項が２つあるのは、応力の対称性を入れるためである。これを循環的対称モデルと呼ぶ。
【０１２４】
循環的対称モデルについて説明する。基本境界では対称条件を入れて解く。これを図３（ａ）で説明する。陽的モデル化項は部分積分項の符号を替えて外力として機械的に加えればよいモデルである。対称かつ不連続な境界では、隣の系の等価集中力もあるのでこれを加えるモデルである。これは、図３(ｂ)の周期境界条件（離れたａ端とｃ端を連続させる条件）を入れた対象な系を解くことと同じである。基本境界のΓの積分項は、Ｃ0-連続な偏導関数要素では、ノード未知数は連立方程式に変数として残ることとなる。
【０１２５】
外力としての移動壁の値は、ステップ１で求めた基本変数要素{βi}の値を通じて、制約条件式Rccで偏導関数要素（渦度）に伝わる。従来の、非特許文献１や非特許文献６では図３(ｃ)の基本変数要素値の勾配を境界上で求めて、図３(ｄ)のように基本変数要素の境界値として与えていた。循環的対称モデルで解いた結果は図３(ｅ）のようになる。図３(ｄ)はアンダーシュートしているが、図３(ｅ)は、対称条件によってアンダーシュートしない。これが発明の循環的対称モデルの特長である。偏導関数要素法は循環的対称モデルをＣ0-連続な要素で実行する数値計算装置である。
【０１２６】
次に、繰返し手段について説明する。解析対象物の物理量を表すベクトル場をHelmholtz分解した変数で解析対象物の変化を規定する微分方程式と制約条件とから最小二乗法により求めた極小条件式を設定する方程式設定手段１と、解析対象物の境界条件を設定する境界条件設定手段２と、解析対象物の初期条件を設定する初期条件設定手段３とで組立てられた連立方程式を、極小条件式を準ニュートン法により有限要素法で解く数値計算手段41と、その解を補正する補正手段42とで解くことで、第１近似解を得る。次いで、第１近似解を初期値として、第２近似解を計算する。これを繰り返す。第２近似解と第１近似解の差が小さくなれば、収束したと判断して、繰返し計算を終了し、Δｔ後の速度の状態を、加速度から計算する。
【０１２７】
Δｔを小さく設定しているので、Δｔ後の速度の状態は、次のΔｔ後の初期値の資格がある。よって、時間を進めて計算し、所望の時間まで、または所望の速度分布の状態になるまで計算する。計算が完了すれば結果を出力する。
【０１２８】
図９は、Ｃ0-連続な偏導関数要素法で計算した結果の速度要素の値である。プロポーション１×１×６の３Ｄキャビティフローを、対称条件を入れず、格子数13×13×23（要素数12×12×22）として、計算した結果を図９に示す。等間隔の格子11×11×21の端をそれぞれ２つ割りして不等間隔格子を設定したものである。Ｚ方向は1/6に縮尺して作図している。図９(a-1)は、対称断面の流速プロファイルであり、○は非特許文献１の計算結果である。多少のずれは、格子が粗いことのほか、渦度のソレノイダル条件を直接的には考慮していないこと、圧力には偏導関数要素を適用していないこと、などが原因として挙げられる。速度の偏導関数要素およびそれへの循環的対称モデルの適用、および速度要素のソレノイダル条件の補正技法は適用して、安定的に作動することは確認できた。さらに、渦度のソレノイダル条件を、補正技法を適用して精度高く満足させること、圧力にも偏導関数要素を適用して、分解能を向上させること、など発明の技法を他の面へも応用することで、精度を高めることができる。図９(a-2)は圧力コンターであり、偏導関数要素の値を与えて計算して、定性的に正しい結果を得ることが確認できた。図９(ｂ)は速度の各断面のコンターであるが、対称な結果を得ることが確認できた。
【０１２９】
本発明の実施例３では、数値解析装置を４ステップから成るハイブリッド法とし、本発明の実施例２の循環的対称モデルをＣ0-連続な有限要素に適用することが特徴である。Ｃ0-連続な要素は基本変数とその偏導関数に適用し、前者を基本変数要素と呼び、後者を偏導関数要素と呼んだ。循環的対称モデルは偏導関数要素に適用した。従来のハイブリッド法が、速度要素の境界上での勾配を偏導関数要素の境界値として与えて固定するのに対し、循環的対称モデルでは偏導関数要素の境界上の値は変数として取り扱うところが、従来のハイブリッド法にない、新しい特徴である。
【０１３０】
本発明の実施例３は、解析対象物の物理量を表すベクトル場を基本変数（変位、速度）レベルでHelmholtz分解して変数表示し、解析対象物の変化を規定する微分方程式に渦度レベル（運動方程式の勾配のレベル）の高階の微分方程式を選び、分解した変数で循環的対称モデルを適用して表した微分方程式と緩和係数を乗じた制約条件とから最小二乗法により求めた極小条件式を設定し、解析対象物の境界条件と初期条件を設定し、極小条件式を準ニュートン法によりＣ0-連続な有限要素法で解く構成とした。緩和係数はυ＝１でも安定して解けることと、簡素な要素をハイブリッドすることで、循環的対称モデルによって、適合要素として要求される微分階数よりも高い、応力レベルの階数の微分にＣ0-連続な有限要素を適用して、境界ノード値を内的に決定できるようにし、高階の微分方程式を対象として解析できることで、よりミクロな変化を捉えることができるようになったことが、本発明の実施例３の大きな特徴である。
【産業上の利用可能性】
【０１３１】
本発明の数値解析装置は、流体や弾性体の動きを解析する数値解析装置として最適である。
【図面の簡単な説明】
【０１３２】
【図１】本発明の実施例における数値解析装置の機能ブロック図と動作手順を示す流れ図である。
【図２】本発明の実施例における数値解析装置の基本原理を説明する図である。
【図３】本発明の実施例における数値解析装置の循環的対称モデルと偏導関数要素法を説明する図である。
【図４】本発明の実施例における数値解析装置の修正なしの計算結果を示す図である。
【図５】本発明の実施例における数値解析装置の修正ありの計算結果を示す図である。
【図６】本発明の実施例における数値解析装置の計算結果を示す図である。
【図７】本発明の実施例における数値解析装置の計算結果を示す図である。
【図８】Ｃ1-連続な速度要素で計算した結果である。
【図９】Ｃ0-連続な速度要素で計算した結果である。
【図１０】従来の有限要素法を説明する図である。
【図１１】従来例の計算結果を示すグラフである。
【符号の説明】
【０１３３】
１  方程式設定手段
２  境界条件設定手段
３  初期条件設定手段
４  数値解法手段
５  出力手段
41  数値計算手段
42  補正手段
43  繰返し手段
【図１】
2007286801.tif 000003
【図２】
2007286801.tif 000004
【図３】
2007286801.tif 000005
【図４】
2007286801.tif 000006
【図５】
2007286801.tif 000007
【図６】
2007286801.tif 000008
【図７】
2007286801.tif 000009
【図８】
2007286801.tif 000010
【図９】
2007286801.tif 000011
【図１０】
2007286801.tif 000012
【図１１】
2007286801.tif 000013