JP2005050063A - 流体の数値計算方法及び装置及び制御プログラム及び記憶媒体 - Google Patents
流体の数値計算方法及び装置及び制御プログラム及び記憶媒体 Download PDFInfo
- Publication number
- JP2005050063A JP2005050063A JP2003205055A JP2003205055A JP2005050063A JP 2005050063 A JP2005050063 A JP 2005050063A JP 2003205055 A JP2003205055 A JP 2003205055A JP 2003205055 A JP2003205055 A JP 2003205055A JP 2005050063 A JP2005050063 A JP 2005050063A
- Authority
- JP
- Japan
- Prior art keywords
- difference
- fluid
- flux
- equation
- numerical calculation
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Withdrawn
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
- Complex Calculations (AREA)
Abstract
【課題】効率的で数値拡散が少なく、状態方程式の関数形によらない汎用的な数値計算手法として利用することができる流体の数値計算方法及び装置及び制御プログラム及び記憶媒体を提供する。
【解決手段】ＣＰＵ１１６により保存形式の方程式を計算する際に、流束をその物理的性質に応じて分離し、風上差分と中心差分型の蛙跳び差分法を併用して時間発展を行う。
【選択図】 図８
【解決手段】ＣＰＵ１１６により保存形式の方程式を計算する際に、流束をその物理的性質に応じて分離し、風上差分と中心差分型の蛙跳び差分法を併用して時間発展を行う。
【選択図】 図８
Description
【０００１】
【発明の属する技術分野】
本発明は、圧縮性を有する流体等の流体の数値計算方法及び装置及びその計算方法をコンピュータにて実行するための制御プログラム及びその制御プログラムを格納した記憶媒体に関する。
【０００２】
【従来の技術】
一般に、圧縮性を有する流体の数値解析手法は、大きくは、風上化を意識する解法と、格子系を工夫することにより風上化を意識しない解法とに分類することができる。
【０００３】
近年、効率的で安定性の高い圧縮性流体の風上化解法として、ＡＵＳＭ（Ａｄｖｅｃｔｉｏｎ Ｕｐｓｔｒｅａｍ Ｓｐｌｉｔｔｉｎｇ Ｍｅｔｈｏｄ）スキームが提案されている（非特許文献１）。このＡＵＳＭスキームにおいては、非線形性の原因となる圧力を別扱いし、マッハ数に応じた圧力評価を行うために、状態方程式の関数形に応じた音速の評価を行っている。
【０００４】
また、風上化を意識しない解法としては、中心差分スキーム（セントラルスキーム）（非特許文献２）や、ＣＥ／ＳＥスキームが提案されている（非特許文献３）。中心差分スキームにおいては双対格子系を、ＣＥ／ＳＥスキームにおいては時空間における保存則を考慮するために、複数の格子系を採用し計算を行う。
【０００５】
【非特許文献１】
Ｍ． Ｓ． Ｌｉｕ ａｎｄ Ｃ． Ｊ． Ｓｔｅｆｆｅｎ， Ｊｒ：Ａ Ｎｅｗ Ｆｌｕｘ Ｓｐｌｉｔｔｉｎｇ Ｓｃｈｅｍｅ， Ｊｏｕｒｎａｌ ｏｆ Ｃｏｍｐｕｔａｔｉｏｎａｌ Ｐｈｙｓｉｃｓ， １０７， ２３−３９（１９９３）
【非特許文献２】
Ｈ． Ｎｅｓｓｙａｈｕ， Ｅ． Ｔａｄｍｏｒ， “Ｎｏｎ−ｏｓｃｉｌｌａｔｏｒｙ ｃｅｎｔｒａｌ ｄｉｆｆｅｒｅｎｃｉｎｇ ｆｏｒ ｈｙｐｅｒｂｏｌｉｃ ｃｏｎｓｅｒｖａｔｉｏｎ ｌａｗｓ” Ｊｏｕｒｎａｌ ｏｆ Ｃｏｍｐｕｔａｔｉｏｎａｌ Ｐｈｙｓｉｃｓ， ８７， ４０８−４６３（１９９０）
【非特許文献３】
Ｓ． Ｃ． Ｃｈａｎｇ “Ｔｈｅ Ｍｅｔｈｏｄ ｏｆ Ｓｐａｃｅ−Ｔｉｍｅ Ｃｏｎｓｅｒｖａｔｉｏｎ Ｅｌｅｍｅｎｔ ａｎｄ Ｓｏｌｕｔｉｏｎ Ｅｌｅｍｅｎｔ −Ａ Ｎｅｗ Ａｐｐｒｏａｃｈ ｆｏｒ Ｓｏｌｖｉｎｇ ｔｈｅ Ｎａｖｉｅｒ−Ｓｔｏｋｅｓ ａｎｄ Ｅｕｌｅｒ Ｅｑｕａｔｉｏｎｓ” Ｊｏｕｒｎａｌ ｏｆ Ｃｏｍｐｕｔａｔｉｏｎａｌ Ｐｈｙｓｉｃｓ， １１９， ２９５−３２４（１９９５）
【０００６】
【発明が解決しようとする課題】
ＡＵＳＭスキームは、状態方程式の関数形に応じた音速の評価を行う必要があり、中心差分スキームやＣＥ／ＳＥスキームは、計算負荷が大きく、数値拡散が生じ易い複数の格子系を採用する必要がある。
【０００７】
本発明は、このような従来技術の有する事情に鑑みてなされたもので、その目的は、効率的で数値拡散が少なく、状態方程式の関数形によらない汎用的な数値計算手法として利用することができる流体の数値計算方法及び装置及び制御プログラム及び記憶媒体を提供することである。
【０００８】
【課題を解決するための手段】
上記目的を達成するために本発明の流体の数値計算方法は、保存形式の方程式を計算する際に、流束をその物理的性質に応じて分離し、風上差分と中心差分型の蛙跳び差分法を併用して時間発展を行うことを特徴とする。
【０００９】
また、上記目的を達成するために本発明の流体の数値計算方法は、保存形式の圧縮性流体方程式を計算する際に、流束を移流部分と圧力依存部分とに分離し、前記移流部分に関しては流速に基づく風上差分を用いて時間発展を行い、前記圧力依存部分については中心差分型の蛙跳び差分法を用いて時間発展を行うことを特徴とする。
【００１０】
また、上記目的を達成するために本発明の流体の数値計算装置は、保存形式の方程式を計算する際に、流束をその物理的性質に応じて分離し、風上差分と中心差分型の蛙跳び差分法を併用して時間発展を行うように制御する制御手段を有することを特徴とする。
【００１１】
また、上記目的を達成するために本発明の流体の数値計算装置は、保存形式の圧縮性流体方程式を計算する際に、流束を移流部分と圧力依存部分とに分離し、前記移流部分に関しては流速に基づく風上差分を用いて時間発展を行い、前記圧力依存部分については中心差分型の蛙跳び差分法を用いて時間発展を行うように制御する制御手段を有することを特徴とする。
【００１２】
また、上記目的を達成するために本発明の制御プログラムは、前記流体の数値計算方法をコンピュータにて実行するためのプログラムコードから成ることを特徴とする。
【００１３】
更に、上記目的を達成するために本発明の記憶媒体は、前記制御プログラムを格納したことを特徴とする。
【００１４】
【発明の実施の形態】
以下、本発明の各実施の形態を図面に基づき説明する。
【００１５】
（第１の実施形態）
まず、本発明の第１の実施形態を図１乃至図６に基づき説明する。
【００１６】
下記式▲１▼は、一次元の場合の圧縮性を有する流体の基本方程式であり、Ｕは保存量の密度・運動量密度・単位体積当たりの全エネルギー、Ｆは流束である。ここで流束とは、下記式▲３▼で表される単位時間・単位体積当たりの量である。この方程式▲１▼は、保存量の時間変化と流束の収支の和が保存することを表している。
【００１７】
また、下記式▲２▼は、保存量の密度・運動量密度・単位体積当たりの全エネルギーＵを求めるための計算式である。
【００１８】
更に、下記式▲３▼は、流束Ｆを求めるための計算式である。
【００１９】
図１は、本実施形態の特徴を示す図である。ここでは中心差分格子を採用しているが、流束を評価する際に、上記式▲３▼のように移流部分と圧力依存部分と散逸部分とに分割し、移流部分に関しては評価点における流速の方向に応じた風上差分を適用し、圧力依存部分に関しては中心差分によって計算した値を用いる。また、散逸部分については陽解法若しくは陰解法によって評価する。
【００２０】
図２は、一次元の場合の通常格子と中心差分格子とを示す図であり、同図においてｆｌｕｘはフラックス（流束）である。
【００２１】
図２の通常格子においてはセル中心とセル境界の位置は固定しているのに対し、中心差分格子においては双対格子系であり、セル中心とセル境界の位置が、１時間ステップ毎に半メッシュ移動するようになっている。また、通常格子においてはセル中心に保存量を定義し、セル境界において流束を評価する。また、中心差分格子においては保存量の定義位置と流束の評価位置とが同じである。
【００２２】
以下、本実施形態に係る圧縮性流体の数値計算方法による計算手順について、図３及び図４に基づき説明する。
【００２３】
図３及び図４は、本実施形態に係る圧縮性流体の数値計算方法による計算手順を示すフローチャートである。
【００２４】
まず、ステップＳ３０１において、保存量Ｕの初期設定を時刻ｎ及び時刻ｎ−１について行う（図４のステップＳ４０１）。
【００２５】
次に、ステップＳ３０２において、時刻ｎの内部エネルギー密度εを計算する（図４のステップＳ４０２）。
【００２６】
次に、ステップＳ３０３において、時刻ｎの圧力ｐを計算する（図４のステップＳ４０３）。
【００２７】
次に、ステップＳ３０４において、時刻ｎにおける流束Ｋと時刻ｎ＋１における流束の移流部分を計算する（図４のステップＳ４０４）。
【００２８】
移流部分の計算のための風上差分には、一次の風上差分だけでなく、高解像度スキーム（ＴＶＤ（Ｔｏｔａｌ Ｖａｒｉａｔｉｏｎ Ｄｉｍｉｎｉｓｈｉｎｇ）、ＥＮＯ（Ｅｓｓｅｎｔｉａｌｌｙ Ｎｏｎ Ｏｓｃｉｌｌａｔｏｒｙ）、ＰＰ（Ｐｏｓｉｔｉｖｅｌｙ Ｐｒｅｓｅｒｖｉｎｇ）スキーム等）を用いることができる。
【００２９】
次に、ステップＳ３０５（図４のステップＳ４０５）において、保存量Ｕの更新を行い、時刻ｎ＋１の値を求める。ここで、流束の圧力依存部分に関しては中心差分を、移流部分については風上差分を、散逸部分に関しては時刻ｎ−１（陽解法）の値か、時刻ｎ＋１（陰解法）の値のいずれかを用いる。
【００３０】
本手法は、流束の圧力依存部分に関しては中心差分を適用するため、ＡＵＳＭスキームのように圧力評価のための音速を計算する必要がなく、状態方程式によらない汎用的な計算手法として利用することができる。
【００３１】
また、計算の安定条件が厳しく、数値振動が生じるような場合には、人工粘性項等を付与させることにより、計算の安定条件を満たすことが可能である。
【００３２】
図８は、本実施形態に係る圧縮性流体の数値計算方法を実行する計算機（流体の数値計算装置）の構成を示すブロック図である。
【００３３】
図８に示すように、この計算機は、複数のＣＰＵ１１６やメモリ１１８等からなる演算処理部、ハードディスク１１７やフロッピー（登録商標）ディスク１１０からなる記憶装置部、キーボード１１５やマウス１１４等からなる入力部、そしてディスプレイ１１２等からなる出力表示部等で構成されているシステムである。
【００３４】
この計算機は、図３及び４に示したフローチャートに対応するプログラムをメモリ１１８上にロードし、必要な計算領域を確保し、適当な方法で入力された領域形状や格子点座標及び物理量に基づき所定の演算処理を行い、その結果得られた領域上での物理量の値をハードディスク１１７等に書き込んで記憶させたり、ディスプレイ１１２上に表示させる。
【００３５】
図５は、圧縮性流体の代表的なベンチマークテストである衝撃波管問題を、本手法と中心差分スキームの一つであるＣＷＥＮＯスキームで計算した結果の一例を示す図である。
【００３６】
図６は、図５の一部を拡大した図である。
【００３７】
図５及び図６において、実線は解析解を、黒丸実線は本手法の計算結果を、三角点線はＣＷＥＮＯスキームの計算結果を示している。
【００３８】
図５及び図６にて明確なように、本手法の方が数値拡散が少なく、高い解像度が得られていることが分る。
【００３９】
以上詳述したように本実施の形態によれば、風上化解法の利点と非風上化解法の利点を採用し、効率的で数値拡散が少なく、状態方程式の関数形によらない汎用的な数値計算手法として利用することができる。
【００４０】
また、本実施の形態によれば、計算スキームの変更無しに、理想気体の状態方程式のみならず、ｖａｎ ｄｅｒ Ｗａａｌｓ状態方程式等を扱うことができる。
【００４１】
（第２の実施形態）
次に、本発明の第２の実施形態を図７に基づき説明する。
【００４２】
下記式▲４▼は、二次元の場合の圧縮性流体の基本方程式であり、Ｕは保存量の密度・ｘ方向運動量密度・ｙ方向運動量密度・単位体積当たりの全エネルギー、Ｆはｘ方向流束、Ｇはｙ方向流束である。ここでｘ方向流束とは、下記式で表される単位時間・単位体積当たりの量である。この方程式▲４▼は、保存量の時間変化と流束の収支の和が保存することを表している。
【００４３】
また、下記式▲５▼は、保存量の密度・運動量密度・単位体積当たりの全エネルギーＵを求めるための計算式である。
【００４４】
また、下記式▲６▼は、ｘ方向流束Ｆを求めるための計算式である。
【００４５】
更に、下記式▲７▼は、ｙ方向流束Ｇを求めるための計算式である。
【００４６】
図７は、二次元の場合の通常格子と中心差分格子の一例を示す図であり、同図（ａ）は通常格子を、同図（ｂ）は中心差分格子を示す。
【００４７】
この二次元は、上述した第１の実施形態における一次元の場合と同様に、通常格子においてはセル中心とセル境界の位置は固定しているのに対し、中心差分格子ではセル中心とセル境界の位置が、１時間ステップ毎に半メッシュ（ｘ方向だけでなく、ｙ方向にも）移動するようになっている。
【００４８】
また、この二次元は、一次元の場合と同様に、通常格子においてはセル中心に保存量を定義し、セル境界において流束を評価するが、中心差分格子においては保存量の定義位置と流束の評価位置とが同じである。
【００４９】
本発明を二次元の場合に適用するには、流束を評価する際に、上記式▲６▼、▲７▼のようにｘ方向流束とｙ方向流束とを移流部分と圧力依存部分と散逸部分とに分割し、移流部分に関しては評価点における流速の方向に応じた風上差分を適用し、圧力依存部分に関しては中心差分によって計算した値を用いれば良い。また、散逸部分については陽解法若しくは陰解法によって評価する。
【００５０】
これらの手続きは、ｙ方向の成分が付け加わったことを除けば、上述した第１の実施形態における一次元の場合と同様である。また、流体方程式系に、重力等の外力が付け加わった場合や、三次元の場合にも本手法は利用できる。
【００５１】
【発明の効果】
以上説明したように本発明によれば、効率的で数値拡散が少なく、状態方程式の関数形によらない汎用的な数値計算手法として利用することができるという効果を奏する。
【図面の簡単な説明】
【図１】本発明の第１の実施形態に係る流体の数値計算方法を示す概念図である。
【図２】本発明の第１の実施形態に係る流体の数値計算方法における通常格子とセントラル格子の概念図である。
【図３】本発明の第１の実施形態に係る流体の数値計算方法における計算手順を示すフローチャートである。
【図４】本発明の第１の実施形態に係る流体の数値計算方法における計算手順を示すフローチャートである。
【図５】本発明の第１の実施形態に係る流体の数値計算方法による計算結果と中心差分スキームによる計算結果の比較例を示す図である。
【図６】図５の一部を拡大した図である。
【図７】本発明の第２の実施形態に係る流体の数値計算方法における二次元の通常格子と中心差分格子の概念図である。
【図８】本実施形態に係る圧縮性流体の数値計算方法を実行する計算機（流体の数値計算装置）の構成を示すブロック図である。
【符号の説明】
１１０ フロッピー（登録商標）ディスク
１１２ ディスプレイ
１１４ マウス
１１５ キーボード
１１６ ＣＰＵ
１１７ ハードディスク
１１８ メモリ
【発明の属する技術分野】
本発明は、圧縮性を有する流体等の流体の数値計算方法及び装置及びその計算方法をコンピュータにて実行するための制御プログラム及びその制御プログラムを格納した記憶媒体に関する。
【０００２】
【従来の技術】
一般に、圧縮性を有する流体の数値解析手法は、大きくは、風上化を意識する解法と、格子系を工夫することにより風上化を意識しない解法とに分類することができる。
【０００３】
近年、効率的で安定性の高い圧縮性流体の風上化解法として、ＡＵＳＭ（Ａｄｖｅｃｔｉｏｎ Ｕｐｓｔｒｅａｍ Ｓｐｌｉｔｔｉｎｇ Ｍｅｔｈｏｄ）スキームが提案されている（非特許文献１）。このＡＵＳＭスキームにおいては、非線形性の原因となる圧力を別扱いし、マッハ数に応じた圧力評価を行うために、状態方程式の関数形に応じた音速の評価を行っている。
【０００４】
また、風上化を意識しない解法としては、中心差分スキーム（セントラルスキーム）（非特許文献２）や、ＣＥ／ＳＥスキームが提案されている（非特許文献３）。中心差分スキームにおいては双対格子系を、ＣＥ／ＳＥスキームにおいては時空間における保存則を考慮するために、複数の格子系を採用し計算を行う。
【０００５】
【非特許文献１】
Ｍ． Ｓ． Ｌｉｕ ａｎｄ Ｃ． Ｊ． Ｓｔｅｆｆｅｎ， Ｊｒ：Ａ Ｎｅｗ Ｆｌｕｘ Ｓｐｌｉｔｔｉｎｇ Ｓｃｈｅｍｅ， Ｊｏｕｒｎａｌ ｏｆ Ｃｏｍｐｕｔａｔｉｏｎａｌ Ｐｈｙｓｉｃｓ， １０７， ２３−３９（１９９３）
【非特許文献２】
Ｈ． Ｎｅｓｓｙａｈｕ， Ｅ． Ｔａｄｍｏｒ， “Ｎｏｎ−ｏｓｃｉｌｌａｔｏｒｙ ｃｅｎｔｒａｌ ｄｉｆｆｅｒｅｎｃｉｎｇ ｆｏｒ ｈｙｐｅｒｂｏｌｉｃ ｃｏｎｓｅｒｖａｔｉｏｎ ｌａｗｓ” Ｊｏｕｒｎａｌ ｏｆ Ｃｏｍｐｕｔａｔｉｏｎａｌ Ｐｈｙｓｉｃｓ， ８７， ４０８−４６３（１９９０）
【非特許文献３】
Ｓ． Ｃ． Ｃｈａｎｇ “Ｔｈｅ Ｍｅｔｈｏｄ ｏｆ Ｓｐａｃｅ−Ｔｉｍｅ Ｃｏｎｓｅｒｖａｔｉｏｎ Ｅｌｅｍｅｎｔ ａｎｄ Ｓｏｌｕｔｉｏｎ Ｅｌｅｍｅｎｔ −Ａ Ｎｅｗ Ａｐｐｒｏａｃｈ ｆｏｒ Ｓｏｌｖｉｎｇ ｔｈｅ Ｎａｖｉｅｒ−Ｓｔｏｋｅｓ ａｎｄ Ｅｕｌｅｒ Ｅｑｕａｔｉｏｎｓ” Ｊｏｕｒｎａｌ ｏｆ Ｃｏｍｐｕｔａｔｉｏｎａｌ Ｐｈｙｓｉｃｓ， １１９， ２９５−３２４（１９９５）
【０００６】
【発明が解決しようとする課題】
ＡＵＳＭスキームは、状態方程式の関数形に応じた音速の評価を行う必要があり、中心差分スキームやＣＥ／ＳＥスキームは、計算負荷が大きく、数値拡散が生じ易い複数の格子系を採用する必要がある。
【０００７】
本発明は、このような従来技術の有する事情に鑑みてなされたもので、その目的は、効率的で数値拡散が少なく、状態方程式の関数形によらない汎用的な数値計算手法として利用することができる流体の数値計算方法及び装置及び制御プログラム及び記憶媒体を提供することである。
【０００８】
【課題を解決するための手段】
上記目的を達成するために本発明の流体の数値計算方法は、保存形式の方程式を計算する際に、流束をその物理的性質に応じて分離し、風上差分と中心差分型の蛙跳び差分法を併用して時間発展を行うことを特徴とする。
【０００９】
また、上記目的を達成するために本発明の流体の数値計算方法は、保存形式の圧縮性流体方程式を計算する際に、流束を移流部分と圧力依存部分とに分離し、前記移流部分に関しては流速に基づく風上差分を用いて時間発展を行い、前記圧力依存部分については中心差分型の蛙跳び差分法を用いて時間発展を行うことを特徴とする。
【００１０】
また、上記目的を達成するために本発明の流体の数値計算装置は、保存形式の方程式を計算する際に、流束をその物理的性質に応じて分離し、風上差分と中心差分型の蛙跳び差分法を併用して時間発展を行うように制御する制御手段を有することを特徴とする。
【００１１】
また、上記目的を達成するために本発明の流体の数値計算装置は、保存形式の圧縮性流体方程式を計算する際に、流束を移流部分と圧力依存部分とに分離し、前記移流部分に関しては流速に基づく風上差分を用いて時間発展を行い、前記圧力依存部分については中心差分型の蛙跳び差分法を用いて時間発展を行うように制御する制御手段を有することを特徴とする。
【００１２】
また、上記目的を達成するために本発明の制御プログラムは、前記流体の数値計算方法をコンピュータにて実行するためのプログラムコードから成ることを特徴とする。
【００１３】
更に、上記目的を達成するために本発明の記憶媒体は、前記制御プログラムを格納したことを特徴とする。
【００１４】
【発明の実施の形態】
以下、本発明の各実施の形態を図面に基づき説明する。
【００１５】
（第１の実施形態）
まず、本発明の第１の実施形態を図１乃至図６に基づき説明する。
【００１６】
下記式▲１▼は、一次元の場合の圧縮性を有する流体の基本方程式であり、Ｕは保存量の密度・運動量密度・単位体積当たりの全エネルギー、Ｆは流束である。ここで流束とは、下記式▲３▼で表される単位時間・単位体積当たりの量である。この方程式▲１▼は、保存量の時間変化と流束の収支の和が保存することを表している。
【００１７】
また、下記式▲２▼は、保存量の密度・運動量密度・単位体積当たりの全エネルギーＵを求めるための計算式である。
【００１８】
更に、下記式▲３▼は、流束Ｆを求めるための計算式である。
図１は、本実施形態の特徴を示す図である。ここでは中心差分格子を採用しているが、流束を評価する際に、上記式▲３▼のように移流部分と圧力依存部分と散逸部分とに分割し、移流部分に関しては評価点における流速の方向に応じた風上差分を適用し、圧力依存部分に関しては中心差分によって計算した値を用いる。また、散逸部分については陽解法若しくは陰解法によって評価する。
【００２０】
図２は、一次元の場合の通常格子と中心差分格子とを示す図であり、同図においてｆｌｕｘはフラックス（流束）である。
【００２１】
図２の通常格子においてはセル中心とセル境界の位置は固定しているのに対し、中心差分格子においては双対格子系であり、セル中心とセル境界の位置が、１時間ステップ毎に半メッシュ移動するようになっている。また、通常格子においてはセル中心に保存量を定義し、セル境界において流束を評価する。また、中心差分格子においては保存量の定義位置と流束の評価位置とが同じである。
【００２２】
以下、本実施形態に係る圧縮性流体の数値計算方法による計算手順について、図３及び図４に基づき説明する。
【００２３】
図３及び図４は、本実施形態に係る圧縮性流体の数値計算方法による計算手順を示すフローチャートである。
【００２４】
まず、ステップＳ３０１において、保存量Ｕの初期設定を時刻ｎ及び時刻ｎ−１について行う（図４のステップＳ４０１）。
【００２５】
次に、ステップＳ３０２において、時刻ｎの内部エネルギー密度εを計算する（図４のステップＳ４０２）。
【００２６】
次に、ステップＳ３０３において、時刻ｎの圧力ｐを計算する（図４のステップＳ４０３）。
【００２７】
次に、ステップＳ３０４において、時刻ｎにおける流束Ｋと時刻ｎ＋１における流束の移流部分を計算する（図４のステップＳ４０４）。
【００２８】
移流部分の計算のための風上差分には、一次の風上差分だけでなく、高解像度スキーム（ＴＶＤ（Ｔｏｔａｌ Ｖａｒｉａｔｉｏｎ Ｄｉｍｉｎｉｓｈｉｎｇ）、ＥＮＯ（Ｅｓｓｅｎｔｉａｌｌｙ Ｎｏｎ Ｏｓｃｉｌｌａｔｏｒｙ）、ＰＰ（Ｐｏｓｉｔｉｖｅｌｙ Ｐｒｅｓｅｒｖｉｎｇ）スキーム等）を用いることができる。
【００２９】
次に、ステップＳ３０５（図４のステップＳ４０５）において、保存量Ｕの更新を行い、時刻ｎ＋１の値を求める。ここで、流束の圧力依存部分に関しては中心差分を、移流部分については風上差分を、散逸部分に関しては時刻ｎ−１（陽解法）の値か、時刻ｎ＋１（陰解法）の値のいずれかを用いる。
【００３０】
本手法は、流束の圧力依存部分に関しては中心差分を適用するため、ＡＵＳＭスキームのように圧力評価のための音速を計算する必要がなく、状態方程式によらない汎用的な計算手法として利用することができる。
【００３１】
また、計算の安定条件が厳しく、数値振動が生じるような場合には、人工粘性項等を付与させることにより、計算の安定条件を満たすことが可能である。
【００３２】
図８は、本実施形態に係る圧縮性流体の数値計算方法を実行する計算機（流体の数値計算装置）の構成を示すブロック図である。
【００３３】
図８に示すように、この計算機は、複数のＣＰＵ１１６やメモリ１１８等からなる演算処理部、ハードディスク１１７やフロッピー（登録商標）ディスク１１０からなる記憶装置部、キーボード１１５やマウス１１４等からなる入力部、そしてディスプレイ１１２等からなる出力表示部等で構成されているシステムである。
【００３４】
この計算機は、図３及び４に示したフローチャートに対応するプログラムをメモリ１１８上にロードし、必要な計算領域を確保し、適当な方法で入力された領域形状や格子点座標及び物理量に基づき所定の演算処理を行い、その結果得られた領域上での物理量の値をハードディスク１１７等に書き込んで記憶させたり、ディスプレイ１１２上に表示させる。
【００３５】
図５は、圧縮性流体の代表的なベンチマークテストである衝撃波管問題を、本手法と中心差分スキームの一つであるＣＷＥＮＯスキームで計算した結果の一例を示す図である。
【００３６】
図６は、図５の一部を拡大した図である。
【００３７】
図５及び図６において、実線は解析解を、黒丸実線は本手法の計算結果を、三角点線はＣＷＥＮＯスキームの計算結果を示している。
【００３８】
図５及び図６にて明確なように、本手法の方が数値拡散が少なく、高い解像度が得られていることが分る。
【００３９】
以上詳述したように本実施の形態によれば、風上化解法の利点と非風上化解法の利点を採用し、効率的で数値拡散が少なく、状態方程式の関数形によらない汎用的な数値計算手法として利用することができる。
【００４０】
また、本実施の形態によれば、計算スキームの変更無しに、理想気体の状態方程式のみならず、ｖａｎ ｄｅｒ Ｗａａｌｓ状態方程式等を扱うことができる。
【００４１】
（第２の実施形態）
次に、本発明の第２の実施形態を図７に基づき説明する。
【００４２】
下記式▲４▼は、二次元の場合の圧縮性流体の基本方程式であり、Ｕは保存量の密度・ｘ方向運動量密度・ｙ方向運動量密度・単位体積当たりの全エネルギー、Ｆはｘ方向流束、Ｇはｙ方向流束である。ここでｘ方向流束とは、下記式で表される単位時間・単位体積当たりの量である。この方程式▲４▼は、保存量の時間変化と流束の収支の和が保存することを表している。
【００４３】
また、下記式▲５▼は、保存量の密度・運動量密度・単位体積当たりの全エネルギーＵを求めるための計算式である。
【００４４】
また、下記式▲６▼は、ｘ方向流束Ｆを求めるための計算式である。
【００４５】
更に、下記式▲７▼は、ｙ方向流束Ｇを求めるための計算式である。
図７は、二次元の場合の通常格子と中心差分格子の一例を示す図であり、同図（ａ）は通常格子を、同図（ｂ）は中心差分格子を示す。
【００４７】
この二次元は、上述した第１の実施形態における一次元の場合と同様に、通常格子においてはセル中心とセル境界の位置は固定しているのに対し、中心差分格子ではセル中心とセル境界の位置が、１時間ステップ毎に半メッシュ（ｘ方向だけでなく、ｙ方向にも）移動するようになっている。
【００４８】
また、この二次元は、一次元の場合と同様に、通常格子においてはセル中心に保存量を定義し、セル境界において流束を評価するが、中心差分格子においては保存量の定義位置と流束の評価位置とが同じである。
【００４９】
本発明を二次元の場合に適用するには、流束を評価する際に、上記式▲６▼、▲７▼のようにｘ方向流束とｙ方向流束とを移流部分と圧力依存部分と散逸部分とに分割し、移流部分に関しては評価点における流速の方向に応じた風上差分を適用し、圧力依存部分に関しては中心差分によって計算した値を用いれば良い。また、散逸部分については陽解法若しくは陰解法によって評価する。
【００５０】
これらの手続きは、ｙ方向の成分が付け加わったことを除けば、上述した第１の実施形態における一次元の場合と同様である。また、流体方程式系に、重力等の外力が付け加わった場合や、三次元の場合にも本手法は利用できる。
【００５１】
【発明の効果】
以上説明したように本発明によれば、効率的で数値拡散が少なく、状態方程式の関数形によらない汎用的な数値計算手法として利用することができるという効果を奏する。
【図面の簡単な説明】
【図１】本発明の第１の実施形態に係る流体の数値計算方法を示す概念図である。
【図２】本発明の第１の実施形態に係る流体の数値計算方法における通常格子とセントラル格子の概念図である。
【図３】本発明の第１の実施形態に係る流体の数値計算方法における計算手順を示すフローチャートである。
【図４】本発明の第１の実施形態に係る流体の数値計算方法における計算手順を示すフローチャートである。
【図５】本発明の第１の実施形態に係る流体の数値計算方法による計算結果と中心差分スキームによる計算結果の比較例を示す図である。
【図６】図５の一部を拡大した図である。
【図７】本発明の第２の実施形態に係る流体の数値計算方法における二次元の通常格子と中心差分格子の概念図である。
【図８】本実施形態に係る圧縮性流体の数値計算方法を実行する計算機（流体の数値計算装置）の構成を示すブロック図である。
【符号の説明】
１１０ フロッピー（登録商標）ディスク
１１２ ディスプレイ
１１４ マウス
１１５ キーボード
１１６ ＣＰＵ
１１７ ハードディスク
１１８ メモリ
Claims (6)
- 保存形式の方程式を計算する際に、流束をその物理的性質に応じて分離し、風上差分と中心差分型の蛙跳び差分法を併用して時間発展を行うことを特徴とする流体の数値計算方法。
- 保存形式の流体方程式を計算する際に、流束を移流部分と圧力依存部分とに分離し、前記移流部分に関しては流速に基づく風上差分を用いて時間発展を行い、前記圧力依存部分については中心差分型の蛙跳び差分法を用いて時間発展を行うことを特徴とする流体の数値計算方法。
- 保存形式の方程式を計算する際に、流束をその物理的性質に応じて分離し、風上差分と中心差分型の蛙跳び差分法を併用して時間発展を行うように制御する制御手段を有することを特徴とする流体の数値計算装置。
- 保存形式の流体方程式を計算する際に、流束を移流部分と圧力依存部分とに分離し、前記移流部分に関しては流速に基づく風上差分を用いて時間発展を行い、前記圧力依存部分については中心差分型の蛙跳び差分法を用いて時間発展を行うように制御する制御手段を有することを特徴とする流体の数値計算装置。
- 請求項１または２記載の流体の数値計算方法をコンピュータにて実行するためのプログラムコードから成ることを特徴とするコンピュータ読み取り可能な制御プログラム。
- 請求項５記載の制御プログラムを格納したことを特徴とする記憶媒体。
Priority Applications (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2003205055A JP2005050063A (ja) | 2003-07-31 | 2003-07-31 | 流体の数値計算方法及び装置及び制御プログラム及び記憶媒体 |
| US10/887,849 US7516054B2 (en) | 2003-07-31 | 2004-07-12 | Method and apparatus for numerical calculation related to fluids, program for implementing the method, and storage medium storing the program |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2003205055A JP2005050063A (ja) | 2003-07-31 | 2003-07-31 | 流体の数値計算方法及び装置及び制御プログラム及び記憶媒体 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| JP2005050063A true JP2005050063A (ja) | 2005-02-24 |
Family
ID=34100692
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2003205055A Withdrawn JP2005050063A (ja) | 2003-07-31 | 2003-07-31 | 流体の数値計算方法及び装置及び制御プログラム及び記憶媒体 |
Country Status (2)
| Country | Link |
|---|---|
| US (1) | US7516054B2 (ja) |
| JP (1) | JP2005050063A (ja) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102708233A (zh) * | 2012-04-24 | 2012-10-03 | 北京航空航天大学 | 一种适用于表面安装工艺的可靠性实现能力评估方法 |
| CN111274682A (zh) * | 2020-01-15 | 2020-06-12 | 桂林电子科技大学 | 基于蛙跳算法的数字微流控芯片在线测试路径优化方法 |
Families Citing this family (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103065039B (zh) * | 2012-12-10 | 2016-02-24 | 北京航空航天大学 | 一种基于欧拉公式的高精度正弦/余弦函数计算方法 |
| US11493066B2 (en) | 2016-01-20 | 2022-11-08 | Soliton Holdings | Generalized jet-effect and enhanced devices |
| US11705780B2 (en) | 2016-01-20 | 2023-07-18 | Soliton Holdings Corporation, Delaware Corporation | Generalized jet-effect and generalized generator |
| US20180266394A1 (en) * | 2016-01-20 | 2018-09-20 | Soliton Holdings Corporation, Delaware Corporation | Generalized Jet-Effect and Generalized Generator |
| US11499525B2 (en) | 2016-01-20 | 2022-11-15 | Soliton Holdings Corporation, Delaware Corporation | Generalized jet-effect and fluid-repellent corpus |
| GB201601094D0 (en) * | 2016-01-20 | 2016-03-02 | Soliton Holdings Corp Delaware Corp | Method for computational fluid dynamics and apparatuses for jet-effect use |
| US20180286972A1 (en) * | 2017-02-13 | 2018-10-04 | Ahmad Tarakji | Aluminum-rich field-plated nitride transistors for record high currents |
| CN110457806A (zh) * | 2019-08-02 | 2019-11-15 | 南京航空航天大学 | 基于交错网格的中心五阶weno格式的全流场模拟方法 |
| CN113222410B (zh) * | 2021-05-12 | 2023-05-19 | 西安建筑科技大学 | 双向式布局模式下货位分配模型建立方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP3034433B2 (ja) * | 1994-10-31 | 2000-04-17 | 核燃料サイクル開発機構 | 構造物の熱設計方法およびその設計に最適な数値計算装置 |
-
2003
- 2003-07-31 JP JP2003205055A patent/JP2005050063A/ja not_active Withdrawn
-
2004
- 2004-07-12 US US10/887,849 patent/US7516054B2/en not_active Expired - Fee Related
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102708233A (zh) * | 2012-04-24 | 2012-10-03 | 北京航空航天大学 | 一种适用于表面安装工艺的可靠性实现能力评估方法 |
| CN111274682A (zh) * | 2020-01-15 | 2020-06-12 | 桂林电子科技大学 | 基于蛙跳算法的数字微流控芯片在线测试路径优化方法 |
| CN111274682B (zh) * | 2020-01-15 | 2024-01-05 | 桂林电子科技大学 | 基于蛙跳算法的数字微流控芯片在线测试路径优化方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| US7516054B2 (en) | 2009-04-07 |
| US20050027498A1 (en) | 2005-02-03 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Hötzer et al. | The parallel multi-physics phase-field framework Pace3D | |
| Chen et al. | A local adaptive sampling method for reliability-based design optimization using Kriging model | |
| Kwon et al. | Skew-resistant parallel processing of feature-extracting scientific user-defined functions | |
| Ipek et al. | An approach to performance prediction for parallel applications | |
| Almgren et al. | CASTRO: A new compressible astrophysical solver. I. Hydrodynamics and self-gravity | |
| Mirin et al. | Very high resolution simulation of compressible turbulence on the IBM-SP system | |
| Lively et al. | Energy and performance characteristics of different parallel implementations of scientific applications on multicore systems | |
| Radtke et al. | On efficient simulations of multiscale kinetic transport | |
| Huang et al. | Reliability analysis of structures using stochastic response surface method and saddlepoint approximation | |
| JP2005050063A (ja) | 流体の数値計算方法及び装置及び制御プログラム及び記憶媒体 | |
| Pal et al. | An overview of the numerical approaches to water hammer modelling: The ongoing quest for practical and accurate numerical approaches | |
| Huang et al. | Hybrid multiple-relaxation-time lattice-Boltzmann finite-difference method for axisymmetric multiphase flows | |
| JP2014013424A (ja) | 気液二相流シミュレーションの計算方法、プログラム、および装置 | |
| CN101887474A (zh) | 基于有限元法与广义傅里叶级数法的结构振动分析方法 | |
| Lu et al. | Data‐driven two‐level performance evaluation of eddy‐current tuned mass damper for building structures using shaking table and field testing | |
| Chen et al. | Towards high-accuracy deep learning inference of compressible turbulent flows over aerofoils | |
| Song et al. | Dynamic reliability analysis of the stability of the 3D slope of a rockfill dam based on the copula function and generalized probability density evolution method | |
| Xu et al. | Modeling techniques for strain-range-dependent hardening behavior of low-yield-point steel shear panel dampers | |
| Beck et al. | Application and development of the high order discontinuous Galerkin spectral element method for compressible multiscale flows | |
| Zhang et al. | Extended objective-oriented sequential sampling method for robust design of complex systems against design uncertainty | |
| Hadimlioglu et al. | Visualization of flooding using adaptive spatial resolution | |
| Soroushian et al. | Simultaneous Reduction in Time and Space for Dynamic Finite Element Models of Beam–Column Assemblies | |
| Ding et al. | Parameter identification of airfoil systems using an elite-based clustering Jaya algorithm and incremental vibration responses | |
| Jirasek et al. | Study of Benchmark Super Critical Wing at Mach 0.8 | |
| Ying et al. | Numerical simulation of aerostatic force coefficients of bridge decks using continuous torsional motion technique |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| RD03 | Notification of appointment of power of attorney |
Free format text: JAPANESE INTERMEDIATE CODE: A7423 Effective date: 20060417 |
|
| A300 | Withdrawal of application because of no request for examination |
Free format text: JAPANESE INTERMEDIATE CODE: A300 Effective date: 20061003 |