(57)【要約】【課題】GSMAC有限要素法において、解析領域が複雑形状を有する場合や、流体と構造の相互作用が強い場合の非圧縮性流体解析を、低記憶容量の計算機で計算速度をできるだけ損なわずに高精度に行えるようにする。【解決手段】３角形要素と４角形要素を使用し、速度の補間関数として２次関数を使用し、圧力の補間関数として１次関数を使用する。３角形要素と４角形要素の混合メッシュを用いて計算することも可能である。離散化Poisson方程式を優対角近似し、圧力の節点の衛星要素を考慮して、速度と圧力を同時緩和法で収束するまで計算して、収束した結果について時間進行する。【選択図】図１
2006099549.tif 000002
【特許請求の範囲】
【請求項１】
GSMAC有限要素法により非圧縮性流体の流れを解析する数値計算方法おいて、任意の形状の要素を使用して解析対象をメッシュ化し、任意の次数の補間関数を使用して流体の速度と圧力を表わし、離散化Poisson方程式を優対角近似し、圧力の節点の衛星要素を考慮して速度と圧力を同時緩和法で収束するまで計算して、収束した結果について時間進行することを特徴とする数値計算方法。
【請求項２】
速度の補間関数として２次関数を用い、圧力の補間関数として１次関数を用いることを特徴とする請求項１記載の数値計算方法。
【請求項３】
前記要素として、４角形要素を使用し、速度の補間関数として２次関数を使用し、圧力の補間関数として１次関数を使用することを特徴とする請求項１記載の数値計算方法。
【請求項４】
前記要素として、３角形要素を使用し、速度の補間関数として２次関数を使用し、圧力の補間関数として１次関数を使用することを特徴とする請求項１記載の数値計算方法。
【請求項５】
前記要素として、４角形要素と３角形要素を組み合わせて使用し、速度の補間関数として２次関数を使用し、圧力の補間関数として１次関数を使用することを特徴とする請求項１記載の数値計算方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、GSMAC有限要素法による高次混合要素Poissonソルバーの数値計算手法に関し、特に、解析領域の形状が複雑である場合や、流体と構造の相互作用が強い場合における非圧縮性流体の動きを、低容量の計算機で高精度に解析するためのGSMAC有限要素法による高次混合要素Poissonソルバーの数値計算手法に関する。
【背景技術】
【０００２】
非圧縮性流体の非定常解析では、差分法のMAC(Marker and Cell)系解法(非特許文献１参照)から発展した分離型有限要素法がよく用いられる。分離型有限要素法では、運動方程式を時間についてのみ離散化した半離散化式から、非圧縮拘束条件を利用して分離式を求め、それを解いて近似解を求める。このとき、分離式中に圧力または修正速度ポテンシャルのPoisson方程式が導入される。このPoisson方程式を解くことによって、非圧縮拘束条件を満足する解を得ることができる。
【０００３】
GSMAC有限要素法(Generalized - Simplified MAC - FEM)(非特許文献２参照)は、差分法のHSMAC法(Highly Simplified MAC method)(非特許文献３参照)を拡張した分離型有限要素法である。GSMAC有限要素法では、優対角近似した離散化Poisson方程式を計算する、速度と圧力の同時緩和法(非特許文献２参照)が利用される。これによって、離散化Poisson方程式を解くための全体係数行列を作成する手順を省くことができる。また、要素係数行列を要素形状に依存する部分と要素形状に依存しない部分に分解する、係数行列の解析的表示を用いる(非特許文献２参照)。要素係数行列自身を記憶せずに、要素形状に依存しない部分のみを記憶し、要素形状に依存する部分は毎時間ステップで計算を行う。これにより、Gaussの数値積分のような計算時間を必要とせずに、記憶容量が低減できる。
【０００４】
流体解析で重要となるのは、高Reynolds数における対流項の計算である。rotating cone問題において、１次の補間関数を使う４角形要素を使用した場合と、２次の補間関数を使う４角形要素を使用した場合について、移流方程式の解の精度が比較されている(非特許文献４参照)。解析領域内の節点数が少なくても、２次の補間関数を使う４角形要素の方が、精度が良い。１次の補間関数を使う要素で同程度の精度を得るには、２次の補間関数を使う要素を利用する場合より、計算機記憶容量と計算時間が多く必要である。このことは、３角形要素でも同様である。
【０００５】
流体と構造の相互作用が強い場合に用いられる強連成法では、要素の適合条件を満たすためには、界面上における流体速度と構造変位の形状関数が一致する必要がある(非特許文献５参照)。構造変位に高次の形状関数を用いる場合には、流体速度にもそれを用いる必要がある。高次の形状関数を使う要素を用いる場合には、速度と圧力の組み合わせが重要となる。定常状態の流れ場解析では、Lagrange未定乗数法やペナルティ法に代表される混合型有限要素法で解く場合が多く、速度と圧力の補間関数の組み合わせが研究されている。これは、線形解析におけるinf-sup条件(非特許文献６参照)によって、圧力の次数を速度の次数より下げる必要があるからである。非定常状態の流れ場解析では、速度と圧力の補間関数の組み合わせについての研究例がほとんど見られず、速度に対して１次の関数の要素が一般に用いられている。
【非特許文献１】Harlow, F. H. and Welch, J. E., Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface, Physics of Fluids, 8-12, 1965, pp.2182-2189.
【非特許文献２】棚橋隆彦，「流れの有限要素法解析I」, 朝倉書店, 1997.
【非特許文献３】Hirt, C. W. and Cook, J. L., Calculating three-dimensional flows around structures and over rough terrain, Journal of Computational Physics, 10, 1972, pp.324-340.
【非特許文献４】鄭忠孝, 小川裕司, 棚橋隆彦, 有限要素法における移流項の高精度計算(双2次要素を用いたGSMACFEM),日本機械学会論文集, B編, 67-653, 2001, pp.1-8.
【非特許文献５】Ghattas, O. and Li, X., A variational finite element method for stationary nonlinear fluid-solid iteration, Journal of Computational Physics, 121, 1995, pp.347-356.
【非特許文献６】菊池文雄,「有限要素法の数理」, 培風館, 1994.
【非特許文献７】棚橋隆彦,「流れの有限要素法解析II」, 朝倉書店, 1997.
【非特許文献８】中山司, 森峰男, 剛体着水現象の数値解析のためのマーカー粒子法を併用した有限要素法, 日本機械学会論文集, B 編, 61-583, 1995, pp.169-176.
【非特許文献９】名古屋靖一郎, 非圧縮流体計算におけるソレノイダル空間への射影法, 日本応用数理学会論文誌, 10-3, 2000, pp.273-282.
【非特許文献１０】Ghia, U., Ghia, K. H. and Shin, C. T., High-Resolutions for incompressible flow using the Navier-Stokes equations and a multi-grid method, Journal of Computational Physics, 48, 1982, pp.387-411.
【発明の開示】
【発明が解決しようとする課題】
【０００６】
しかし、従来のGSMAC有限要素法では、用いる要素の種類が制限されているという問題がある。２次元では、４角形要素のみが用いられ、速度に１次の補間関数が用いられ、圧力に要素内一定の補間関数が用いられる。これは、離散化Poisson方程式を優対角近似する時に、４角形要素を仮定し、かつ要素内一定の圧力補間関数を仮定しているからである。
【０００７】
解析領域が複雑な形状の場合の非圧縮性流体解析では、４角形要素より３角形要素が有利である。混合型有限要素法では、速度に２次以上の補間関数を使う３角形要素が用いられる。しかし、分離型有限要素法では、速度に２次以上の補間関数を使う３角形要素を用いることができず、解析領域が複雑形状を有する場合の非圧縮性流体解析が困難であった。
【０００８】
流体と構造の相互作用が強い場合には、強連成法が用いられる。しかし、構造が曲げ変形する場合の解析に１次の形状関数を用いると、shear lockingが生じるために、２次以上の形状関数を用いる必要がある。さらに、界面上における流体速度と構造変位の形状関数を一致させる必要があるので、速度に２次以上の形状関数を用いる必要がある。しかし、従来の方法では、この条件で計算することはできなかった。
【０００９】
本発明の目的は、上記従来の問題を解決して、GSMAC有限要素法において、解析領域が複雑形状を有する場合や、流体と構造の相互作用が強い場合の非圧縮性流体解析を、低容量の計算機で計算速度をできるだけ損なわずに高精度に行えるようにすることである。
【課題を解決するための手段】
【００１０】
上記の課題を解決するために、本発明では、GSMAC有限要素法により非圧縮性流体の流れを解析する数値計算方法を、任意の形状の要素を使用して解析対象をメッシュ化し、任意の次数の補間関数を使用して流体の速度と圧力を表わし、離散化Poisson方程式を優対角近似し、圧力の節点の衛星要素を考慮して速度と圧力を同時緩和法で収束するまで計算して、収束した結果について時間進行する方法とした。
【発明の効果】
【００１１】
本発明では、上記のような計算方法としたことにより、以下のことが可能となる。要素形状を限定しないことで、複雑形状周りの解析に対して有効である３角形要素を用いて、GSMAC有限要素法で解析できる。特に、４角形要素のＱ2Ｑ1要素と３角形要素のＰ2+Ｐ1要素は、速度と圧力の同時緩和反復の収束性が良い補間関数の組みであり、高次要素を用いることによる計算時間の増加を少なくすることができる。さらに、３角形要素と４角形要素の混合メッシュを用いて計算することも可能になる。また、係数行列の解析的表示を用いることで、高次要素を用いても、要素係数行列を記憶した場合より記憶容量を少なくすることができる。
【発明を実施するための最良の形態】
【００１２】
以下、本発明を実施するための最良の形態について、図１〜図19を参照しながら詳細に説明する。
【実施例】
【００１３】
本発明の実施例は、GSMAC有限要素法において、３角形要素と４角形要素を使用し、速度の補間関数として２次関数を使用し、圧力の補間関数として１次関数を使用し、離散化Poisson方程式を優対角近似し、圧力の節点の衛星要素を考慮して速度と圧力を同時緩和法で計算して時間進行する数値計算方法である。
【００１４】
最初に、本発明の実施例における数値計算方法を従来方法と対比した図１を参照しながら、数値計算方法の概略を説明する。GSMAC有限要素法において、要素の形状と補間関数の次数を限定せずに、離散化Poisson方程式を優対角近似する。すなわち、要素形状と補間関数の制限を取り除いて、任意の形状の要素と任意の次数の補間関数を使用可能にする。例えば、２次元の場合であれば、３角形要素を使用可能にするとともに、速度について２次以上の補間関数を使用可能にし、圧力について１次以上の補間関数を使用可能にする。さらに、圧力の節点の衛星要素を考慮しながら、速度と圧力を同時緩和法で計算する。また、Poisson方程式の行列表現を経由しないで、係数行列の解析的表示を用いることで、高次要素を用いても、要素係数行列を記憶した場合より記憶容量を少なくして、計算速度をできるだけ損なわないようにする。
【００１５】
次に、基礎方程式について説明する。非圧縮Newton流体の基礎方程式は、連続の式、Cauchyの運動方程式、Newton流体の構成方程式、変形速度-速度こう配関係式
∇・ｖ＝０   (1a)
ρ(∂v/∂t)＋ρｖ・∇ｖ＝∇・Ｔ+ρｂ   (1b)
Ｔ＝−ｐＩ＋２μＤ   (1c)
Ｄ＝(1/2)(∇ｖ＋ｖ[∇←])   (1d)
である。ただし、∇は空間座標のナブラ、ｖは速度、ρは密度、∂/∂tは空間時間微分、ＴはCauchy応力テンソル、ｂは単位質量当たりの体積力、ｐは圧力、Ｉは恒等テンソル、μは粘性係数、Ｄは変形速度テンソルである。[∇←]は、後形（右形）ナブラである。本実施例では、上述の基礎方程式を、
∇・ｖ＝０   (2a)
ρ(∂v/∂t)＋ρｖ・∇ｖ＝−∇ｐ＋∇・τ＋ρｂ   (2b)
τ＝μ(∇ｖ＋ｖ[∇←])   (2c)
のように変形する。
【００１６】
次に、離散化計算法について説明する。基礎方程式をGSMAC有限要素法によって離散化する。時間に関しては、MAC系解法の時間進行に従って離散化する。また、空間に関しては、Galerkin有限要素法により離散化する。時間と空間に関して離散化した後に、離散化Poisson方程式を優対角近似する。
【００１７】
次に、MAC系解法の時間進行を説明する。基礎方程式において、速度を陽的に、圧力を陰的に、時間的に離散化すると、
∇・ｖn+1＝０   (3a)
ρ(vn+1−ｖn)/Δｔ＝−ρｖn・∇ｖn＋ρ∇・(Δｔ/2)(ｖnｖn)・∇ｖn
−∇ｐn+1＋∇τn＋ρｂn   (3b)
τn＝μ(∇ｖn＋ｖn[∇←])   (3c)
となる。ただし、右辺第２項は、BTD(Balancing Tensor Diffusivity)項(非特許文献７)であり、１次の陽的な時間進行を行う有限要素法解析において生じる数値振動を抑制する効果がある。なお、肩付きの記号は、冪指数ではなく、時間を示す添え字である。さらに、MAC系解法の時間進行より、式(3a)、式(3b)を
(予測子ステップの式)
ρ([ｖ~]−ｖn)/Δt＝−ρｖn・∇ｖn＋ρ∇・(Δｔ/2)(ｖnｖn) ・∇ｖn
−∇ｐn＋∇τn＋ρｂn   (4a)
(Poisson 方程式)
∇2φp＝ρ∇・[ｖ~]   (4b)
(修正子ステップの式)
ρvn+1＝ρ[ｖ~]−∇φp   (4c)
ｐn+1＝ｐn＋φp/Δt   (4d)
のように分離する。ただし、φpは修正速度ポテンシャルである。
【００１８】
次に、Galerkin有限要素法による離散化を説明する。得られた時間的な離散化式を空間的に離散化する。ｖ、ｂは形状関数Ｎi(１≦ｉ≦ｎN)で補間され、ｐは形状関数Ｎlp(1≦l≦ｎNp)で補間されるとする。要素Ωe(１≦ｅ≦ｎE)において、節点ｉが局所節点αに対応し、節点ｌが局所節点ａに対応する。各要素に対する離散化式を以下に示す。
【００１９】
(予測子ステップの式)
ρＭ-αβ([ｖ~]β−ｖβn)/Δｔ
＝−ρｖen・Ａαβｖβn−(Δｔ/２)(ｖenｖen):Ｄαβｖβn＋Ｃαbｐbn
−μ{(trＤαβ)Ｉ＋(Ｄαβ)T}・ｖβn＋ρＭ-αβｂβn
＋ρ∫ΓeＮα[n]・(Δt/２)(ｖnｖn)・∇ｖndΓ＋∫ΓeＮα[t]ndΓ   (5a)
【００２０】
(Poisson方程式)
Ｃβa・<∇φp>β＝ρＣβa・[ｖ~]β   (5b)
<∇φp>i＝−{Σe=1nE(Ｃαbφb)}/{Σe=1nE∫ΩeＮαdΩ}   (5c)
【００２１】
(修正子ステップの式)
ρＭ-αβｖn+1β＝ρＭ-αβ[ｖ~]β＋Ｃαbφpb−∫ΓeＮα(φp[n])dΓ   (5d)
ｐan+1＝ｐan＋φpa/Δt   (5e)
【００２２】
ただし，係数行列は
(質量行列)
Ｍαβ＝∫ΩeＮαＮβdΩ   (6a)
(移流行列)
Ａαβ＝∫ΩeＮα∇ＮβdΩ   (6b)
(こう配行列)
Ｃαb＝∫Ωe(∇Ｎα)ＮbpdΩ   (6c)
(拡散行列)
Ｄαβ＝∫Ωe∇Ｎα∇ＮβdΩ   (6d)
で表され，Ｍ-αβは質量の集中化
ρＭ-αβ＝ρ∫ΩeＮαdΩ   （for α＝β）
＝０   （for α≠β）   (7)
から得られる集中質量行列である。また、Γe＝∂Ωeは要素境界であり、[n]はその外向き法線ベクトルであり、[t]は応力ベクトルである。
【００２３】
次に、優対角近似を説明する。式(5c)を式(5b)に代入して、全要素についての和を考えると、離散化Poisson方程式
−Σe=1nE[Ｃβa・{Σe=1nE(Ｃβcφpc)/Σe=1nE∫ΩeＮβdΩ}]
＝Σe=1nE(ρＣβa・[ｖ~]β)   (8)
を得る。式(8)の左辺においてφplの係数に着目すると、優対角近似におけるPoisson方程式の係数λlは
Σe=1nE{Ｃβa・(Σe=1nEＣβa/Σe=1nE∫ΩeＮβdΩ)}
＝(Σe=1nE∫ΩeＮapdΩ)λl  (a: no sum)   (9)
より求まる。Ｎαが１次関数、Ｎapが要素内一定の関数の場合には、従来のGSMAC有限要素法のλl(非特許文献２参照)とは一致せず、中山らが剛体着水現象の有限要素法解析で用いたλl(非特許文献８参照)に一致する。
【００２４】
次に、速度と圧力の同時緩和法を説明する。優対角近似より、速度と圧力の同時緩和法を適用した各要素に対するPoisson方程式と修正子ステップの式は、
−(∫ΩeＮapdΩ)λa[φ~](k)pa＝ρＣβa・[ｖ~](k)β  (a: no sum)   (10a)
ρＭ-αβ[ｖ~](k+1)β＝ρＭ-αβ[ｖ~](k)β＋Ｃαb[φ~]pb(k)
−∫ΓeＮα([φ~]p(k)[n])dΓ   (10b)
[ｐ~](k+1)a＝[ｐ~](k)a＋[φ~](k)pa/Δｔ   (10c)
のように書き直される。差分法のHSMAC法における速度と圧力の同時緩和法が、離散化Poisson方程式のSOR法による反復計算と等価であることは証明されている(非特許文献９参照)。以下では、この証明を有限要素法に拡張することによって、本実施例で用いるGSMAC有限要素法においても、同様のことが成り立つことを示す。
【００２５】
まず、
(速度の節点平均値)
<ｆ>i≡Σe=1nE∫ΩeＮαｆdΩ/Σe=1nE∫ΩeＮαdΩ   (11a)
(圧力の節点平均値)
<ｇ>l≡Σe=1nE∫ΩeＮapｇdΩ/Σe=1nE∫ΩeＮapdΩ   (11b)
のように定義する。次に、式(10b)の全要素について和をとり、両辺に形状関数Ｎiを乗じると
[ｖ~](k+1)＝[ｖ~](k)−∇[φ~](k)p   (12)
となる。式(12)における両辺の発散を節点平均すると
<∇・[ｖ~](k+1)>l＝<∇・[ｖ~](k)>l-<∇・∇[φ~](k)p>l   (13)
となり、
<∇Ｎβ・[ｖ~](k+1)β>l＝<∇Ｎβ・[ｖ~](k)β>l−<∇Ｎβ・<∇[φ~](k)p>β>l   (14)
が成り立つ。
【００２６】
式(10a)の全要素についての和を
(−λl)[φ~](k)pl＝<∇Ｎβ・[ｖ~](k)β>l   (15)
のように表し、これを式(14)に代入すると
(−λl)[φ~](k+1)pl
＝−{<∇Ｎβ・<∇[φ~](k)p>β>l−(−λl)[φ~](k)pl}   (16)
のように変形できる。一方、式(8)は
<∇Ｎβ・<∇φp>β>l＝<∇Ｎβ・[ｖ~]β>l   (17)
のように表すことができる。式(17)の左辺におけるφplの係数が−λlであることを考慮して、式(17)にJacobi法を適用すると、
(−λl)φ(k+1)pl＝<∇Ｎβ・[ｖ~]β>l
−{<∇Ｎβ・<∇φ(k)p>β>l−(−λl)φ(k)pl}   (18)
となる。
【００２７】
ここで、φ(0)p＝０、φ(k+1)p−φ(k)p＝[φ~](k)pであるから、式(18)から式(16)を得ることができる。よって、同時緩和法がJacobi法による反復計算と等価になるようにλlを決定することが示される。圧力が要素内一定の関数で補間されている場合には、式(10)を要素ごとに順次計算すると、同時緩和法はGauss-Seidel法による反復計算と等価になる。さらに、式(10a)の右辺に加速係数を乗じると、SOR法による反復計算と等価になる。本手法では、図２に示す圧力の節点l（図の中央の節点）が有する衛星要素（ハッチング部分）を考慮することで、式(10)を、圧力の節点ごとに順次計算する。以上のGSMAC有限要素法による数値計算方法の概略手順の流れ図を、図３に示す。速度と圧力の同時緩和法による反復計算の流れ図を、図４に示す。
【００２８】
以下に、その具体的な計算手順を示す。
１．圧力の節点ｌ＝１に着目する。
２．圧力の節点ｌが有する衛星要素内の速度から、速度の発散の圧力の節点平均値Ｄ~(k+1)lを求める。
３．Ｄ~(k+1)lより、修正速度ポテンシャル[φ~](k+1)plを求める。
４．ｌ＝２,３,…,ｎpNのように、２．と３．の計算を繰り返す。
５．ｌ＝１,２,…,ｎpNについてＤ~(k+1)lを求め、その最大値max(1≦l≦ｎpN)Ｄ~(k+1)lを求める。
６．max(1≦l≦ｎpN)Ｄ~(k+1)lが、収束判定基準ε1より小さいならば、速度と圧力を時間進行する。
【００２９】
次に、計算方法の検証として、Reynolds数5000のCavity内強制対流を扱い、様々な種類の要素について精度、記憶容量、計算速度、同時緩和反復回数の評価を行った結果について説明する。Cavity内強制対流を十分解析でき、本実施例の離散化Poisson方程式の優対角近似の方法と、圧力の節点に対する速度と圧力の同時緩和法の妥当性が、以下のように確認できた。このとき、要素係数行列自身を記憶せずに、係数行列の解析的表示を利用することで、要素内節点数の増大による記憶容量の増大をできるだけ抑えた。
【００３０】
２次元cavity内強制対流の解析について説明する。検証問題として、Reynolds数5000の２次元cavity内強制対流を扱う。cavity内の中心線に沿った速度分布の精度に関して、Ghiaらの結果(非特許文献10参照)と比較する。図５に、解析モデルと境界条件を示す。ただし、Lrは代表長さ、Vrは代表速度である。このとき、Reynolds数はRe＝ρVrLr/μで定義される。また、初期条件は解析領域内でｖ＝０、ｐ＝０である。
【００３１】
解析領域は、図６と図７に示す４角形要素または３角形要素によって分割される。Ｑ0要素は、補間関数が定数の４角形要素である。Ｑ1要素は、１次の補間関数の４角形要素である。Ｑ2要素は、２次の補間関数の４角形要素である。Ｐ0要素は、補間関数が定数の３角形要素である。Ｐ1要素は、１次の補間関数の３角形要素である。Ｐ2+要素は、気泡関数を含む２次の補間関数の３角形要素である。速度ｖをＱ1要素で表わし、圧力ｐをＱ0要素で表わす場合の要素を、Ｑ1Ｑ0要素とよぶ。速度ｖをＰ2+要素で表わし、圧力ｐをＰ1要素で表わす場合の要素を、Ｐ2+Ｐ1要素とよぶ。その他についても同様である。Ｑ2要素の要素内節点がないＱ(8)2要素では、式(7)の対角成分のα＝１，２，３，４が負の値となる。Ｐ2+要素の要素内節点がないＰ2要素では、式(7)の対角成分のα＝１，２，３が零となる。そのため、本実施例では、Ｑ(8)2要素とＰ2要素を用いない。様々な要素を用いて、要素数が同じ場合の速度分布の精度、計算機記憶容量、計算時間、速度と圧力の同時緩和反復回数について比較する。
【００３２】
４角形要素のメッシュを用いた解析について説明する。４角形要素のメッシュを用いて解析する。４角形要素の組み合わせとして、従来のGSMAC有限要素法(非特許文献２参照)、Ｑ1Ｑ0要素、Ｑ1Ｑ1要素、Ｑ2Ｑ0要素、Ｑ2Ｑ1要素を使用した手法の５種類によって解析を行う。ただし、従来のGSMAC有限要素法とＱ1Ｑ0要素の違いは、離散化Poisson方程式を優対角近似する方法のみである。また、計算データを、表１に示す。ただし、√(Se)は要素面積の平方根である。要素数30×30の解析メッシュと要素数60×60の解析メッシュを図８に示す。
【表１】
2006099549.tif 000003
【００３３】
図８に示すように、cavity内中央から壁面へ向けてメッシュを等比数列的に生成する。同時緩和反復計算の収束性が悪くなるように公比を大きく設定する。cavity内中心線に沿った速度分布を、図９に示す。同じ空間精度であるために、Ｑ1Ｑ0要素の場合の速度は従来のGSMAC有限要素法の解と一致する。また、Ｑ1Ｑ1要素の場合の速度も従来のGSMAC有限要素法の解とほとんど一致するが、cavity内全域に圧力振動が生じる。その振動は要素分割を細かくしてもあまり抑制されないことから、Ｑ1Ｑ1要素による離散化式の形式が大きく影響していることがわかる。この不安定性の原因は、離散化Poisson方程式において隣の節点同士の関係が疎になるためである。差分法のMAC系解法では、この不安定性を生じさせないように、速度と圧力の格子点を一致させないスタッガード格子が用いられる。
【００３４】
Ｑ2Ｑ0要素の場合の速度は、用いた４角形要素の中で精度が最も悪い。Ｑ1Ｑ0要素に比べて、自由度数が多いにもかかわらず、中央での流れが小さくなる。この場合の圧力を見ると、cavity内中央での圧力差が非常に小さいことがわかる。Ｑ2Ｑ1要素の場合の速度は、図６に示す４角形要素の中で、最もGhiaの解に近いことがわかる。解の分布が最大値または最小値となる近くで精度の差が大きい。Ｑ2Ｑ1要素の場合には、Ｑ1Ｑ1要素の場合のような圧力振動も見られない。
【００３５】
図10と図11に、同時緩和反復回数の時刻歴を示す。表２に、ｔ＝０からｔ＝100までの同時緩和反復回数の総和を示す。図10と図11を見ると、Ｑ1Ｑ0要素の場合の反復回数は従来のGSMAC有限要素法と比べて収束前では少ないが、収束後では多いことがわかる。この収束性の違いは、優対角近似の方法の違いのみによって現れるものであり、メッシュを細かくすることによってその違いは顕著に現れる。図10と図11より、Ｑ2Ｑ1要素の場合の反復回数は収束前でも収束後でも比較的少ないことがわかる。表２より、Ｑ2Ｑ1要素の場合の反復回数の総和は４角形要素の中で最も少なく、メッシュが細かくなると他の要素との反復回数の差が顕著になる。
【表２】
2006099549.tif 000004
【００３６】
表３に、計算機記憶容量の最大値と計算時間を示す。要素数30×30のメッシュを用いると、Ｑ2Ｑ1要素の場合の計算は従来のGSMAC有限要素法の場合と比べて1.66倍の記憶容量と4.66倍の計算時間が必要である。要素数を４倍にすると、Ｑ2Ｑ1要素の計算は従来のGSMAC有限要素法の場合と比べて2.01倍の記憶容量と2.65倍の計算時間が必要である。メッシュを細かくすることによって計算時間の比が小さくなることは、同時緩和反復の収束性が大きく影響している。これは、計算が大規模になるほど計算時間に大きく影響する。Ｑ2Ｑ1要素を用いると、自由度が多いために精度の向上が見られ、同時緩和反復の収束性も良い。Ｑ2Ｑ1要素は速度と圧力の補間関数の組み合わせが良い要素であると考えられる。
【表３】
2006099549.tif 000005
【００３７】
３角形要素のメッシュを用いた解析について説明する。３角形要素のメッシュを用いて解析する。３角形要素の組み合わせとして、Ｐ1Ｐ1要素、Ｐ2+Ｐ0要素、Ｐ2+Ｐ1要素を使用した手法の３種類によって解析を行う。また、計算データを、表４に示す。要素数30×30×２の解析メッシュと要素数60×60×２の解析メッシュを、図12に示す。これらのメッシュは、図８に示すメッシュにおける４角形要素を３角形要素に等分したものである。cavity内中心線に沿った速度分布を、図13に示す。Ｐ1Ｐ1要素の場合の速度では、ｘ＝0.5上で数値的な振動が見られる。また、Ｐ1Ｐ1要素の場合では、cavity内全域で圧力振動が生じる。Ｐ2+Ｐ0要素の場合の速度はcavity内中央で小さく、精度が悪い。Ｐ2+Ｐ0要素の場合には、Ｑ2Ｑ0要素の場合と同様に、cavity内中央で圧力差が非常に小さい。Ｐ2+Ｐ1要素の場合の速度は、用いた３角形要素の組み合わせの中で最もGhiaらの解に近い。特に、中心線上の速度分布が最大値または最小値となる近くでＰ1+Ｐ1要素の解の精度が良いことがわかる。Ｐ2+Ｐ1要素の場合、要素数が30×30×２のメッシュでは圧力振動が生じる部分があるが、要素数が60×60×２のメッシュではその圧力振動は非常に小さくなる。
【表４】
2006099549.tif 000006
【００３８】
図14と図15に、同時緩和反復回数の時刻歴を示す。表５に、ｔ＝０からｔ＝100までの同時緩和反復回数の総和を示す。Ｐ2+Ｐ1要素の反復回数は、収束前でも収束後でも最も少ない。表２より、Ｐ2+Ｐ1要素の反復回数の総和は３角形要素の中で最も少ない。表６に、計算機記憶容量の最大値と計算時間を示す。Ｐ1Ｐ1要素の場合と比べて，Ｐ2+Ｐ1要素の場合の記憶容量として要素数30×30×２のメッシュを用いた計算では1.53倍、要素数60×60×２のメッシュを用いた計算では1.71倍が必要である。しかし、計算時間として要素数30×30×２を用いた計算では0.21倍、要素数60×60×２を用いた計算では0.05倍しか必要でない。Ｐ2+Ｐ1要素の場合の計算では、メッシュを細かくするとその計算時間の比が非常に小さくなる。解の精度と同時緩和反復の収束性を考慮すると、Ｐ2+Ｐ1要素は速度と圧力の補間関数の組み合わせが良い要素であると考えられる。
【表５】
2006099549.tif 000007
【表６】
2006099549.tif 000008
【００３９】
混合メッシュを用いた解析について説明する。４角形要素と３角形要素の混合メッシュを用いて解析する。同時緩和反復の収束性と要素の適合条件を考慮して、４角形のＱ2Ｑ1要素と３角形のＰ2+Ｐ1要素によって解析を行う。計算データを、表７に示す。要素数1350の解析メッシュと要素数5400の解析メッシュを図16に示す。メッシュは、４角形要素と３角形要素の接する中心線で速度分布を検討できるように生成される。cavity内中心線に沿った速度分布を図17に示し、圧力分布を図18に示す。図17を見ると、混合メッシュを用いても速度の分布は数値的に振動することもなく、Ｑ2Ｑ1要素の解とＰ2+Ｐ1要素の解と比べて速度の精度も損なわれないことがわかる。図18を見ると、要素数が1350のメッシュでは圧力振動が生じる部分があるが、要素数が5400のメッシュではその圧力振動は非常に小さくなる。図19に同時緩和反復回数の時刻歴を示す。図19より、反復回数は収束前でも収束後でも少なく、混合メッシュを用いたことによって反復回数の増加も起こらないことがわかる。
【表７】
2006099549.tif 000009
【００４０】
高次要素を適用したGSMAC有限要素法によって、cavity内強制対流を十分解析できる。要素係数行列自身を記憶せずに、係数行列の解析的表示を利用することで、要素内節点数の増大による記憶容量の増大を抑えることができる。
【００４１】
４角形要素ではＱ2Ｑ1要素が、同時緩和反復の収束性が良い。３角形要素では、Ｐ2+Ｐ1要素が、同時緩和反復の収束性が良い。要素数が増えると、記憶容量の比は増えるが、同時緩和反復の回数が計算速度に影響するために、計算時間の比は減少した。同じ精度の解を得る場合、Ｐ2+Ｐ1要素は、３角形要素の中で最も計算時間が少ない。４角形のＱ2Ｑ1要素と３角形のＰ2+Ｐ1要素による混合メッシュを用いることも可能である。
【００４２】
上記のように、本発明の実施例では、GSMAC有限要素法における数値計算方法を、３角形要素と４角形要素を使用し、速度の補間関数として２次関数を使用し、圧力の補間関数として１次関数を使用し、離散化Poisson方程式を優対角近似し、圧力の節点の衛星要素を考慮して速度と圧力を同時緩和法で計算して時間進行する方法としたので、対象領域が複雑形状の場合や、流体と構造の相互作用が強い場合における非圧縮性流体の流れを、低記憶容量の計算機で計算速度をできるだけ損なわずに高精度に解析できる。
【産業上の利用可能性】
【００４３】
本発明のGSMAC有限要素法における数値計算方法は、解析領域が複雑形状を有する場合の非圧縮性流体解析に最適である。また、流体と構造の相互作用が強い場合に強連成法で解析する方法としても適当である。
【図面の簡単な説明】
【００４４】
【図１】本発明の実施例における数値計算方法を従来方法と比較した説明図である。
【図２】本発明の実施例における数値計算方法で使用する４角形要素と３角形要素における衛星要素を示す図である。
【図３】本発明の実施例におけるGSMAC有限要素法の数値計算方法の流れ図である。
【図４】本発明の実施例における数値計算方法で用いる速度と圧力の同時緩和法による反復計算の流れ図である。
【図５】本発明の実施例における数値計算方法の解析例の解析モデルと境界条件を示す図である。
【図６】本発明の実施例における数値計算方法で用いる解析領域を分割する４角形要素を示す図である。
【図７】本発明の実施例における数値計算方法で用いる解析領域を分割する３角形要素を示す図である。
【図８】本発明の実施例における数値計算方法の解析例での要素数30×30の解析メッシュと要素数60×60の解析メッシュを示す図である。
【図９】本発明の実施例における数値計算方法の解析例でのcavity内中心線に沿った速度分布を示す図である。
【図１０】本発明の実施例における数値計算方法の解析例での同時緩和反復回数の時刻歴を示す図である。
【図１１】本発明の実施例における数値計算方法の解析例での同時緩和反復回数の時刻歴を示す図である。
【図１２】本発明の実施例における数値計算方法の解析例での要素数30×30×２の解析メッシュと要素数60×60×２の解析メッシュを示す図である。
【図１３】本発明の実施例における数値計算方法の解析例でのcavity内中心線に沿った速度分布を示す図である。
【図１４】本発明の実施例における数値計算方法の解析例での同時緩和反復回数の時刻歴を示す図である。
【図１５】本発明の実施例における数値計算方法の解析例での同時緩和反復回数の時刻歴を示す図である。
【図１６】本発明の実施例における数値計算方法の解析例での要素数1350の解析メッシュと要素数5400の解析メッシュを示す図である。
【図１７】本発明の実施例における数値計算方法の解析例でのcavity内中心線に沿った速度分布を示す図である。
【図１８】本発明の実施例における数値計算方法の解析例での圧力分布を示す図である。
【図１９】本発明の実施例における数値計算方法の解析例での同時緩和反復回数の時刻歴を示す図である。
【図１】
2006099549.tif 000010
【図２】
2006099549.tif 000011
【図３】
2006099549.tif 000012
【図４】
2006099549.tif 000013
【図５】
2006099549.tif 000014
【図６】
2006099549.tif 000015
【図７】
2006099549.tif 000016
【図８】
2006099549.tif 000017
【図９】
2006099549.tif 000018
【図１０】
2006099549.tif 000019
【図１１】
2006099549.tif 000020
【図１２】
2006099549.tif 000021
【図１３】
2006099549.tif 000022
【図１４】
2006099549.tif 000023
【図１５】
2006099549.tif 000024
【図１６】
2006099549.tif 000025
【図１７】
2006099549.tif 000026
【図１８】
2006099549.tif 000027
【図１９】
2006099549.tif 000028