JP4604260B2

JP4604260B2 - ＧＳＭＡＣ有限要素法による高次混合要素Ｐｏｉｓｓｏｎソルバーの数値計算手法

Info

Publication number: JP4604260B2
Application number: JP2004286430A
Authority: JP
Inventors: ▲隆▼彦棚橋; 橋本　　学
Original assignee: Keio University
Current assignee: Keio University
Priority date: 2004-09-30
Filing date: 2004-09-30
Publication date: 2011-01-05
Anticipated expiration: 2024-09-30
Also published as: JP2006099549A

Description

本発明は、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有限要素法で解析できる。特に、４角形要素のＱ₂Ｑ₁要素と３角形要素のＰ₂ ⁺Ｐ₁要素は、速度と圧力の同時緩和反復の収束性が良い補間関数の組みであり、高次要素を用いることによる計算時間の増加を少なくすることができる。さらに、３角形要素と４角形要素の混合メッシュを用いて計算することも可能になる。また、係数行列の解析的表示を用いることで、高次要素を用いても、要素係数行列を記憶した場合より記憶容量を少なくすることができる。

以下、本発明を実施するための最良の形態について、図１〜図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系解法の時間進行を説明する。基礎方程式において、速度を陽的に、圧力を陰的に、時間的に離散化すると、
∇・ｖⁿ⁺¹＝０ (3a)
ρ(vⁿ⁺¹−ｖⁿ)/Δｔ＝−ρｖⁿ・∇ｖⁿ＋ρ∇・(Δｔ/2)(ｖⁿｖⁿ)・∇ｖⁿ
−∇ｐⁿ⁺¹＋∇τⁿ＋ρｂⁿ (3b)
τⁿ＝μ(∇ｖⁿ＋ｖⁿ[∇^←]) (3c)
となる。ただし、右辺第２項は、BTD(Balancing Tensor Diffusivity)項(非特許文献７)であり、１次の陽的な時間進行を行う有限要素法解析において生じる数値振動を抑制する効果がある。なお、肩付きの記号は、冪指数ではなく、時間を示す添え字である。さらに、MAC系解法の時間進行より、式(3a)、式(3b)を
(予測子ステップの式)
ρ([ｖ~]−ｖⁿ)/Δt＝−ρｖⁿ・∇ｖⁿ＋ρ∇・(Δｔ/2)(ｖⁿｖⁿ) ・∇ｖⁿ
−∇ｐⁿ＋∇τⁿ＋ρｂⁿ (4a)
(Poisson 方程式)
∇²φ_p＝ρ∇・[ｖ~] (4b)
(修正子ステップの式)
ρvⁿ⁺¹＝ρ[ｖ~]−∇φ_p (4c)
ｐⁿ⁺¹＝ｐⁿ＋φ_p/Δt (4d)
のように分離する。ただし、φ_pは修正速度ポテンシャルである。

次に、Galerkin有限要素法による離散化を説明する。得られた時間的な離散化式を空間的に離散化する。ｖ、ｂは形状関数Ｎ_i(１≦ｉ≦ｎ_N)で補間され、ｐは形状関数Ｎ_l ^p(1≦l≦ｎ_N ^p)で補間されるとする。要素Ω_e(１≦ｅ≦ｎ_E)において、節点ｉが局所節点αに対応し、節点ｌが局所節点ａに対応する。各要素に対する離散化式を以下に示す。

(予測子ステップの式)
ρＭ^- _αβ([ｖ~]_β−ｖ_β ⁿ)/Δｔ
＝−ρｖ_e ⁿ・Ａ_αβｖ_β ⁿ−(Δｔ/２)(ｖ_e ⁿｖ_e ⁿ):Ｄ_αβｖ_β ⁿ＋Ｃ_αbｐ_b ⁿ
−μ{(trＤ_αβ)Ｉ＋(Ｄ_αβ)^T}・ｖ_β ⁿ＋ρＭ^- _αβｂ_β ⁿ
＋ρ∫_ΓeＮ_α[n]・(Δt/２)(ｖⁿｖⁿ)・∇ｖⁿdΓ＋∫_ΓeＮ_α[t]ⁿdΓ (5a)

(Poisson方程式)
Ｃ_βa・<∇φ_p>_β＝ρＣ_βa・[ｖ~]_β (5b)
<∇φ_p>_i＝−{Σ_e=1 ^nE(Ｃ_αbφ_b)}/{Σ_e=1 ^nE∫_ΩeＮ_αdΩ} (5c)

(修正子ステップの式)
ρＭ^- _αβｖⁿ⁺¹ _β＝ρＭ^- _αβ[ｖ~]_β＋Ｃ_αbφ_pb−∫_ΓeＮ_α(φ_p[n])dΓ (5d)
ｐ_a ⁿ⁺¹＝ｐ_a ⁿ＋φ_pa/Δt (5e)

ただし，係数行列は
(質量行列)
Ｍ_αβ＝∫_ΩeＮ_αＮ_βdΩ (6a)
(移流行列)
Ａ_αβ＝∫_ΩeＮ_α∇Ｎ_βdΩ (6b)
(こう配行列)
Ｃ_αb＝∫_Ωe(∇Ｎ_α)Ｎ_b ^pdΩ (6c)
(拡散行列)
Ｄ_αβ＝∫_Ωe∇Ｎ_α∇Ｎ_βdΩ (6d)
で表され，Ｍ^- _αβは質量の集中化
ρＭ^- _αβ＝ρ∫_ΩeＮ_αdΩ （for α＝β）
＝０（for α≠β） (7)
から得られる集中質量行列である。また、Γ_e＝∂Ω_eは要素境界であり、[n]はその外向き法線ベクトルであり、[t]は応力ベクトルである。

次に、優対角近似を説明する。式(5c)を式(5b)に代入して、全要素についての和を考えると、離散化Poisson方程式
−Σ_e=1 ^nE[Ｃ_βa・{Σ_e=1 ^nE(Ｃ_βcφ_pc)/Σ_e=1 ^nE∫_ΩeＮ_βdΩ}]
＝Σ_e=1 ^nE(ρＣ_βa・[ｖ~]_β) (8)
を得る。式(8)の左辺においてφ_plの係数に着目すると、優対角近似におけるPoisson方程式の係数λ_lは
Σ_e=1 ^nE{Ｃ_βa・(Σ_e=1 ^nEＣ_βa/Σ_e=1 ^nE∫_ΩeＮ_βdΩ)}
＝(Σ_e=1 ^nE∫_ΩeＮ_a ^pdΩ)λ_l (a: no sum) (9)
より求まる。Ｎ_αが１次関数、Ｎ_a ^pが要素内一定の関数の場合には、従来のGSMAC有限要素法のλ_l(非特許文献２参照)とは一致せず、中山らが剛体着水現象の有限要素法解析で用いたλ_l(非特許文献８参照)に一致する。

次に、速度と圧力の同時緩和法を説明する。優対角近似より、速度と圧力の同時緩和法を適用した各要素に対するPoisson方程式と修正子ステップの式は、
−(∫_ΩeＮ_a ^pdΩ)λ_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=1 ^nE∫_ΩeＮ_αｆdΩ/Σ_e=1 ^nE∫_ΩeＮ_αdΩ (11a)
(圧力の節点平均値)
<ｇ>_l≡Σ_e=1 ^nE∫_ΩeＮ_a ^pｇdΩ/Σ_e=1 ^nE∫_ΩeＮ_a ^pdΩ (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)
となる。

ここで、φ⁽⁰⁾ _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を求める。
４．ｌ＝２,３,…,ｎ^p _Nのように、２．と３．の計算を繰り返す。
５．ｌ＝１,２,…,ｎ^p _NについてＤ~^(k+1) _lを求め、その最大値max(1≦l≦ｎ^p _N)Ｄ~^(k+1) _lを求める。
６．max(1≦l≦ｎ^p _N)Ｄ~^(k+1) _lが、収束判定基準ε₁より小さいならば、速度と圧力を時間進行する。

次に、計算方法の検証として、Reynolds数5000のCavity内強制対流を扱い、様々な種類の要素について精度、記憶容量、計算速度、同時緩和反復回数の評価を行った結果について説明する。Cavity内強制対流を十分解析でき、本実施例の離散化Poisson方程式の優対角近似の方法と、圧力の節点に対する速度と圧力の同時緩和法の妥当性が、以下のように確認できた。このとき、要素係数行列自身を記憶せずに、係数行列の解析的表示を利用することで、要素内節点数の増大による記憶容量の増大をできるだけ抑えた。

２次元cavity内強制対流の解析について説明する。検証問題として、Reynolds数5000の２次元cavity内強制対流を扱う。cavity内の中心線に沿った速度分布の精度に関して、Ghiaらの結果(非特許文献10参照)と比較する。図５に、解析モデルと境界条件を示す。ただし、Lrは代表長さ、Vrは代表速度である。このとき、Reynolds数はRe＝ρVrLr/μで定義される。また、初期条件は解析領域内でｖ＝０、ｐ＝０である。

解析領域は、図６と図７に示す４角形要素または３角形要素によって分割される。Ｑ₀要素は、補間関数が定数の４角形要素である。Ｑ₁要素は、１次の補間関数の４角形要素である。Ｑ₂要素は、２次の補間関数の４角形要素である。Ｐ₀要素は、補間関数が定数の３角形要素である。Ｐ₁要素は、１次の補間関数の３角形要素である。Ｐ₂ ⁺要素は、気泡関数を含む２次の補間関数の３角形要素である。速度ｖをＱ₁要素で表わし、圧力ｐをＱ₀要素で表わす場合の要素を、Ｑ₁Ｑ₀要素とよぶ。速度ｖをＰ₂ ⁺要素で表わし、圧力ｐをＰ₁要素で表わす場合の要素を、Ｐ₂ ⁺Ｐ₁要素とよぶ。その他についても同様である。Ｑ₂要素の要素内節点がないＱ⁽⁸⁾ ₂要素では、式(7)の対角成分のα＝１，２，３，４が負の値となる。Ｐ₂ ⁺要素の要素内節点がないＰ₂要素では、式(7)の対角成分のα＝１，２，３が零となる。そのため、本実施例では、Ｑ⁽⁸⁾ ₂要素とＰ₂要素を用いない。様々な要素を用いて、要素数が同じ場合の速度分布の精度、計算機記憶容量、計算時間、速度と圧力の同時緩和反復回数について比較する。

４角形要素のメッシュを用いた解析について説明する。４角形要素のメッシュを用いて解析する。４角形要素の組み合わせとして、従来のGSMAC有限要素法(非特許文献２参照)、Ｑ₁Ｑ₀要素、Ｑ₁Ｑ₁要素、Ｑ₂Ｑ₀要素、Ｑ₂Ｑ₁要素を使用した手法の５種類によって解析を行う。ただし、従来のGSMAC有限要素法とＱ₁Ｑ₀要素の違いは、離散化Poisson方程式を優対角近似する方法のみである。また、計算データを、表１に示す。ただし、√(Se)は要素面積の平方根である。要素数30×30の解析メッシュと要素数60×60の解析メッシュを図８に示す。

図８に示すように、cavity内中央から壁面へ向けてメッシュを等比数列的に生成する。同時緩和反復計算の収束性が悪くなるように公比を大きく設定する。cavity内中心線に沿った速度分布を、図９に示す。同じ空間精度であるために、Ｑ₁Ｑ₀要素の場合の速度は従来のGSMAC有限要素法の解と一致する。また、Ｑ₁Ｑ₁要素の場合の速度も従来のGSMAC有限要素法の解とほとんど一致するが、cavity内全域に圧力振動が生じる。その振動は要素分割を細かくしてもあまり抑制されないことから、Ｑ₁Ｑ₁要素による離散化式の形式が大きく影響していることがわかる。この不安定性の原因は、離散化Poisson方程式において隣の節点同士の関係が疎になるためである。差分法のMAC系解法では、この不安定性を生じさせないように、速度と圧力の格子点を一致させないスタッガード格子が用いられる。

Ｑ₂Ｑ₀要素の場合の速度は、用いた４角形要素の中で精度が最も悪い。Ｑ₁Ｑ₀要素に比べて、自由度数が多いにもかかわらず、中央での流れが小さくなる。この場合の圧力を見ると、cavity内中央での圧力差が非常に小さいことがわかる。Ｑ₂Ｑ₁要素の場合の速度は、図６に示す４角形要素の中で、最もGhiaの解に近いことがわかる。解の分布が最大値または最小値となる近くで精度の差が大きい。Ｑ₂Ｑ₁要素の場合には、Ｑ₁Ｑ₁要素の場合のような圧力振動も見られない。

図10と図11に、同時緩和反復回数の時刻歴を示す。表２に、ｔ＝０からｔ＝100までの同時緩和反復回数の総和を示す。図10と図11を見ると、Ｑ₁Ｑ₀要素の場合の反復回数は従来のGSMAC有限要素法と比べて収束前では少ないが、収束後では多いことがわかる。この収束性の違いは、優対角近似の方法の違いのみによって現れるものであり、メッシュを細かくすることによってその違いは顕著に現れる。図10と図11より、Ｑ₂Ｑ₁要素の場合の反復回数は収束前でも収束後でも比較的少ないことがわかる。表２より、Ｑ₂Ｑ₁要素の場合の反復回数の総和は４角形要素の中で最も少なく、メッシュが細かくなると他の要素との反復回数の差が顕著になる。

表３に、計算機記憶容量の最大値と計算時間を示す。要素数30×30のメッシュを用いると、Ｑ₂Ｑ₁要素の場合の計算は従来のGSMAC有限要素法の場合と比べて1.66倍の記憶容量と4.66倍の計算時間が必要である。要素数を４倍にすると、Ｑ₂Ｑ₁要素の計算は従来のGSMAC有限要素法の場合と比べて2.01倍の記憶容量と2.65倍の計算時間が必要である。メッシュを細かくすることによって計算時間の比が小さくなることは、同時緩和反復の収束性が大きく影響している。これは、計算が大規模になるほど計算時間に大きく影響する。Ｑ₂Ｑ₁要素を用いると、自由度が多いために精度の向上が見られ、同時緩和反復の収束性も良い。Ｑ₂Ｑ₁要素は速度と圧力の補間関数の組み合わせが良い要素であると考えられる。

３角形要素のメッシュを用いた解析について説明する。３角形要素のメッシュを用いて解析する。３角形要素の組み合わせとして、Ｐ₁Ｐ₁要素、Ｐ₂ ⁺Ｐ₀要素、Ｐ₂ ⁺Ｐ₁要素を使用した手法の３種類によって解析を行う。また、計算データを、表４に示す。要素数30×30×２の解析メッシュと要素数60×60×２の解析メッシュを、図12に示す。これらのメッシュは、図８に示すメッシュにおける４角形要素を３角形要素に等分したものである。cavity内中心線に沿った速度分布を、図13に示す。Ｐ₁Ｐ₁要素の場合の速度では、ｘ＝0.5上で数値的な振動が見られる。また、Ｐ₁Ｐ₁要素の場合では、cavity内全域で圧力振動が生じる。Ｐ₂ ⁺Ｐ₀要素の場合の速度はcavity内中央で小さく、精度が悪い。Ｐ₂ ⁺Ｐ₀要素の場合には、Ｑ₂Ｑ₀要素の場合と同様に、cavity内中央で圧力差が非常に小さい。Ｐ₂ ⁺Ｐ₁要素の場合の速度は、用いた３角形要素の組み合わせの中で最もGhiaらの解に近い。特に、中心線上の速度分布が最大値または最小値となる近くでＰ₁ ⁺Ｐ₁要素の解の精度が良いことがわかる。Ｐ₂ ⁺Ｐ₁要素の場合、要素数が30×30×２のメッシュでは圧力振動が生じる部分があるが、要素数が60×60×２のメッシュではその圧力振動は非常に小さくなる。

図14と図15に、同時緩和反復回数の時刻歴を示す。表５に、ｔ＝０からｔ＝100までの同時緩和反復回数の総和を示す。Ｐ₂ ⁺Ｐ₁要素の反復回数は、収束前でも収束後でも最も少ない。表２より、Ｐ₂ ⁺Ｐ₁要素の反復回数の総和は３角形要素の中で最も少ない。表６に、計算機記憶容量の最大値と計算時間を示す。Ｐ₁Ｐ₁要素の場合と比べて，Ｐ₂ ⁺Ｐ₁要素の場合の記憶容量として要素数30×30×２のメッシュを用いた計算では1.53倍、要素数60×60×２のメッシュを用いた計算では1.71倍が必要である。しかし、計算時間として要素数30×30×２を用いた計算では0.21倍、要素数60×60×２を用いた計算では0.05倍しか必要でない。Ｐ₂ ⁺Ｐ₁要素の場合の計算では、メッシュを細かくするとその計算時間の比が非常に小さくなる。解の精度と同時緩和反復の収束性を考慮すると、Ｐ₂ ⁺Ｐ₁要素は速度と圧力の補間関数の組み合わせが良い要素であると考えられる。

混合メッシュを用いた解析について説明する。４角形要素と３角形要素の混合メッシュを用いて解析する。同時緩和反復の収束性と要素の適合条件を考慮して、４角形のＱ₂Ｑ₁要素と３角形のＰ₂ ⁺Ｐ₁要素によって解析を行う。計算データを、表７に示す。要素数1350の解析メッシュと要素数5400の解析メッシュを図16に示す。メッシュは、４角形要素と３角形要素の接する中心線で速度分布を検討できるように生成される。cavity内中心線に沿った速度分布を図17に示し、圧力分布を図18に示す。図17を見ると、混合メッシュを用いても速度の分布は数値的に振動することもなく、Ｑ₂Ｑ₁要素の解とＰ₂ ⁺Ｐ₁要素の解と比べて速度の精度も損なわれないことがわかる。図18を見ると、要素数が1350のメッシュでは圧力振動が生じる部分があるが、要素数が5400のメッシュではその圧力振動は非常に小さくなる。図19に同時緩和反復回数の時刻歴を示す。図19より、反復回数は収束前でも収束後でも少なく、混合メッシュを用いたことによって反復回数の増加も起こらないことがわかる。

高次要素を適用したGSMAC有限要素法によって、cavity内強制対流を十分解析できる。要素係数行列自身を記憶せずに、係数行列の解析的表示を利用することで、要素内節点数の増大による記憶容量の増大を抑えることができる。

４角形要素ではＱ₂Ｑ₁要素が、同時緩和反復の収束性が良い。３角形要素では、Ｐ₂ ⁺Ｐ₁要素が、同時緩和反復の収束性が良い。要素数が増えると、記憶容量の比は増えるが、同時緩和反復の回数が計算速度に影響するために、計算時間の比は減少した。同じ精度の解を得る場合、Ｐ₂ ⁺Ｐ₁要素は、３角形要素の中で最も計算時間が少ない。４角形のＱ₂Ｑ₁要素と３角形のＰ₂ ⁺Ｐ₁要素による混合メッシュを用いることも可能である。

上記のように、本発明の実施例では、GSMAC有限要素法における数値計算方法を、３角形要素と４角形要素を使用し、速度の補間関数として２次関数を使用し、圧力の補間関数として１次関数を使用し、離散化Poisson方程式を優対角近似し、圧力の節点の衛星要素を考慮して速度と圧力を同時緩和法で計算して時間進行する方法としたので、対象領域が複雑形状の場合や、流体と構造の相互作用が強い場合における非圧縮性流体の流れを、低記憶容量の計算機で計算速度をできるだけ損なわずに高精度に解析できる。

本発明のGSMAC有限要素法における数値計算方法は、解析領域が複雑形状を有する場合の非圧縮性流体解析に最適である。また、流体と構造の相互作用が強い場合に強連成法で解析する方法としても適当である。

本発明の実施例における数値計算方法を従来方法と比較した説明図である。本発明の実施例における数値計算方法で使用する４角形要素と３角形要素における衛星要素を示す図である。本発明の実施例におけるGSMAC有限要素法の数値計算方法の流れ図である。本発明の実施例における数値計算方法で用いる速度と圧力の同時緩和法による反復計算の流れ図である。本発明の実施例における数値計算方法の解析例の解析モデルと境界条件を示す図である。本発明の実施例における数値計算方法で用いる解析領域を分割する４角形要素を示す図である。本発明の実施例における数値計算方法で用いる解析領域を分割する３角形要素を示す図である。本発明の実施例における数値計算方法の解析例での要素数30×30の解析メッシュと要素数60×60の解析メッシュを示す図である。本発明の実施例における数値計算方法の解析例でのcavity内中心線に沿った速度分布を示す図である。本発明の実施例における数値計算方法の解析例での同時緩和反復回数の時刻歴を示す図である。本発明の実施例における数値計算方法の解析例での同時緩和反復回数の時刻歴を示す図である。本発明の実施例における数値計算方法の解析例での要素数30×30×２の解析メッシュと要素数60×60×２の解析メッシュを示す図である。本発明の実施例における数値計算方法の解析例でのcavity内中心線に沿った速度分布を示す図である。本発明の実施例における数値計算方法の解析例での同時緩和反復回数の時刻歴を示す図である。本発明の実施例における数値計算方法の解析例での同時緩和反復回数の時刻歴を示す図である。本発明の実施例における数値計算方法の解析例での要素数1350の解析メッシュと要素数5400の解析メッシュを示す図である。本発明の実施例における数値計算方法の解析例でのcavity内中心線に沿った速度分布を示す図である。本発明の実施例における数値計算方法の解析例での圧力分布を示す図である。本発明の実施例における数値計算方法の解析例での同時緩和反復回数の時刻歴を示す図である。

Claims

コンピュータを用いてGSMAC有限要素法により非圧縮性流体の流れを解析する数値計算方法おいて、任意の形状の要素を使用して解析対象をメッシュ化した構造データを用い、任意の次数の補間関数を使用して流体の速度と圧力を表わした計算式を用い、離散化Poisson方程式を優対角近似した計算式を用いて計算する計算手段によって、圧力のデータがある節点の衛星要素内の速度から該圧力のデータがある節点における速度の発散の平均値を求め、該平均値より修正速度ポテンシャルを求める計算式を用いて速度と圧力を同時緩和法で収束するまで計算して、収束した結果について時間進行計算手段により時間進行することを特徴とする数値計算方法。
速度の補間関数として２次関数を用い、圧力の補間関数として１次関数を用いることを特徴とする請求項１記載の数値計算方法。
前記要素として、４角形要素を使用し、速度の補間関数として２次関数を使用し、圧力の補間関数として１次関数を使用することを特徴とする請求項１記載の数値計算方法。
前記要素として、３角形要素を使用し、速度の補間関数として２次関数を使用し、圧力の補間関数として１次関数を使用することを特徴とする請求項１記載の数値計算方法。
前記要素として、４角形要素と３角形要素を組み合わせて使用し、速度の補間関数として２次関数を使用し、圧力の補間関数として１次関数を使用することを特徴とする請求項１記載の数値計算方法。
請求項１〜５のいずれかに記載の数値計算方法をコンピュータに実行させるための処理手順を記述したコンピュータプログラム。
請求項６記載のコンピュータプログラムを記録したコンピュータ読取可能な記録媒体。