JPH01307826A - プログラム生成方法 - Google Patents

Abstract

(57)【要約】本公報は電子出願前の出願データであるた
め要約のデータは記録されません。

Description

【発明の詳細な説明】 〔産業上の利用分野〕 本発明は、プログラム生成方法に関し、特に数値計算プ
ログラムを作成する前の段階で、偏微分方程式と境界条
件で表わされる物理問題の数値シミュレーションを行う
際に、その偏微分方程式で表わされた物理法則の意味す
る保存性を維持させることにより、高精度な数値計算プ
ログラムを自動的に生成できるようにしたプログラム生
成方法に関する。

〔従来の技術〕

従来、差分法を実行するための離散式をコード化してプ
ログラムを作成するプログラム自動生成方法についでは
、例えば、情報処理学会「コンピュータシステム」シン
ポジウム、佐用他「数値シミュレーション語ＤＥＱＳＯ
Ｌｊに記載されている。これによると、内点では偏微分
方程式を、また境界では境界条件式を、それぞれテーラ
−展開により求めた離散化公式に従って離散式を導いて
、それらの離散式をコード化することによりプログラム
を生成していた。

ここで、離散化とは、解析領域を複数の離散点と呼ばれ
る代表点により表わし、各微分演算項を離散点での値か
らなる近似式で表わすことである。

この場合、偏微分方程式および境界条件式は、各離散点
における物理量どうしの関係式となるので、これを離散
式と呼んでいる。また、その離散化の公式を離散化規則
と呼ぶ。

第３図および第４図は、それぞれ従来の方法による離散
化規則を示す図である。

３次元空間を想定した場合、微分演算子は５次の１０個
、すなわち、ａ　／　ａ　ｘ　、　ａ　／　ａ　ｙ　、
　ａ　／ａ　ｚ　、ａ２／　ａ　ｘ２．ａ２／　ａ　ｙ
２．ａ２／　ａ　ｚ２．　ｄｉｖ。

ｇｒａｄ、　ｒｏｔ、　ｖ”である。これらのうち、基
本となるものは、ａ　／　ａ　ｘ　、　ａ　／　ａ　ｙ
　、　ａ　／　ａ　ｚ　（７）３　ツであって、他の７
つは数学的にこれら３つの演算の組み合わせで示される
。第３図では、これらの基本となる３つの微分演算子ａ
　／　ａ　ｘ　＋　ａ　／　ａ　ｙ　。

ａ　／　ａ　Ｚに対する内点の離散式と、境界（下限）
離散式と、境界（上限）の離散式とが示されている。

また、第４図では、他の７つの微分演算子の定義、つま
り３つの基本演算子を用いた式が示されている。

従来の方法では、この数学的な規則に従って、ａ’／ａ
ｘ２〜ｖ２の式を展開していた。また、基本の演算子ａ
　／　ａ　ｘ　、　ａ　／　ａ　ｙ　、　ａ　／　ａ　
ｚに対して、従来の方法では、第３図に示すように、テ
ーラ−展開により求めた離散式を使用していた。すなわ
ち、各々、内点では中心差分を、境界では片側差分を、
それぞれ用いている。一般に、片側差分よりも中心差分
の方が近似精度は高いが、境界において中心差分を用い
ると離散参照点が領域外にはみ出るため、上記仕様とな
っていた。

なお、本例では、３次元のものを示しているが、１次元
、２次元の場合も同じように考える。

〔発明が解決しようとする課題〕

前述の従来方法では、次のような２つの課題が生じる。

その第１は、従来の方式では、偏微分方程式の物理的な
意味をとらえずに離散化を行っている点である。基本的
に偏微分方程式は、その物理現象の保存性を満たしてお
り、本来はその離散式も同じでなくではならない。しか
し、従来の方法では、微分演算子を全て数学的にとらえ
て、第３図に示すａ　／　ａ　ｘ　、　ａ　／　ａ　ｙ
　、　ａ　／　ａ　ｚ　（７）合成トシテ離散化してい
るので、物理現象の保存性を近似した離散式とはならな
い。すなわち、領域の境界部分では、正確な物理現象を
表わしていないことになる。微分方程式には、物理的に
その微分演算の領域的な範囲（これを、本発明ではコン
トロールボリュームと呼ぶ）の演算を意味するものと勾
配の演算を意味するものの２種類がある。そして、各々
の意味を考えて、適切な離散化を行う必要がある。

他の１つの課題は、境界上の点における離散式の生成で
ある。第３図に示したように、境界では、片側差分によ
る離散化規則を用いて、境界条件をｍｈ化している。例
えば、ａφ／　ａ　Ｘ　＝φが点ｉの境界条件であると
すると、点ｉでは（φ、−φ□−０）／△Ｘ＝φ工が離
散式となる。ところが、この離散式は、点ｉ＋（１／２
）でａφ／　ａ　Ｘ　＝φを内点と同じ中心差分により
離散化した式と同一である。中心差分は片側差分より近
似精度が高く、従って境界条件を片側差分することは、
本来の成立式が半メツシュ内側にずれることになる。つ
まり、これもまた偏微分方程式の持つ保存則を近似した
離散式とはならない。

本発明の目的は、これらの課題を解決し、偏微分方程式
が持つ保存則を満たした離散式の生成が可能であり、高
精度な数値計算用のシミュレーションプログラムを作成
することができるプログラム生成方法を提供することに
ある。

〔課題を解決するための手段〕

上記目的を達成するため、本発明のプログラム生成方法
は、物理現象を表わす解析領域の内点では偏微分方程式
を、該解析領域の境界では偏微分方程式及び境界条件式
を、それぞれ成立する微少領域においてテーラ−展開に
より求めた離散化公式に従って離散式を導くことにより
、該離散式をコード化してシミュレーションプログラム
を生成するプログラム生成方法において、上記物理現象
を表わす偏微分方程式と該物理現象を解析する解析領域
と該解析領域の境界条件と、該解析領域の離散化に必要
な離散化規則と、該解析領域の代表点である離散点の座
標値とを入力し、先ず該解析領域を上記座標値の各離散
点を含む互いに重複しない微小領域に分割し、上記偏微
分方程式の中から上記微小領域の演算を意味する微分演
算子を抽出し、該微分演算子を上記離散化規則に従って
離散化することにより中間コードである展開式を得、次
に上記解析領域の境界部分に対してのみは特別な演算を
行わせるため、該展開式の中の境界勾配の演算を意味す
る項を抽出し、これを上記境界条件で置換し、次に上記
展開式中の勾配演算を意味する微分演算子を上記離散化
規則に従って離散化することにより離散式を生成し、最
後に該離散式をコード化してプログラムを自動生成する
ことに特徴がある。

〔作　　用〕

本発明においでは、偏微分方程式が持つ保存性を精度よ
く近似した離散式を生成するために、解析領域を離散点
の各座標値を中心とした微小領域（これを、本発明では
コントロールボリュームと呼ぶ）に分割し、この微小領
域を用いて処理する方式（以下、この方式をコントロー
ルボリューム方式と呼ぶ）によって自動離散化を実現し
ている。

数式処理において、コントロールボリューム方式を用い
た離散化を自動化するためには、次の３つの動作原理を
用いる。

ｉ）第１の原理は、偏微分方程式にある微分演算子の物
理的な意味の把握、つまり微分演算子がコントロールボ
リュームの演算を意味するものか、あるいは勾配の演算
を意味するものかを自動的に把握することである。

ｉｉ）第２の原理は、第１の原理により把握した２種類
の微分演算に対して、各々物理的意味のある離散化規則
を適用し、自動的に離散化することである。

これら第１と第２の原理により、課題の１つである偏微
分方程式の物理的な意味をとらえた離散化が実現できる
。

１ｉｉ）第３の原理は、境界上の点の離散化において、
偏微分方程式を離散化する際に、その境界の勾配の演算
を意味する項に対し境界条件を自動的に取り入れること
である。

これによって、もう１つの課題である境界上における保
存性を満たす離散式の生成が可能となる。

以下、上記第１．第２および第３の各原理に従った作用
を述べる。

第１の原理、つまり偏微分方程式中の微分演算子の物理
的意味を把握する方法は、次のような動作を行う。例え
ば、第６図に示す偏微分方程式。

ｄｉｖ（ｇｒａｄ（φ））＝α（式５０３）についでは
、　ｄｉｖをコントロールボリュームの演算を意味する
微分演算子、ｇｒａｄを勾配の演算を意味する微分演算
子として把握する。他の場合にも各々の微分演算子につ
いて、式の形式および式中に現われた位置によりいずれ
の意味を示すのかを判断する。

第２の原理、つまり第１の原理で判断した各微分演算子
に対して、各々物理的に意味のある離散化を行う方法は
次のような動作を行う。

例えば、第６図に示した偏微分方程式ｄｉｖ（ｇｒａｄ
（φ））＝αにおいて、ｄｉｖはコントロールボリュー
ムの演算を意味している（第１の原理により）。そこで
、第６図に示すように、内点（ｉｊ）でのコントロール
ボリュームでは、Ｘ方向がｉ＋（１／２）から１−（１
／２）、Ｘ方向がｊ−（１／２）からｊ＋（１／２）と
すると、点（ｉｊ）の収支はその面積の平均として、次
のように離散化する。

（ａφ／ａ　Ｘ）１弓、、−（ａφ／　ａ　ｘｌ＝−シ
、Ｊ△Ｘ Δｙ ・　・　・　・　・　・　・　・　・　（１）これが、
コントロールボリュームの演算を意味する微分演算子に
対する離散化である。

さらに、ｇｒａｄは勾配の演算を意味するものであるた
め、その離散化、例えば中心差分（ａφ／ａ×）＝（φ
ｉ＋ｉ−φｔ４）　／△ｘ、（ａφ／む）＝（φ司−φ
＝−Ｌ）／△ｙを用いて離散化し、次の離散式を得る。

・　・　・　・　・　・　・　・　・　・　・　（２）
このように、第１および第２の原理を用いて得られた内
点の離散化は、コントロールボリューム方式による離散
化となる。

第３の原理、つまり境界条件の自動的な取り込みに対し
では、次のような処理を行う。

境界上の点の離散化の際に、第１の原理によりその微分
演算子の意味を把握し、第２の原理によりコントロール
ボリュームの演算を意味する微分演算子の離散化を行う
。次に、第３の原理によりその離散式中にある境界の勾
配演算を意味する項を自動抽出して、さらにこれを境界
条件に置換える。この第３の原理により、境界における
離散化もコントロールボリューム方式によるものと同一
となる。

〔実施例〕

以下、本発明の原理および実施例を、図面により詳細に
説明する。

第５図は１本発明におけるコントロールボリューム方式
の原理説明図である。

最初に、偏微分方程式の持つ保存則を近似した離散化方
式であるコントロールボリューム方式について詳述する
。

前述のように、解析領域をコントロールボリュームと呼
ぶ小領域に分割し、各小領域における収支を考えた偏微
分方程式の離散化を行う方法をコントロールボリューム
方式と呼ぶ。その場合に、各離散化点を１つのコントロ
ールボリュームが取り囲むように、かつ解析領域全体を
互いに重なり合わない複数個のコントロールボリューム
に分割する。

第５図において、破線で囲まれた領域が２次元の解析領
域におけるコントロールボリューム４０３である。点Ａ
−Ｄ（４０４〜４０７０：対す７．ｍｌントロールボリ
ュームが示されている。例えば、内点Ａ（４０４）、Ｂ
（４０５）では、各点を囲むΔＸ、Δｙを１辺とするコ
ントロールボリュームを、また境界上の点Ｃ（４０６）
では、境界を１辺とする内点の半分の面積のコントロー
ルボリュームを。

また角点Ｄ（４０７）では、境界を２辺とする内点の１
／４の面積のコントロールボリュームを、それぞれ定義
する。これにより、全離散点のコン１〜ロールボリユー
ムの総和が、解析領域に一致する。

そして、偏微分方程式を各々のコントロールボリューム
にわたって積分することにより離散式を算出する。積分
により得られた離散式は、各コントロールボリュームに
おける偏微分方程式の示す物理量の収支を示しているこ
とになる。コントロールボリュームの総領域は、解析領
域自体であるため、全コントロールボリュームから得ら
れた離散式は、解析領域に対する偏微分方程式が示す収
支を表わしている。つまり、このことが偏微分方程式の
持つ保存則を満たしているということである。

第６図および第７図は、それぞれコントロールボリュー
ム方式における離散化の考え方を示す図である。

第６図は、第５図の内点Ａ（４０４）に対するコントロ
ールボリューム（５０１）での離散化例を示すものであ
る。第５図に示した解析領域（４０１）で、偏微分方程
式ｄｉｖ（ｇｒａｄ（φ））＝α（式５０３）が成立す
る場合、内点Ａを含むコントロールボリューム（５０１
）においてもこの式（５０３）が成立する。いま、φを
変数、αを定数とし、φの収支を考える。点Ａに対する
コントロールボリューム（５０１）は、Ｘ方向に△ｘ　
（５０４）、Ｘ方向に△ｙ（５０５）の大きさを持ち、
座標上、Ｘ方向はｘ　ｉ−４からＸＬ＋奇まで、Ｘ方向
はｙニー春からｙ工→までがそれぞれ定義される。

このコントロールボリューム内で、偏微分方程式（５０
３）を積分する。以下、変形式で表わす。

ｄｉｖ　（ｇｒａｄ　（φ））＝α（５０３）→／　ｃ
ｖｄｉｖ　（ｇｒａｄ　（φ））ｄｓ＝ｆｏｖα・・・
・・　（３）ガウスの発散定理より、 →　ｆ　　　　ｇｒａｄ（φ　）　−ｎｄｕ＝　　ｒｘ
　　＊　　Δ　Ｘ　Δ　ｙ　　・　拳　−（４）ここで
、−ｎはノーマル方向のベクトルであり、→（ａφ／　
ａ　Ｘ）ｉ＋４Ｊ・Δｙ−（ｃｌφ／ａＸ）□−÷４・
△ｙ＋（ａφ／ａｙ）０．、＋１・△Ｘ−（ａφ／ａｙ
）、４−１・△ｘ＝α・ΔＸ△ｙ・・・・ｌｌ・・・・
拳・　（５）上式（５）の左辺の各項は、φの各辺ｄｆ
、　ｂｈ、　ｂｄ。

ｈｆにおける収支（５０２）を示していることがわかる
（第６図参照）。

上式（５）はコントロールボリュームに対する収支式で
あり、点Ａはこのコントロールボリュームの平均と考え
ると、離散式は上式を面積ΔＸ・△ｙで割った式となる
。

△Ｘ 各微分演算子をテーラ−展開による中心差分により離散
化し、次式（７）を得る。

・・・・・・・・・・・　（７） これが、コントロールボリューム方式により得られる内
点の離散式となる。

次に、第７図において、境界点Ｂ（４０５）に対するコ
ントロールボリューム（６０１）での離散式例を示す。

コントロールボリューム（６０１）内では、偏微分方程
式ｄｉｖ　−（ｇｒａｄ　（φ））＝α（６０２）が、
ｘ＝ｉの境界上では境界条件ｇｒａｄ　（φ）＝β（６
０３）が成立するものとして、内点と同じようにφを変
数、αを定数として、φの収支を考える。

点Ｂに対するコントロールボリューム（６０１）は、Ｘ
方向にΔｘ／２（６０４）、Ｘ方向にΔｙ（６０５）の
大きさを持ち、座標上、Ｘ方向はＸ工からＸ　ｉ＋４ま
で、Ｘ方向はｙ、、−圭から：！／ｊ＋４までが、それ
ぞれ定義されている。

このコントロールボリューム（６０１）内で偏微分方程
式（６０２）を積分する。内点の場合と同じように式を
変形すると、次式（８）を得る。

（δφ／　ａ　Ｘ）ｉ＋Ｌｊ　ｊΔｙ−（ａ　φ／ａｘ
）ｔ−４ａ・△ｙ＋＜ａφ／　ａ　ｙ）ｉ＝丹・Δｘ／
２−（ａφ／　ａｙ　）　ｘ　Ｊ−ｉ・△ｘ／２＝α、
・・・・・・・・・・・・・　（８）これもまた、左辺
の各項はφの各辺ｄｆ、ａｃ。

ｄｃ、ｆａにおける各収支（６０７）を表わしている。

これらのうち、辺ａｃにおける収支（ａφ／ａＸ）＝Δ
ｙは、境界条件に他ならない。従って、上式（８）は、
次のように書き換えることができる。

Ｃａ　φ／　ａ　ｘ）＝、ｈΔｙ−β・Δｙ　＋　Ｃａ
　φ／　ａ　ｙ）ｉ＝４・△ｘ／２−（ａφ／δｙ）よ
Ｊ−１・△ｘ　／　２　＝　（α□ΔＸΔｙ）／２　・
・・・・・・・・・・・・・　（９）従って、点Ｂに対
する離散式は、式（９）をさらに面積で割り、各微分演
算子をテーラ−展開による中心差分により離散化した式
となる。すなわち、・・・・・・・・・・・・　（１０
） 以上がコントロールボリューム方式による自動離散化で
あり、またその生成された離散式を使用した精度のよい
シミュレーションプログラムの自動生成方法である。

第２図は、本発明の一実施例を示すプログラム生成方法
の全体の処理フローチャートである。

入力情報（１０１）を与えると、処理装置により計算精
度の高い数値計算プログラムを自動生成する。その手順
として３つの流れがあって、第１の流れで、入力情報（
１０１）の中の解析領域（１０６）とメツシュ情報（１
０８０）を入力とした領域分割処理（１０２１）により
、離散点座標値データ（１０８）を求める。第２の流れ
では、この離散点座標値データ（１０８）と入力情報（
１０１）の偏微分方程式（１０４）、境界条件（１０５
）および離散化規則（１０７）を入力として離散式生成
処理（１０２）により、離散式１（１２３）を生成する
。ここでは、コントロールボリューム方式により計算精
度の高い離散式を生成する。第３の流れでは、この離散
式１（１２３）と、入力情報（１０１）中の解析対象変
数（１０９）を入力としたプログラム作成処理（１２５
）により、生成プログラム（１０３）を出力する。これ
らのうち、領域分割処理（１０２１）とプログラム作成
処理（１２５）は、公知の方法を用いて処理されるが、
離散式生成処理（１０２）は本発明による新しい方法で
処理される。

従って、以下は主として１、本発明の特徴とする部分で
ある離散式生成処理プログラム（１０２）について、詳
述する。

最初に、第２図における入力情報（１０１）について、
説明する。

第１３図は、第２図におけるプログラム入力情報を示す
図である。

２次元矩形領域で、第１３図に示す物理モデル（１３０
）を解く場合、その下に記載された入力情報（１３１）
を与える。先ず、メツシュ情報（ＭＥＳ　Ｈ）（１０８
０）は、メツシュ刻みを示すものであり、ＭＥＳＨ文で
指示する。ここでは、Ｘ方向はＯから５まで１刻み、Ｙ
方向はＯから４まで１刻みが示されている。解析領域情
報（１０６）は、解析領域がＸ方向が１〜５、Ｙ方向が
１〜４であることをＦＲＥＧＩＯＮ文で指示する。境界
条件（１０５）は、ＲＥＧＩＯＮ文とＢ　ＣＯＮ　Ｄ文
を用いて、式とその成立境界部分を示す。例えば、第２
種境界条件δＵ　／　ａ　Ｘ　＝αが境界座標Ｘ＝１゜
Ｙ＝１〜４の範囲で成立する場合、先ずこの領域範囲に
ＲＥＧＩＯＮ文でＬＥＦＴという名称を付け、次にＢＣ
ＯＮＤ文でＬＥＦＴ領域にＤＸ（ｕ）＝α、つまりａＵ
／δＸ＝αが成立することを指示する。偏微分方程式（
１０４）と解析対象変数（１０９）は、５ＯＬＶＥ文で
指示される。離散化規則（１０７）は、第１３図では例
示されていないが、式で与えることができる。また、離
散化規則（１０７）が与えられていない場合には、シス
テムのデフォルトで規定した離散化規則を用いる。

次に、本発明の出力情報について、詳述する。

第１４図は、本発明の出力情報である生成プログラム例
を示す図である。

この例では、第１３図に示す入力情報例に対する出力情
報の生成プログラム（１０３）を示している。生成ＦＯ
ＲＴＲＡＮ（１４０）は、５つのＤ○ループから構成さ
れる。これは、第１３図に示す物理モデルに示した４辺
の境界および内点に対するコードである。例えば、第１
４図のＤｏループ（１４１）は、ｙ＝１の境界に対する
ＦＯＲＴＲＡＮコーＴＲ間じ＜Ｄｏループ（１４２）は
内点に対するＦＯＲＴＲＡＮコーＴＲ間る。￥￥Ｃ０Ｎ
５には、離散式の定数、￥￥ＭＴＲＸＩから￥￥ＭＴＲ
Ｘ５には、離散式のφ０，１−□、φト、＋Ｊ＋φ□、
。

φ、ヤ０，１．φｉ＋Ｊ＋□の係数を入れる。全雛散点
について、上述のような離散式の定数と係数を入力する
コードを生成する。

これらの￥￥Ｃ０Ｎ５および￥￥ＭＴＲＸ１〜￥￥ＭＴ
ＲＸ５を、第１４図に示す連立１次方程式（１４５）に
して、これを解くことによりφの解を求めることができ
る。以上が、第２図における出力の生成プログラム（１
０３）の概要である。

次に、第２図に示す各処理の流れについて、詳述する。

第２図における領域分割処理（１０２１）は、入力情報
（１０１）である解析領域（１０６）およびメツシュ情
報（１０８０）より離散点を求め、離散点座標値データ
（１０８）を生成するものである。

領域分割処理（１０２１）の次に、第２図の離散式生成
処理（１０２）を実行する。

第１図は、本発明の一実施例を示す離散式生成処理のフ
ローチャートである。

前述のように、領域分割処理（１０２１）で生成された
離散点座標値データ（１０８）と、入力情報のうちの偏
微分方程式（１０４）、境界条件（１０５）および離散
化規則（１０７）を入力として、最初にコントロールボ
リューム離散化処理（１１１）を行う。

この処理は、偏微分方程式の中からコントロールボリュ
ームを示す微分演算子を自動的に抽出し、離散化する。

目的の微分演算子の抽出は、偏微分方程式の物理的な意
味を自動的に推定することによって可能であり、この詳
細なフローは後述する。

また、離散化処理は、全離散点について、各々与えられ
た離散化規則（１０７）がある場合にはその規則に従い
、ない場合にはシステムがデフォルトとして規定した離
散化規則に従う。なお、各離散点は、各コントロールボ
リュームの代表点として離散化が行われる。コントロー
ルボリュームの大きさ、位置についても、離散化規則中
に規定される。これらの離散化規則の詳細についでは、
後述する。

コントロールボリューム離散化処理（１１１）により生
成された中間コードを、以下展開式１％式％） 次に、全離散点の座標値（１０８）に対して、盾散点振
り分は処理（１１５）を行う。すなわち、全離散点の中
から、境界をコントロールボリュームに含まれる離散点
と内点とに振り分ける。そして、前者の境界の点におけ
る展開式１（１１３）についでは、境界条件置換処理（
１１６）を行う。この処理（１１６）は、展開式１（１
１３）中にある境界の勾配を示す項を自動的に抽出し、
与えられた境界条件（１０５）に置き換える処理である
。この置き換え処理（１１６）により生成された中間コ
ードを展開式２（１１８）と呼ぶ。

次に、境界を含むコントロールボリュームを代表する離
散点からは展開式２（１１８）を入力とし、また内点に
ついでは展開式１（１１３）を入力として、勾配離散化
処理（１２１）を行う。この処理（１２１）は、入力さ
れた展開式（１１３，１１８）から微分演算子を抽出し
、コントロールボリューム離散化処理（１１１）と同じ
ように、与えられた離散化規則（１０７）、またはデフ
ォルトとして規定した離散化規則に従って離散化を行う
（１２１）。

なお、ここで展開される微分演算子は、全て勾配を示す
ものである。

以上により、第１図に示すコントロールボリューム方式
を用いた離散化処理、つまり離散式生成処理（１１０）
が終了する。これによって得られた中間コードは、離散
式１（１２３）と呼ばれる。

第１Ｓ図は、第１図における中間コード・離散式１の一
例を示す図である。

入力情報（１０１）の偏微分方程式（１０４）に対する
離散式（１５０）は、離散式生成処理（１１０）により
、次のようになる。

φ、ヤ□−φ、　　φ、−φ１−１ △Ｘエヤ、＋ΔＸ工 そして、これが図示したテキスト形式の中間コード・離
散式１（１２３）となる。

再び第２図に戻って、離散式生成処理（１０２）の次の
処理であるプログラム作成処理（１２５）について、詳
述する。この処理（１２５）は、離散式１（１２３）と
入力情報（１０１）の１つである解析対象変数（１０９
）を入力として、生成プログラム（Ｌ　０３）を出力す
る処理である。

第１６図は、第２図におけるプログラム作成処理の概要
を示す図である。

解析対象変数（１０９）と離散式１（１２３）とを入力
として、始めに、括弧展開処理（１６１）を行う。この
処理（１６１）は、入力された離散式１（１２３）に対
して、対象変数（１０９）について括弧展開する。この
括弧展開処理（１６１）が済んだ離散式を離散式２（１
６２）と呼ぶ。

第１７図は、第１６図における離散式１と離散式２の一
例を示す図である。

第１７図に示す離散式１　（１２０）の内容（１５０）
は、第１５図に示したものと同じであるが、これに対す
る離散式２（１６２）は第１７図の式（１７０）となる
。対象変数（１０９）φについて、括弧展開した結果と
なっている。

次に、第１６図に示すように、この離散式２（１６２）
に対して、定数抽出処理（１６３）を行う。

これは、離散式２（１６２）の中から対象変数（１０９
）を含まない項を抽出する処理である。例えば、第１７
図の例の場合、離散式２（１７０）中の項−ａｔ（１７
１）がこれに該当する。さらに、抽出した項をプログラ
ム出力処理１（１６４，）により、定数項としてプログ
ラムに出力する。最後に、第１６図の処理（１６５）〜
（１６７）に示すように、解析対象変数（１０９）の各
添字ごとに（１６５）、係数抽出処理（１６６）を行い
、さらにプログラム出力処理２（１６７）を行う。例え
ば、第１７図の離散式２の例（１７０）では、φ□ヤ、
の係数として（１７２）を、φ□の係数として（１７３
）を、φ１−□の係数として（１７４）をそれぞれ抽出
し、プログラム（１０３）に出力する。

第１７図の生成ＦＯＲＴＲＡＮ（１７００）は、離散式
２（１７０）に対する出力である。￥￥Ｃ０Ｎ５（Ｉ）
に定数を、￥￥ＭＴＲｘにφ、、０．φ□。

φ、−０の係数をそれぞれ入力するコードを、生成プロ
グラムとして出力する。

以上で、第２図の第３の処理であるプログラム作成処理
（１２５）の概要を、第１６図と第１７図を用いて説明
した。

次に、第２図に示す３つの処理（１０２１，１０２，１
２５）のうち１本発明の特徴的部分である離散式生成処
理（１１０）について、詳述する。

初めに、微分演算子の抽出、つまりコントロールボリュ
ームまたは勾配を示す微分演算子の自動抽出処理につい
て、述べる。この処理は、第１図におけるコントロール
ボリューム離散化処理（１１１）にて、コントロールボ
リュームを示す微分演算子を選択することにより実施し
ているり第８図は、本発明における判定規則を示す図で
ある。本発明の選択方法は、第８図の判定規則により自
動判定することを可能にしている。

ｄｉｖ、　ｇｒａｄ、　ｒｏｔの各微分演算子は、それ
ぞれ物理的な性質を示している。そこで、式中に現われ
た場合、全てｄｉｖはコントロールボリュームを、ｇｒ
ａｄとｒｏｔは勾配を示すものと定義した（規則８１）
。

これらより展開されたａ　／　ａ　Ｘ　、δ／ａｙ、ａ
／δ２も同じである（規則８２）。

次に、δ／ａＸ、ａ／δｙ、ａ／ａｚが式中にそのまま
現われた場合には、以下に述べる理由により規則（８２
）に示すような定義を行う。

ｉ）境界条件式に現われる微分演算子δ／　ａ　Ｘ　。

ａ／ａｙ、ａ／ａｚは、全て１階微分項の中に現われる
。例えば、ａφ／δＸ＝α等の形式となる。

これは、境界勾配に他ならない。従って、この微分演算
子ａ／ａｘ、ａ／ａｙ、ａ／ａ　ｚを、勾配と定義する
。

次に、２階の偏微分方程式について示す。これは、一般
によく知られているように、 ｄｉｖ（ｋ　、　ｇｒａｄ（φ）＋Ｖ、φ）＝α・・・
（１２）という形に帰着される。従って、偏微分方程式
％式％ と記述されていても、各々がｄｉｖまたはｚｒａｄに対
応するとして、次のように定義する。

１１）２階微分項の内側のδ／ａＸは、勾配と定義する
。また。

１ｉｉ）２階微分項の外側のａ　／　ａ　Ｘは、コント
ロールボリュームと定義する。

１ｖ）１階微分項のａ　／　ａ　Ｘは、コントロールボ
リュームと定義する。

ここでは、Ｘ方向のみ述べたが、Ｘ方向、Ｚ方向につい
ても同じである。また、ｖ２についでは、ｖ２＝δ”／
ａｘ２と展開して、上記の定義を適用するものとする。

ｖ）１階の偏微分方程式の場合も同じように、これはｄ
ｉｖ（φ）＝αという形式に帰着される。従って、ａφ
／　ａ　Ｘ　＋　ａφ／　ａ　ｙ　＝α、と記載された
場合には、このａ／ａｘ、ａ／δｙは全てコントロール
ボリュームと定義する。

以上述べた規則（８１，８２）に従うことにより、コン
トロールボリュームを示す微分演算子と勾配を示す微分
演算子の自動判別を可能とし、その振り分けを第１図に
示すコントロールボリューム離散化処理（１１１）で実
施する。

次に、離散化規則の詳細について詳述する。なお、離散
化規則は、第１図に示すコントロールボリューム離散化
処理（１１１）で、コン１〜ロールボリユームを示す微
分演算子の離散化を行う際と、勾配離散化処理（１２１
）で勾配を示す微分演算子の離散化を行う際に、参照す
る。

本実施例においでは、ユーザから離散化規則の入力がな
い場合、システムでデフオル１−のに数代規則を規定し
ている。

第９図は、システムにおけるデフォルトの離散化規則を
示す図である。

微分演算子はコントロールボリュームを示すもの（９１
）と、勾配を示すもの（９２）の２種を用意して、離散
化規則をそれぞれ内点（９３）と境界（９４）について
定義する。

コントロールボリュームを示すもの（９１）について、
内点（９３）では（Ｘ方向のみを考えると）次式とする
。

ａφ、　　　φ、＋１−φ、−★ □＝□・・・・・・・・（１３） ａＸ　　　　△Ｘ工＋Δｘ２ これは、第９図において（図では、２次元で示されてい
る）、点ｉｊを取り囲む点線で示した通り、離散点から
各方向の半メツシュずつの大きさを持つものである。ま
た、境界上の点（９４）では、その境界を１辺とするコ
ントロールボリュームを定義し、離散式は１次式として
いる。

また、勾配（９２）について、内点では、テーラ−展開
で求めた中心差分公式 ％式％ を使用する。また、境界上の点（９４）についでは、境
界勾配ａφ／ａｙは離散化せずに、境界条件で置き換え
るようにする。ただし、置き換えるべき境界条件が存在
しない場合に限り、従来の方法と同じく、テーラ−展開
で求めた片側差分公式を使用するように定義°シている
。

また、コントロールボリュームにおいでは、離散化規則
を変更することにより、シミュレーションの対象となる
物理現象の範囲が広がり、また勾配についでは、離散化
規則を導出する際のテーラ−展開式の数を増加すること
により、より高精度な離散化規則を与えることができる
。

次に、第２図における歴数式生成処理（プログラム生成
処理）（１０２）の詳細について、内点の場合と境界上
の点の場合とに分けて、説明する。

第１０図は、内点の離散点における処理およびデータの
流れを示す図である。

ここでは、均一メツシュを用いて説明するが、不均一メ
ツシュでも処理は同じである。−例として、偏微分方程
式ｖ２φ十ｋ（ａφ／ａＸ）＝ａ（１０００）が１次元
領域の内点ｊにおいて成立した場合の処理を示す。

初めに、偏微分方程式（１０００）に対して、第１図に
示すコントロールボリューム離散化処理（１１１）を行
う。この処理（１１１）は、次の２つの処理から成り立
っている。第１の処理は、偏微分方程式（１０００）の
中からコントロールボリュームを示す微分演算子を抽出
する処理（１００１）である。抽出方法は、第８図に示
した定義に従う。

すなわち、第１項ｖ２φは（ａ　／　ａ　ｘ）（ａ　φ
／　ａ　ｘ）と書き換えられ、従って、外側のａφ／　
ａ　Ｘがコントロールボリュームを示すものとする。ま
た、第２項のｋ（ａφ／ａＸ）は、２階の偏微分方程式
中の１階微分項であるため、この微分演算子ａ／ａＸも
コントロールボリュームを表わしている。

コノ処理（１００１）テは、これら２つ（１）ａ／ａｘ
に対して、勾配を示す微分演算子と区別するため。

￥印を付けて次のように表わす。

￥ａＣａφ／ａｘ）／ａｘ＋に￥（ａφ／ａｘ）＝ａ　
　（１００２）・・・・・・・・・・　（１７） 第２の処理では、この式（１００２）に対して、コント
ロールボリュームを示す微分演算子を離散展開する（１
００３）。これは、先に￥印を付けた微分演算子に対し
てのみ、第９図に示したシステムの離散化規則 あるいはユーザの指定した離散化規則を適用し、・・・
・・・・・・・・・　（１８） を得る。これは、第１０図の中間コード・展開式１（１
１３）に該当する。ここでは、均一メツシュを使用する
ため、△Ｘ工＝△ｘ２＝ΔＸとしている。

次に、第１図の離散点振り分は処理（１１５）では、現
在、内点ｉの処理を行っている際中であるため、境界条
件置き換え処理（１１６）は行わずに、次の勾配離散化
処理（１２２）を行う。この処理は、第１０図に示すよ
うに２つの処理から成り立っている。第１の処理は、展
開式１（１１３）の勾配を示す微分演算子を展開する処
理（１００４）である。

すなわち、展開式１（１１３）中の微分演算子は全て勾
配を示すものであるため、これらを第９図に示す離散化
規則 に従って離散化する。ただし、ユーザが離散化規則を指
定してあれば、それに従う。展開式１の離散化の結果１
次式（１００５）を得る。

φ、や□−φ□　φ□−φ、−０ 次に、第２の処理として、上式（１９）つまり第１０図
中の式（１００５）における添字１／２を解消し、式を
まとめる（１００６）。その結果、離散式１　（１２３
）として、次式を得る。

以上が、内点における離散式の生成処理手順である。

第１１図および第１２図は、第１図における境界点に対
する処理の詳細説明図である。

先ず、第１１図により、処理の手順を述べる。

第１番目に、内点の場合と同じように、コントロールボ
リューム離散化処理（１１１）を行う。すなわち、偏微
分方程式の中からコントロールボリュームを示す微分演
算子を抽出して、これを展開し、中間コード展開式１を
生成する。次に、第１図に示す離散点振り分は処理（１
１５）により、境界条件置き換え処理（１１６）を行う
。この詳細を次に述べる。

先ず、展開式１の中の境界勾配を示す項ａφ／ａｘｅ（
点ｉは境界を示す）を抽出し、変数φをテーブル１にＲ
１＆する（１１０１）。テーブル１は、境界勾配が現わ
れている変数、例えばａφ／　ａ　Ｘｌが存在するなら
ばφを登録するテーブルであり。

後の境界条件の置き換え処理で使用する。そして、テー
ブル１には、式中に現われる全ての境界勾配から各々の
変数を登録する。ただし、同一変数が現われた場合には
、登録は１回だけにする。例えば、展開式１が、 （ｇｒａｄ（φｓ））　、　　（ｇｒａｄ（φ２））　
ｉ＋（ｇｒａｄ（φ□））ｔ”ｋ＝　（ｇｒａｄ（φｘ
　））　ｉ＋ｘである場合、φ、とφ２を１個ずつ格納
する。なお、この場合、ｊが境界、ｉ＋１は内点てあり
、φ。

は登録しない。

以下、第１１図の処理（１１，０３）〜（１１０８）を
、テーブル１に登録した全ての境界勾配のチエツクが終
了するまで繰り返して行い、中間コード展開式２を生成
する（１１０２）。これらの処理（１１０３）　〜（１
１０８）ハ、テーブル１に登録した境界勾配を境界条件
に置き換える処理である。

これについて、詳述する。先ず、テーブル１に登録した
変数の１つ、例えばφ、に注目して、境界条件の中から
φ□よの境界勾配、この場合にはＸ方向の境界勾配を含
むものを見つけ出す（１１０３）。

もし、見つかった場合には（１１０４）、展開式１の中
の境界勾配積ａφ□／δＸを境界条件式で置き換える（
１１０５）。そして、テーブル１のφ１に対して、変換
済チエツクを行う（１１０６）。例えば２展開式１が（
ｇｒａｄ（φ１　）ｉ　−（ｇｒａｄ　（φ２）、＋・
・・＝であり、境界条件が（ｇｒａｄ（φｘ　）　＝　
＝　（ｇｒａｄ　（φ２）０、＋αの場合には、展開式
１の（ｇｒａｄ（φ、））、を境界条件式の右辺に置き
換え、次のようにする。

（（ｇｒａｄ（φ２））Ｌ＋　α）−（ｇｒａｄ（φ２
））、＋・・・；・・・また、展開式１中に複数の（ｇ
ｒａｄ　（φ、））、が存在した場合には、全てを境界
条件に置き換えるようにする。そして、次に置き換えた
境界条件式中に、テーブル１に登録しである他の変数の
境界勾配が存在しているならば、これもまたテーブル１
に変換済チエツクを行う（１１０７）。上記の例の場合
、変数φ２がこれに該当する。

処理（１１０３）で、境界条件が見つかった場合には、
以上の処理（１１０５）〜（１１０７）を行い。

見つからなかった場合には、テーブル１に登録したこの
変数φ１に対して変換なしチエツクを行う（１１０８）
。このように、テーブル１に登録した全ての変数につい
て境界条件を探し、それがあれば置き換えて、変換あり
チエツク、または変換なしチエツクを行う。以上の処理
（１１０２）〜（１１０８）は、第１図における境界勾
配積の置換処理（１１６）である。

この処理（１１６）の終了後、第１図に示す勾配離散化
処理（１２１）を行う。これは、次の３つの処理より成
る。すなわち、先ず得られた式の数式処理を行う（１１
０９）。これは、数式中に相殺する項があれば、これを
相殺する。例えば、先の（（ｇｒａｄ　（φｚ））ｔ　
＋　（ｚ　）　　　（ｇｒａｄ（φ、））ｉ・・・−＝
に対して、ｇｒａｄ　８φ２）□を相殺して、α＋・・
・・＝の式を生成するにこで得られたものが展開式２で
ある。次に、内点の場合と同じく、展開式２に対して、
勾配を示す微分演算子を展開を行い（１１１０）、添字
の１／２の解消を行い、式をまとめて離散式１を生成す
る（１１１１）。以上の処理が、境界点に対する離散式
の生成処理の手順である。

なお、本来、境界勾配を示す項は、全て境界条件で置き
換えられる必要がある。従って、離散式１を生成した時
点で、置き換えがなく、また数式処理後も式中にδφｔ
／ａＸの項が存在する場合には、ワーニング（警報）を
出力する。

最後に、第１２図により、境界上の離散点の処理の一例
を説明する。

欄外に示した（１２０１）の偏微分方程式ｖ２φ□−ｋ
ｖ２φ２＝α、および座標からＸ＝Ｏでの境界条件ａ　
φ、／ａ　ｘ＝β＋ｋｃａ　φ−／ａ　ｘ）が成立する
２次元領域を考える。ｘ＝Ｏの角ではない離散点（ＩＩ
Ｊ）における離散式の生成を例示する。コントロールボ
リュームは、その下に示した（１２０２）となる。生成
手順は、表（１２０３）内に記載されている。ここの例
では、Ｘ方向のみ示している。

■コントロールボリュームを表わす微分演算子を抽出し
、￥マークを付加して、￥ａ（ａφ、／ａＸ）／ａｘ＋
に￥ａｃａφｚ／ａｙｃ）／ａｘ＝ａを得る。

■この￥のマークの付加されたコントロールボリューム
を示す微分演算子を展開して。

（（ａ　φ１／　ａ　ｘ）　ｉ＋）　　ｋ　（ａφ２／
ａ　ｘ）　、＋４）−（（ａφ、／ａｘ）ｉ　　ｋ　（
ａφｚ／ａｘ）＝’ｔ　＝ａ＝を得る。

■上記■で得られた式の中に存在する境界勾配を示す変
数を、テーブル１に登録する。ここでは、変数φ、とφ
２を登録する。

■テーブル１に登録した変数のうち、先ずφ１に注目し
て、境界条件を探す。

境界条件Ｃａ　φ、／　ａ　ｘ）ｒ＝　ｆｌ＝＋　ｋ　
Ｃａ　φ２／ａＸ）ｔが存在したので、この右辺を上記
０式の（ａφ、／ａＸ）□に代入し、 とする。

また、テーブル１の中の変数φ１に対して変換済チエツ
クを行う。

■上記■で置き換えた境界条件の右辺に変数φ２に対す
る境界勾配が含まれていたため、テーブル１中のφ２に
対しても変換済みチエツクを行う。

■テーブル１にＷＢした全ての変数φ０．φ２に対する
チエツクが終了したので、得られた式の数式処理を行い
、次式を得る。

・・・・・・・・・・・・・・・　（２１）■勾配を示
す微分演算子の展開、および添字１／２の解消を行い、
次式を得た後、処理を終了する６以上が、内点および境
界上の点における離散式生成処理（１０２）の手順であ
る。両者とも、このようにして離散式１（１２３）を生
成した後、第２図に示すように、プログラム生成処理（
１２５）により生成プログラム（１０３）を生成する。

その手順は、先の第１６図および第１７図に示した通り
である。

本実施例においでは、コントロールボリューム方式によ
る数値シミュレーションプログラムの自動生成が可能と
なる。従って、保存性のある精度の高いプログラムが生
成される。

特に、本実施例では、境界点の離散化において、あでは
まる境界条件が存在しない場合、ワーニングを出力する
ことにより、ユーザに偏微分方程式と境界条件の不適合
の存在を知らせることができる。さらに、その場合でも
、片側差分により離散化することにより、離散式の生成
を可能にしているので、ユーザの必要とするプログラム
を生成することができる。

第１８図は、本発明による数値結果例を示す図である。

物理問題として、熱拡散間Ｍ（１８１）を考える。

次元は１次元、解析領域はＸ方向に○〜１で、メツシュ
刻みは０．１、解析対象変数φに対して偏微分方程式ａ
φ／ａｔ＝ｖ２φ、そして、境界条件ａφ／δＸ＝φ（
ｘ＝ｏ）、ａφ／ａｘ＝　　’ｔ’（ｘ　＝　１　）、
および初期条件φ＝１を解く。なお、時間刻みをｔ＝ｏ
、ｏ０２５とする。

第１８図内の表には、ｔ＝０．５０００の時の結果が示
されている（２００回目、１８２）。Ｘ＝Ｏからｘ＝０
．５までの解析解と、本発明を用いたコントロールボリ
ューム方式および従来方式で求めた解とが示されている
。また、コントロールボリューム方式と従来方式の各値
の解析解に対する誤差も合わせて９百分率で示しである
。その誤差を比較すると、解析解に比して、従来方式で
は７〜８％以上を示すのに対し、本発明のコントロール
ボリューム方式では０．２％以下である。

このように、−例ではあるが、本発明の方法を用いれば
、高精度のシミュレーションプログラムの生成が可能と
なる。

特に、本実施例においでは、境界点の離散化に対して、
あではまる境界条件が存在しない場合には、ワーニング
を出力することにより、ユーザに偏微分方程式と境界条
件の不適合を通知することができる。そして、その場合
でも、片側差分によって離散化することにより、離散式
を′生成することができるので、ユーザが必要とするプ
ログラムを生成することができる。

また、特に本実施例では、コントロールボリュームを示
す微分演算に対する離散化規則をユーザが定義できるの
で、シミュレーションの対象となる物理現象を広げるこ
とができる。また、勾配を示す微分演算に対する離散化
規則もユーザが定義できるので、より高精度の離散式を
生成することができる。

さらに、偏微分方程式を指定すると、それに適合する境
界条件の形式を自動的に提示することができるので、ユ
ーザが必要な境界条件の形式を簡単に把握することが可
能であり、その結果、適切な数値計算プログラムを作成
することが可能である。

〔発明の効果〕

以上説明したように、本発明によれば、コントロールボ
リューム方式による数値シミュレーションプログラムを
自動的に生成することができるので、保存性のある、つ
まりその式が表わす物理法則の意味を保存する性質を持
った高精度なプログラムを作成することが可能となる。

【図面の簡単な説明】

第１図は本発明の一実施例を示すプログラム生成方法の
中心部分である離散式生成処理のデータの流れ図、第２
図は本発明のプログラム生成方法の全体の処理の流れを
示す図、第３図は従来の方法による離散化規則を示す図
、第４図は各微分演算子の数学的定義を示す図、第５図
はコントロールボリューム例を示す図、第６図は内点に
おけるコントロールボリューム方式の考え方を示す図。 第７図は境界におけるコントロールボリューム方式の考
え方を示す図、第８図は本発明における各微分演算子の
物理的意味の定義を示す図、第９図は本発明の一実施例
を示す離散化規則例の図、第１０図は内点の離散式生成
例を示す図、第１１図は境界における離散式生成手順を
示す図、第１２図は境界点における離散式生成例を示す
図、第１３図はプログラム方法の入力情報例を示す図、
第１４図は同じく出力である生成プログラム例を示す図
、第１５図は中間コードである離散式１の例を示す図、
第１６図はプログラム作成処理の概要を示す図、第１７
図はプログラム作成処理のデータの流れを示す図、第１
８図は本発明の数値結果例を示す図である。 １０２：離散式生成処理、１１１：コントロールボリュ
ーム離散化処理、１１５：離散点振り分は処理、１１６
：境界条件置換え処理、１２１：勾配離散化処理、１０
１，１３１：入力情報、１０８：離散点座標値データ、
１１３：展開式１．１１８：展開式２．１２３，１５０
：離散式１．１７０：離散式２．１０３：生成プログラ
ム、１０２１：領域分割処理、１２５ニブログラム作成
処理、４０２゜４０４．４０５，４０６．４０７：離散
点、４０１：解析領域、４０３，５０１，６０１：コン
トロールボリューム、１０４，５０３，６０２＝偏微分
方程式、６０３：境界条件、１３０：物理モデル、１４
０．１７００、：生成ＦＯＲＴＲＡＮ、１４２：内点に
対するコード、１４５：連立−次方程式、１５０：離散
式の内容。 第　　　１　　図 第　　　４　　　図 い？ 第　　　６　　　図 ８ｘ〜５０牛 偏微分方程式　ａｌｖ（ｇｒａｄ（＄）　）　−α第　
　　７　　　図 φ−−−−ゆ△ｘ／２　〜６０４ 偏微分方程式　ａｉｖ（ｇｒａｄ（φ））−α〜６０２
境界条件式　　ｇｒ已（φ）−β〜６０３第　　８　　
　図 第　　　９　　　図 第　　　１３　　図 〆一− 諌 区　　　　　　　　　　− 〇

Claims (1)

【特許請求の範囲】
1. １、物理現象を表わす解析領域の内点では偏微分方程式
   を、該解析領域の境界では偏微分方程式及び境界条件式
   を、それぞれ成立する微少領域内においてテーラー展開
   により求めた離散化公式に従がい離散式を導くことによ
   り、該離散式をコード化してシミュレーションプログラ
   ムを生成するプログラム生成方法において、上記物理現
   象を表わす偏微分方程式と該物理現象を解析する解析領
   域と該解析領域の境界条件と該解析領域の離散化に必要
   な離散化規則と該解析領域の代表点である離散点の座標
   値とを入力し、先ず該解析領域を上記座標値の各離散点
   を含む互いに重複しない微小領域に分割し、上記偏微分
   方程式の中から上記微小領域の演算を意味する微分演算
   子を抽出し、該微分演算子を上記離散化規則に従って離
   散化することにより中間コードである展開式を得、次に
   上記解析領域の境界部分に対してのみは特別な演算を行
   わせるため、該展開式の中の境界勾配の演算を意味する
   項を抽出し、これを上記境界条件で置換し、次に上記展
   開式中の勾配演算を意味する微分演算子を上記離散化規
   則に従って離散化することにより離散式を生成し、最後
   に該離散式をコード化してプログラムを自動生成するこ
   とを特徴とするプログラム生成方法。