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

Abstract

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

Description

【発明の詳細な説明】 ［産業上の利用分野コ 本発明は対象とする偏微分方程式系中に複数の変数を含
む連立偏微分方程式と境界条件群とで表される物理問題
の数値シミュレーションプログラムを上記方程式と境界
条件とから自動的に生成するプログラム生成方法に関し
、特に流体解析、構造解析のための数値シミュレーショ
ンプログラムの生成に好適なプログラム生成方法に係る
。

［従来の技術］ 特開昭６１−５０３８０１に記載のように、物理現象を
支配する偏微分方程式を入力し、有限要素法を用いて上
記物理現象を数値的にシミュレーションするプログラム
を自動生成する方法において、従来は取り扱える偏微分
方程式が単一変数のみを含む単一偏微分方程式に限定さ
れていた。即ち、例えば熱伝導現象を支配する。

ｄｉｖ　（ｋ−ｇｒａｄ　（ｆ））＝０のような単一の
偏微分方程式を変数ｆについて解く場合にのみ有限要素
法プログラムの自動生成を行っていた。

また、偏微分方程式中に現れる対急変数（解くべき変数
）の位置にも制約があり。

ｄｉｖ　（ｋ−ｇｒａｄ　（ｆ、）＋Ｖ−ｆ、）の形式
の項においてｆｌの位置に解くべき変数（対象変数）が
出現した場合にのみプログラムの自動生成を行っていた
。

さらに、時間依存の熱伝導問題を陽解法で記述した表現 ｆ＝ｆ０＋ｄ　１　ｔ　−ｄ　ｉｖ　　（ｇｒａｄ（ｆ
ｏ））などが入力された場合、ｆおよびｆＯを評価する
際に常に一般の帯行列（Ｃｏｎｓｉｓｔｅｎｔ　Ｍａｔ
ｒｉｘ）を構成するようなプログラムを自動生成してい
た。

［発明が解決しようとする問題点］ 従来は、対象となる偏微分方程式系が複数の変数を含む
連立偏微分方程式である場合、および複数の変数に対し
異なる精度で補則を行う必要のある場合について考慮が
なされていなかった。例えば流速の遅い流体現象を支配
するストークス方程式 ％式％（１） （但しｍは粘性係数、ｄｘ（）、ｄｙ（）はＸ。

ｙ方向の偏微分、ｕ、ｖは流速、ｐは圧力を表す。） は、（１）〜（３）を連立して解く必要があるが、従来
方式ではこれを解くプログラムを自動生成することがで
きなかった。さらに従来は、連立する偏微分方程式に対
し、境界条件群が与えられた場合、どの偏微分方程式と
どの境界条件が対応するかを一意に決定する配慮がなさ
れていた。

また、従来は偏微分方程式中の対象変数の出現位置が限
定されていたために、例えば移流拡散方程式 ％式％（４） などを変数ｆについて解くプロゲラ１１を自動生成する
ことが不可能であった。

さらに、従来は時間依存の熱伝導問題を陽解法で記述し
た表現 ｆ　＝ｆ　ｏ　＋　ｄ　ｌ　ｔ−ｄ　ｉｖ　（ｇ　ｒ　
ａ　ｄ　（ｆ　ｏ））　（５）に対して、ｆおよびｆＯ
を評価する際に一般の帯行列を構成するようなプログラ
ムを生成しており、帯行列を対角化して計算時間の節約
を図る手法（Ｌ　ａｍｐｊ、ｎｇ　）を用いたプログラ
ムを自動生成することはできなかった。

本発明の目的は、物理現象を支配する偏微分方程式系お
よび境界条件群などを入力し、有限要素法による数値シ
ミュレーションプログラムを自動生成する方法において
、より広範囲な偏微分方程式系に対するプログラムの自
動生成方法を提供することにある。具体的には、異なる
精度によって補間する必要のある複数の変数を含む連立
１次方程式と境界条件群に対して妥当なシミュレーショ
ンプログラムを自動生成する方法を提供すること、任意
の位置に対象変数を含む偏微分方程式に対して妥当なシ
ミュレーションプログラムを自動生成する方法を提供す
ること、および入力情報中で与えられた指定に基づいて
行列の対角化を施したシミュレーションプログラムを自
動生成する方法を提供することにある。

［問題点を解決するための手段］ 上記従来技術の問題点を解決するために、本発明のプロ
グラム生成方法は、物理現象を支配する偏微分方程式系
、境界条件群、上記物理現象の生起する領域形状および
Ｒ１１ｌ化の方法に関する指定を入力情報として与え、
上記物理現象を有限要素法を用いて計算機上でシミュレ
ーションするためのＦＯＲＴＲＡＮプログラムを自動生
成するプログラム生成方法において、連立する複数の変
数および偏微分方程式を含む上記偏微分方程式系に対し
、各偏微分方程式に対して最適な重み関数を自動的に選
択して数式的に重み関数を乗じ、数式的に部分積分を実
行し、境界条件群より各偏微分方程式に対して最適な境
界条件を自動的に選択して数式的に境界条件の導入処理
を行い、数式的に対象変数について総和と積分の順序を
交換して変数の離散化処理を行い、係数行列・定数ベク
トルへ成分を代入し線形計算によりこれを解くためのＦ
ＯＲＴＲＡＮプログラムの自動生成処理を行なうことを
特徴とする。

［作用］ 前記重み関数の自動選択処理によれば、偏微分方程式の
重み関数の次数と対象変数の補間次数とが１対１に対応
する。このことは、離散化した際の未知数の個数と連立
１次方程式の連立本数とが一致することを意味しており
、線形計算によって連立１次方程式の解を得ることが可
能となる。また、各偏微分方程式に対してどの対象変数
の寄与が太きいかを考慮して重み関数を配分するので、
得られる結果の精度面でも有利となる。

一方、上記境界条件の導入処理によれば、連立偏微分方
程式に対する物理的に自然な境界条件を境界積分項中に
最適な形で導入することができる。

また、上記変数の離散化処理は、偏微分方程式中に現わ
れる変数、定数をすべてΣＮ１ｆｉの形に展開した後、
解くべき対象変数について総和と積分の順序を交換する
ように数式的に大変形を行なう。このため、対象変数が
偏微分方程式中に出現する位置に制約がなくなる。また
、入力情報中で指定がなされた場合、離散化処理部にお
いて補間関数Ｎｉを１に置き換えることにより、生成さ
れるＦＯＲＴＲＡＮプログラムが構成する係数行列を対
角成分のみを含む対角行列にすることができる。

［実施例］ 本発明の詳細な説明に先立ち、ガレルギン法を用いた有
限要素法による偏微分方程式の正数化手順を説明する。

有限要素法を用いる際には、シミュレーション対象領域
を適当な大きさの要素を用いて分割する。

要素の形状は３辺形または４辺形（３次元では４面体、
３角柱または６面体）が一般的である。図２に要素分割
の１例を示す。各要素にはその頂点あるいは辺の中点に
節点を設ける。節点は変数の値を求めるべき代表点であ
る。即ち、有限要素法を用いた場合には、各節点に対す
る変数値の組の形で偏微分方程式の近似解が得られるこ
とになる。

１つの要素について、要素と節点との関係を第３図、第
４図に示す。

有限要素法を用いて節点上の変数値を得るためには、予
め偏微分方程式を節点上の変数値に関する連立１次方程
式に変形し、これを線形計算によって解くという手順を
経る。この連立１次方程式への変換操作を離散化と呼ぶ
。

離散化の手順を具体例を用いて示す。

偏微分方程式 ％式％（１０） を対象変数ｆについて解く場合、（１０）の両辺に重み
関数Ｎ、をかけ、その領域全体にわたる積分種を０とお
く。即ち、 ／ｖｃｉ　ｉ　ｖ（ｇ　ｒ　ａ　ｄ（ｆ））ＮＪｄ　ｖ
＝ｏ　　　（１１）とする。重み関数の例を第５図、第
６図に示す。

第５図は第３図の要素に対する重み関数、第６図は第４
図の要素に対する重み関数である。重み関数ＮＪは各節
点毎に１つ定まるので、（１１）は節点毎に成立する式
であり、節点数分連立していると考えられる。

（１１）に数学上の公式である部分積分を施し、次式を
得る。

ｆ　ｖ　ｇｒａｄ（ｆ）ｇｒａｄ（Ｎ、＋）ｄｖ−ｆ　
ｓ　ｎｇｒａｄ（ｆ）Ｎｊｄｓ＝０（ｎは境界に対する
法線方向単位ベクトル）第２項は対象領域の境界に関す
る積分であるのでこれに境界条件 ｎ　ｇ　ｒａｄ　（ｆ）　＝ｑ　　　　　　　　　　（
１３）（ｑは境界で与えられた流束） を取り込み ｆ　ｖ　ｇｒａｄ（ｆ）ｇｒａｄ（ＮＪ）ｄｖ−ｆ　ｓ
　ｑＮＪｄｓ＝ｏ　　　（１４）と変形し式中に境界条
件を導入する。

ここで、未知量ｆを既知の基底関数Ｎ、の線形結合とし
て ｆ＝Ｚ　　Ｎ五ｆ　Ｉ（１５） のように近似する。但しｆ、は節点におけるｆの値であ
る。この近似により、連続量ｆは節点における離散的な
値ｆ、に置き換えられることになる。

即ち、第５図あるいは第６図に示した関数の重ね合わせ
として全体のｆの分布を近似するわけである。（１５）
を（１４）に代入すれば、ｆ　ｖ　ｇｒａｄ（ΣＮｔｆ
＋）ｇｒａｄ（Ｎｔ）ｄｖ−ｆ　ｓ　ｑＮａｄｓ”ｏ　
　（１６）ｆ＋は領域について定数となっているので、
Σとｆの順序を交換し、かつｆ、の積分の外へ出すこと
ができて、結局 ΣＣｆｖ　ｇｒａｄ（Ｎｔ）ｇｒａｄ（Ｎｔ）ｄｖ）ｆ
ｔ＝ｆｓ　ｑＮＪｄｓ　　（１７）と変形できる。（１
７）において各被積分項は既知であるから、（１７）は
ｆ、に関する連立１次方程式となっている。（１７）を
行列形式で書けば、Ａ（ｆ）＝（ｃ）　　　　　　　　
　　　　　　　　　（１ｇ）（ａ、、（行列Ａの（ｉ、
ｊ）成分） ＝　ｆ　ｖ　ｇｒａｄ（Ｎｔ）ｇｒａｄ（Ｎｔ）ｄｖｆ
　Ｊ＝　ｆ　ｓ　ｑＮａｄｓ　　　　　）となる。（１
８）を線形計算によって解けば、求める値（ｆ）を得る
ことができる。以上がガレルギン法を用いた有限要素法
による偏微分方程式の離散化手順の説明である。

次に、変数の補間の精度について説明する。有限要素法
によって離散化を行う場合、変数の値は（１５）によっ
て近似される。換言すれば、節点以外の場所における変
数値は節点における変数値から（１５）によって補間さ
れることにぼる。従って基底関数（＝重み関数）Ｎのと
り方にとって補間の精度に差が生ずるわけである。具体
的には、第３図に示す要素３１を用い、従って第５図に
示す基底関数を用いた場合には、節点間の値は接点値か
ら２次関数的に補間されることになる。このような意味
で第３図に示すような頂点にのみ節点を持つ要素３１を
線形要素、第４図に示すような頂点および辺上の中間点
に節点を持つ要素４１を２次要素という。同様に第５図
の重み関数を線形の重み関数、第６図の重み関数を２次
の重み関数という。

一般に、要素数が等しい場合、補間の次数の高い２次要
素の方が線形要素より精度がよい。以上が有限要素法に
おける変数の補間の精度の説明である。

次に、係数行列の対角化を行うための手段について説明
する。前記（９）に対し上記ガレルキシ法を適用すれば
次式を得る。

Σ（／　ｖＮ＋Ｎａｄｖ）ｆｔ＝　・・・　　　　（２
０）ここで、左辺の補間関数Ｎ１を恒等的に１とおけば
、 Σ（／ｖＮｊｄｖ）ｆｔ＝　・・・　　　　　（２１）
となる。このように変形すれば、左辺から生ずる係数行
列はｉ＝ｊの場合のみ非零となり、ｌ≠ｊの場合には零
となる。換言すれば、係数行列が対角化されたことにな
る。以上が係数行列を対角化するための手段である。

以下、本発明の詳細な説明する。

第１図にプログラム生成処理１２の構成と、これによっ
て数値シミュレーションを行なう場合のデータの流れを
示す。

プログラム生成処理１２は、物理現象を支配する偏微分
方程式系、境界条件群、対象とする領域の形状、および
離散化の方法の指定などの情報を記述したシミュレーシ
ョン情報１１を入力し、その現象を有限要素法を用いて
シミュレ−１−するためのＦＯＲＴＲＡＮプログラム２
６を出力する。

同時に要素・節点の情報を格納する節点、要素データ２
５を出力する。このプログラムをＦＯＲＴＲＡＮ翻訳機
構２７を通して機械語プログラム２８に変換する。計算
装置に機械語プログラム２８を入力し、節点要素データ
２５をデータとして付与すれば、シミュレーション結果
１１０を得ることができる。

第８図に流路中の遅い流れをシミュレートする場合のシ
ミュレーション情報１１の例を示す。シミュレーション
情報１１は領域情報８２、境界条件群８３偏微分方程式
８４よりなる。領域情報８２は領域を特徴づける点Ａ、
Ｂ、Ｃ・・・等の座標情報８５とそれらの点を連結した
辺の構成情報８６．８７、面の構成情報８８．８９より
成る。

８６は辺ＡＢが点Ａと点Ｂを結ぶ直線であること、８７
は辺ＢＤが点Ｂ、Ｃ，Ｇ、Ｄを通る曲線であることを示
す。同時に８８は面ＡＢＣＦが点Ａ。

Ｂ、Ｃ，Ｆを頂点とする４辺形であることを示す。

境界条件群８３は、領域の境界で成立する境界条件より
成る。条件８１０は辺ＢＤ上とＡＥ上で変数ＵとＶの値
が０であること、即ち壁面で流素が０であることを規定
する。同様に条件８１１は流入口ＡＢでの流速を、条件
８１２は流出口ＤＥでの院力と流速勾配を規定している
。偏微分方程式系８４は、現象を支配する偏微分方程式
系と解くべき変数名を保持する。８１３　ｒｓＯＬＶＥ
ｕ：２　　ｖ：２　　ｐ：１　０ＦＪは続く偏微分方程
式系を変数ｕ＋Ｖ＋Ｐについて解くべき事、およびｕ、
ｖは２次の補間関数（第６図）によって補間し、ｐは線
形の補間関数（第５図）によって補間すべき事を示す。

続＜８１４，８１５，８１６は遅い流れを支配するスト
ークス方程式である。

もう１つの例として、第２５図に係数行列の対角化を指
定を含むシミュレーション情報１１（７）例を示す。こ
のうち、領域情報と境界条件群については係数行列対角
化の説明上必要がないので割愛しである。偏微分方程式
２５０１は時間依存の熱伝導問題の記述例であり、次の
ことを指定している。２５０１と２５０５は組をなす指
定であり、間に挾まれた２５０３．２５０４を２５０２
のｒＵＮＴＩ　ＬＪ以下が成立するまで繰り返すことを
指示する。ここで、ＮＴはループ１回ごとに１ずつ繰り
上がるカウンタであり１本例ではループが１００回繰り
返されることになる。２５０３は時間依存熱伝導問題の
時間に関する陽的な代入計算であり、こ話を繰り返すこ
とにより、順次各時刻での温度ｆの分布が求められる。

ここで係数行列の対角化は次のように指示する。即ち、
対角化の対象となる変数を［］で囲み、その直後に対角
化を表わすキーワード（ＬＭ）を〈　〉で挾んで記述す
る。２５０１の例では、左辺の変数ｆに対する係数行列
が対角化されるべきことを指示している。２５０４は変
数の更新を行なう代入計算であり、２５０３によって求
められたｆによって旧値ｆＯを′置き換え、次の時間ス
テップの計算を行なう準備をしている。

メツシュ分割部２４はシミュレーション情報１１を入力
し、要素分割を行って節点要素データ２５を生成する。

第９図に第８図の入力情報１１に対する節点要素データ
２５の例を示す。第８図シミュレーション情報１１中偏
微分方程式系８４に与えられた補間次数の最大値は２で
あるので、要素としては第４図に示した辺上中間節点つ
きのものを用いる。節点要素データ２５は、節点データ
９１．要素データ９２およびサブセットデータ９３より
成る。節点データ９１は各節点につけられだ通し番号（
節点番号）順に節点の座標値を列挙したものである。要
素データ９２は要素に付けられた通し番号（要素番号）
順に要素を構成する節点番号のリストを列挙したもので
ある。サブセットデータ９３は領域情報８２中の各辺２
面の名称に対し、それらの辺、面に含まれる節点の節点
番号に列挙したものである。これらの節点要素情報２５
は、ＦＯＲＴＲＡＮプログラム２６の実行時に計算′！
Ａ置２９によって配列中に読み込まれ、係数行列、定数
ベクトルの成分を数値的に決定する際に参照される。ま
た、シミュレーション結果１１０も、上記節点番号ある
いは要素番号順に出力される。

次に、本発明の中心をなす式処理部第１図１１１につい
て説明する。式処理部１１１は、シミュレーション情報
１１中の境界条件群８３と偏微分方程式系８４を読み込
み、処理に好適な２分本等の形態である穴情報Ａ１９に
変形する人前処理部１３と、穴情報Ａ１９中の補微分方
程式に対する最適な重み関数を決定し、これを偏微分方
程式に乗じて積分形に変換した穴情報Ｂ２０に変換する
重み関数決定部と１式情報Ｂ２０に部分積分を施して穴
情報Ｃ２１に変形する部分積分実行部１５と、穴情報Ｃ
２１中の積分形に変換された偏微分方程式と境界条件群
と１対１に対応付け、境界条件を導入して穴情報Ｄ２２
を作成する境界条件照合・導入部１６と、穴情報Ｄ２２
中の変数を節点における変数値の補間で置き換え、対象
変数については積分と総和の順序を交換して穴情報Ｅ２
３を作成する離散化部１７と、穴情報Ｅ２３中の係数行
列と定数ベクトルの成分にあたる式から実際の成分を数
値的に求めるＦＯＲＴＲＡＮプログラム２６を生成する
プログラム出力部１８より成る。

以下、式処理部第１図１１１の働きを順に追って説明す
る。シミュレーション情報１１の例として、第８図およ
び第２５図に挙げたものを適宜用いて説明を行なう。

人前処理部１３は、シミュレーション情報１１から偏微
分方程式系と境界条件群を順次読み込み、各式を処理が
容易な２分水形式に変形して穴情報Ａ１９を作成する。

その処理の流れを第２４図に示す。まず、シミュレーシ
ョン情報中の偏微分方程式を読み込み、偏微分方程式毎
にヘッダを作成しく２４０２）、偏微分方程式を通常の
コンパイラ等で用いられる構成解析法に基づいて構成解
析して２分木を作成しく２４０３）、この２分木をさぎ
のヘッダにつなぐ（２４０４）。次にシミュレーション
情報中の境界条件式を読み込み、同様にヘッダを作成し
く２４０６）、ヘッダに境界条件の指定された領域の名
称を格納しく２４０７）、境界条件を構成解析して２分
木を作成しヘッダにつなぐ（２４０８，２４０９）。第
８図のシミュレーション情報に対する穴情報Ａ１９の具
体例を第１０図に示す。１００２以下は偏微分方程式系
第８図８４の第１式に対する２分木、１００９以下は同
第２式に対する２分本、１０１０以下は同第３式に対す
る２分木である。同様に１０１１゜１０１２．１０１３
は境界条件群に対する２分木である。各２分木は＋、−
，ｄｉｖ、ｇｒａｄ等の演算子を格納する節と、変数名
ｕ、ｖ、数値Ｏ等を格納する末端より成る。節の具体的
なデータ構造を第２５図に示す。各節には、その節の左
側、右側の節へのポインタ、ならびに演算子名称を格納
する。また、端末の具体的なデータ構造を第２６図に示
す。末端には変数名称または数値と、変数の補間次数お
よび変数に対するオプションを格納する。、（オプショ
ンとしては、その変数に対する係数行列の対角化等の指
定を格納するが、その具体例については後述する。）各
２分木の頂点には、その２分本に対する属性を格納する
ヘッダ（１００１，１００４，１００５，１００６゜１
００７．１００８）を設ける。境界条件に対するヘッダ
には、その境界条件が成立する領域名称および境界条件
の１種、２種の区別を格納する。

（１種、２種の判別については後述する。）偏微分方程
式に対するヘッダは後に重み関数の次数と対応する変数
名称を格納するため空きフィールドとしておく。さらに
、２分木の末端に変数名がある場合、その偶数が何次の
精度で補間されるかを偏微分方程式系８４のｒ、ｕ：２
　　ｖ：２　　ｐ：ＩＪから読み取り、変数名とペアに
して格納する。もう１つの例として、第２０図のシミュ
レーション情報に対して生成される式情報Ａを第２１図
に示す。第２０図中の偏微分方程式系２５０１において
は、２５０３式の左辺変数ｆに対して係数行列の対角化
の指定（ＬＭ）が付されている。これにともない、第２
１図式情報Ａの対応する末端２６０２のオプションフィ
ールドにＬＭが格納され、その変数を離散化する際に係
数行列を対角化すべきことを示す。

次に、重み関数決定部１４は、各偏微分方程式に乗ずる
べき重み関数の次数を前述の重み関数決定方法に基づい
て決定し、第２分水のヘッダに格納する。それをもとに
偏微分方程式に重み関数を乗じて積分を行う。以下に重
み関数決定部の処理の方針を説明する。連立偏微分方程
式を有限要素法によって離散化する際には、連立偏微分
方程式中の各偏微分方程式に重み関数をかける必要があ
る。このとき、どの偏微分方程式にどの（何次の）重み
関数をかけるかを決定する必要がある。そのためにまず
、連立偏微分方程式中の各偏微分方程式と解くべき対象
変数とを１対１に対応付ける。

対応付けの規則は以下による。

■　連立偏微分方程式中の各偏微分方程式についてｄｉ
ｖ　（Ａ−ｇｒａｄ（ｆ））項（Ａは係数テンソル）が
存在し、Ａが対角成分を含む場合、未だ対応のついてい
ないｆが対象変数であればｆとその偏微分方程式とを対
応付ける。１測微分方程式中に該当する項が複数存在す
る場合には、未だ対応のついていない対象変数の内最初
に出現したものとその偏微分方程式を対応付ける。

■　■で対応のつかない対象変数と偏微分方程式が残っ
た場合、それらを対象として再び■と同様の対応付けを
行う。但し今回はテンソルＡ中に対角成分を含まない項
に着目して対応付ける。

■　更に対応のつかない対象変数と偏微分方程式が残っ
た場合には、随れらを対象として次の規則により対応付
ける。対象変数について１階微分のみを含む偏微分方程
式についてはその式に含まれない対象変数と対応をつけ
る。

■　以上で対応の付かない対象変数と偏微分方程式が残
った場合には、あらかじめ与えられた対象変数の並び順
に対応を付ける。

以上によって偏微分方程式と対象変数とを１対１に対応
付けたあと、その対象変数の補間次数と同じ次数を持つ
重み関数をその偏微分方程式の重み関数とする。以上が
連立偏微分方程式に対する重み関数決定部の処理の方針
である。

第２７図にこの処理方針に基づいた重み関数決定部の処
理の流れを示す。第１０図式情報Ａ　１．９を例にとれ
ば、処理は次のようになる。まず第１０図１００２以下
の第１式に対し規則■を適用する。第１式においてはｄ
ｉｖ節の下のｇｒａｄ節の引き数はＵであり、またｇｒ
ａｄには係数が与えられていないので係数テンソルＡは
単位テンソルであると判断できる。従って係数テンソル
Ａは対角成分を含む見なすことができ、規則■により第
１式に対応する変数はＵと判断され、ヘッダ１００１に
は変数名Ｕとその補間次数２を格納する。第２式に対し
ても同様にして対応する変数はＶであると判断され、ヘ
ッダ１００４には名称Ｖと補間次数２が格納される。第
３式はｕ、ｖ以外の変数ｐと対応付けられることになり
、ヘッダ１００５には変数名称ｐと対応付けられること
になり、ヘッダ１００５には変数名称ｐと補間次数１が
格納される。（第１式とＵ、第２式とＶを対をづけた時
点で第３式と変数ｐが残るので、規則■を用いずに直接
第３式とｐとを対応づけて良い。）以上の重み関数決定
結果に基づき、各式に重み関数を乗じ領域に関する積分
を行う。２分水上での操作は、＝節の下の十節に着目し
、その下に領域積分を表すｆｖ節９本節２重み関数を組
み合わせた部分木（第１１図１１０１）を挿入すること
によりなされる。結果として得られる穴情報Ｂを図１１
に示す。ここで、１次の重み関数をＮＪ、２次の重み関
数をＭＪで表す。

続く部分積分実行部１５では、穴情報Ｂ２０中の積分項
に対し部分積分を適用する。部分積分実行部の処理の流
れを第２９図に示す。２分水上での操作は、領域積分に
ついては重み関数の上にｇｒａ’ｄ節を加え（２９０３
）、ｆｖ節の下のｄｉｖ節を削除することによってなさ
れる（２９０５）。同時に、部分積分によって生ずる境
界積分項をヘッダの下につなぐ（２９０３）。

境界積分項は、境界積分を示すｆｓ節、＊節の下に重み
関数、および法線ベクトルとの積を示すｎ節をつなぎ、
その下に領域積分値のｄｉｖｌより下に相当する部分木
をつなぐことにより得られる。

第１１図１１０１に対してこの処理を施した結果として
得られる穴情報Ｃ２１を第１２図に示す。

次いで、境界条件照合、導入部１６において。

前述の境界条件取り込み方法に基づき、各境界領域ごと
に偏微分方程式と境界条件との対応付けを行い境界条件
を導入する。まず、境界条件導入の前処理として、境界
条件の２分本のヘッダ（第１２図１００６．１００７等
）を順次たどって境界条件の定義された領域名称をピッ
クアップし、偏微分方程式のヘッダの下にその情報をつ
ないで境界領域名称ヘッダをつくる。次に、各偏微分方
程式について、各境界領域ごとに導入すべき境界条件を
決定していく。このとき、偏微分方程式と境界条件式を
１対１に対応付ける方針を以下に説明する。

まず、与えられた境界条件を１種、２種に分類する。分
類の基準を第７図に示す。次に２種境界条件を偏微分方
程式と対応付ける。ここでは、偏微分方程式より得られ
る境界積分の被積分項と境界条件とを逐次比較し、被積
分項を全て含む境界条件が存在した場合、その偏微分方
程式と境界条件を対応付ける。どの偏微分方程式にも対
応の付かない２種境界条件が残った場合にはその境界条
件の導入は行わない。次いで、１種境界条件と残った偏
微分方程式の対応を付ける。対応付けの規則を以下に示
す。

（９１種境界条件の対象変数の補間次数と等しい次数の
重み関数がかけられた偏微分方程式と対応付ける。

（Φ　偏微分方程式中で２階微分のかけられた変数と１
種境界条件の対象変数が等しい場合、これらの式を対応
付ける。

■　偏微分方程式中で１階微分のみのかけられた変数と
１種境界条件の対象変数が等しい場合。

これらの式は対応付けない。

規則の優先順位は■、■、■の順とする。またこれらの
規則によって対応が一意に決まらない場合には、予め指
定された偏微分方程式と境界条件の並び順によって対応
を定める。以上の規則によって境界条件と偏微分方程式
の対応を付けた後、２種境界条件は偏微分方程式の境界
積分項と置換し、１種境界条件は偏微分方程式自身と置
換する。但し、１種、２種いずれの境界条件とも対応の
付かない偏微分方程式が残った場合には、その偏微分方
程式が残った場合には、その偏微分方程式に対する境界
積分項を０と置く。以上が、境界条件照合・導入部で連
立偏微分方程式に対して境界条件２を１対１に対応付け
る方針である。

以上の方針に基づいた偏微分方程式と境界条件の対応付
けの処理の流れを第２８図に示す。第１３図の例につい
て具体的に処理の流れを説明する。第１３図第１式１０
０１（７）領域ＢＤ１３０１に対して境界条件の２分水
のヘッダを探索すれば。

２種境界条件は存在しないことがわかる（第８図８３参
照）。そこで１種境界条件を探索すれば、ｕ＝ｏ、ｖ＝
ｏの２条件が指定されている。ここで第１式のヘッダ１
００　’１を参照すれば第１式に対応する変数がＵであ
ることが知れるため、第１式のＢＤに対する境界条件と
してはＵ＝Ｏが選ばれる。そこで、境界条件の２分水中
のＵ＝Ｏに相当する２分水１００６を削除し、その内容
を１３０１にコピーする。同様にして、第１式の２番目
の境界領域ＤＥについて境界条件の２分木のヘッダを探
索すれば、２種境界条件としてｎ（ｇｒａｄ　（ｕ））
＝Ｏ＋　ｎ　（ｇｒａｄ　（ｖ））＝０が指定されてい
る。前述の境界条件取り込み方法に基づき、これらが第
１式の境界積分項と置換可能かを調べる。置換可能性は
２分水上では以下のように判定できる。まず第１式のヘ
ッダから境界積分項をたどりｎ節以下に着目する。（第
１０図１００８参照）此の両者のｎ節以下が完全に一致
すれば２種境界として取り込み可能であると判断できる
。いま、ｎ　（ｇｒａｄ　（ｕ））＝０に石目すればこ
れは第１式の境界積分項に対し上記条件を満たすので取
り込み可能となる。そこで、境界積分項のｎ節より下の
部分木を境界条件２分水の＝節から右側に出た部分木（
第１０図１０１４参照）で置き換えたものを、第１式の
１３０２以下に格納する。また、ここで取り込んだ境界
条件ｎ　ｇ　ｒ　ａ　ｄ　（ｕ）　＝Ｏに対する２分本
は削除する。このような操作をすべて式の全境界領域に
ついて行えば境界条件の導入が完了する。

この状態における穴情報Ｄ２２を第１４図に示す。

ここで、対応する境界条件を満たない境界領域が残った
場合には、その境界では境界積分項二〇の２種境界条件
が成立するものとみなし、境界積分項のｎ以下をＯで置
き換えた部分木を該当する境界条件名称ヘッダにつなぐ
。逆に、取り込むことのできない境界条件が指定された
場合にはその境界条件の２分木は削除されずに最後まで
残る。このような残存する２分木は無視して次の処理に
進むか、何らかのワーニングメツセージを出力すること
も可能である。

続いて、離散化部１７において変数を節点上の変数値の
補間式によって置換し、さらに対象変数については積分
と総和の順序を交換して対象変数を積分の外側に置く処
理を行う。第３０図に離散化部１７における処理の流れ
を示す。変数の補間式による置換は、２分水上では次の
操作により実現できる。２分本の末端を探索し、それが
変数であった場合、総和を表すΣ節と本節、基底関数名
（ＮＪまたはＭＪ）および離散変数名（変数名に添字ｉ
を付したもの）より成る部分木（第１５図１５０１）で
置き換える３００２゜基底関数ＮｔとＭＩとの使い分け
は変数名と対になって格納されている補間次数を参照し
、１の場合にＮ１．２の場合にＭｌを用いる。この時点
の穴情報Ｅ２３を第１５図に示す。また、対象変数を積
分の外へ出す操作は、２分水の末端に（離散化された）
対象変数が存在した場合、それに連結する８節９本節お
よび基底関数をまとめてｆｖ節の直前へ移動することに
より実現できる３００３　（第６図１６０１）、この状
態を第１６図に示す。なお、この離散化操作は１種境界
条件の２分本については特に行う必要はない。また、第
２２図のように行列対角化指定のかけられた項２７０１
が存在する場合には、これを離散化する際に現われるＮ
。

またはＭ、を１で置き換える３００５゜これによって第
２３図のような穴情報Ｅを得ることができる。

以上で式の変形は完了し、プログラム出力部１８が最終
的な大変形結果である穴情報Ｅ２３に基づいてＦＯＲＴ
ＲＡＮプログラムの生成を行う。

プログラム出力部１８の主要な働きは、離散化された偏
微分方程式に対応する連立１次方程式、即ち係数行列と
定数ベクトルを構成するためのＦＯＲＴＲＡＮプログラ
ムを生成することである。プログラム出力部の処理方法
を説明するに先立ち、最終的に構成すべき係数行列の例
を第１９図を用いて説明する。係数行列の大きさは、式
の本数（あるいは対象変数の個数本節点数）を次元数と
する正方行列１９０２となる。この行列中に、さきの離
散化によって得られた対象変数に対する係数式（第１７
図１７１２のｆｖ節以下はｐｉに対する係数式、同様に
１６０３のｆｖ節以下はｕｉに対する係数式）を評価し
て値を代入していくわけである。例えば、第１７図１７
１２以下に現れる重み関数Ｍａのインデクスがｋなる値
をとり、基底関数Ｎｔのインデクスが１なる値をとる場
合を考えると、それを評価した値の格納場所は第１９図
中斜線位置１９０１となる。

ＦＯＲＴＲＡＮプログラム２６中の係数行列への格納は
上記の規則に従うが、１＋ＪをＤｏ−プによって動かし
、同一処理が可能な部分はまとめて処理を行う。このＤ
ｏループ生成と係数行列への格納を行うためのプログラ
ム出力部の処理方法は第１７図、第１８図に示すように
なる。

第１７図は領域積分項に対するプログラム出力方法を示
す。１７０１において係数行列を単位行列化しておく。

これは、例えば辺上中間節点（第４図において各辺の中
点に位置する節点）においては第３式（第８図８１６）
の重み関数が定義されていないので、この節点に対する
行は成分が全てＯとなり、行列が特異となることを避け
る手段である。１７０３，１７０４は全ての式、変数に
対するプログラム生成を行わせるためのループである。

ここで着目式２着目変数が定まるので、以降ではその変
数、式に対するプログラム生成を行う。１７０５におい
て係数行列を評価すべき式を式情報Ｅ２３中から探索す
る。１７０６，１７０７は１７１２と共に基底関数・重
み関数を全節点についてスキャンするためのＤｏループ
を生成する。１７０８．１７０９は変数・式が定義され
ない節点に対し評価のスキップを行っている。

１７１０で実際に数値積分を実行して係数値を評価する
ためのプログラムを生成し、１７１１で第１９図に示し
たごとくに行列の対応位置へ係数値を格納するプログラ
ムを生成する。

第３２図に第１６図式情報Ｅの第１項（ｐの係数）を評
価する際にプログラム出力部１８で生成されるＦＯＲＴ
ＲＡＮプログラムの例を示す。

３２０２は着目式が第１番目の偏微分方程式であること
、３２０２は着目変数が第３番目の対象変数（ｐ）であ
ることを表わすコードである。

３２０４．３２０５．３２０６はそれぞれ第１７図１７
１２．１７０６．１７０７によって生成されたコードで
あり、基底関数、重み関数を全節点についてスキャンす
るためのＤｏループをなす。

Ｍ　Ｅ　Ｓ　ＨＮ　Ｏは全要素数、ＭＯＤＥＳ　（Ｌ）
は第り番目の要素に含まれる節点の数である。

３２０７．３２０８にて着目節点の節点番号を得る。Ｉ
ＡＤＡＴＡ　（Ｉ、Ｊ）は工番目の要素の５番目の節点
の節点番号を格納するテーブルである。

３２０９．３２１０は基底関数、重み関数の定義されて
いない節点をスキップするための判定文である。ＩＮＴ
ＥＲ（Ｉ）は１番目の節点が辺上の中間接点の場合に１
、それ以外の場合に０を格納したテーブルである。３２
１１においてｆｖＭＪ・ｄｘ（Ｎｌ）の数値積分を行な
うサブルーチンエＦＵＮＣをコールする。ここで、５Ｆ
ＮＸはｄ−ｘ（Ｎｌ）を定義するファンクション、ＳＦ
ＭはＭＪを定義するファンクションである。またＶＡＬ
は積分結果を受は取る変数である。３２１２にて積分結
果の係数行列Ｃ０ＦＦの対応する位置への足し込みを行
なう。

第１８図は境界条件用のＦＯＲＴＲＡＮプログラム出力
方法を示す。領域積分の場合と異なり、１本の式に対し
て複数の境界領域が対応しているので、境界領域をスキ
ャンするループ１８０３が加わる。また、１種と２種の
境界条件には取り込み方法に差があるので、１８０４に
よってこれを振り分けている。２種条件に対する処理は
領域積分の場合とほぼ同様であるが、各境界における節
点番号を得るためのプログラムを生成する必要がある。

（１８０８，１８１０）一方、１種境界条件についての
プログラム出力方法は第３１図に示す。領域積分の場合
と同様にＤｏループの生成の行なった後、３１０７．３
１０８にて係数行列の対応する行の対角成分を１．非対
角成分を０とするプログラムを生成し、３１０９で定数
ベクトルの対応する位置に境界値（ｕ＝ｃにおけるＣの
値）を代入するプログラムを生成する。

以上で、本発明の１実施例の説明を終わる。

本実施例ではストークス流体および時間依存熱伝導問題
を取り上げて説明したが、本発明の特徴はむしろ対象と
する偏微分方程式を選ばないところにあり、２階以下の
偏微分方程式系で記述される他の様々な問題にそのまま
適用できるものである。また、説明を用意するために２
次元問題を取り上げたが、３次元問題についてもまった
く同様にプログラムの自動生成が可能となる。

［発明の効果］ 本発明によれば、変数を異なる精度によって補間する必
要のある場合を含む広範囲な連立偏微分方程式系に対し
て適切なＦＯＲＴＲＡＮプログラムを自動生成すること
が可能である。連立偏微分方程式が現れるのは、流体、
不純物拡散、半導体デバイス、構造などきわめて広い分
野にわたっており、しかもこれらの分野は偏微分方程式
の連立によって人手によるプログラミングの困難さも大
きい。このような分野の問題に対してプログラムの自動
生成手段を提供したことにより、プログラム開発工数と
コストの大幅な削減が可能となる。

【図面の簡単な説明】

第１図は本発明の１実施例の全体構成図、第２図は要素
分割の説明図、第３図、第４図は要素の説明図、第５図
は線形の重み関数の説明図、第６図は２字の重み関数の
説明図、第７図は境界条件種別を決定する規則の表、第
８図はシミュレーション情報の説明図、第９図は節点要
素データの説明図、第１０図は穴情報Ａの説明図、第１
１図は穴情報Ｂの説明図、第１２図は穴情報Ｃの説明図
、第１３．１４図は穴情報りの説明図、第１５゜１６図
は穴情報Ｅの説明図、第１７．１８図はプログラム出力
部の処理の流れの説明図、第１９図は係数行列構成の説
明図、第２０図は係数行列の対角化指定を含むシミュレ
ーション情報の説明図。 第２１図は係数行列の対角化指定を含む穴情報Ａの説明
図、第２２図は係数行列の対角化指定を含む穴情報りの
説明図、第２３図は係数行列の対角化指定を含む穴情報
Ｅの説明図、第２４図は大曲処理部の処理の流れの説明
図、第２５図は節のデータ構造の図、第２６図は末端の
データ構造の図、第２７図は重み関数決定部の処理の流
れの説明図、第２８図は境界条件照合部の処理の流れの
説明図、第２９図は部分積分実行部の処理の流れの説明
図、第３０図は離散化部の処理の流れの説明図、第３１
図はプログラム生成部の第１種境界に対する処理の流れ
の説明図、第３２図は生成ＦＯＲＴＲＡＮプログラムの
例。 ［符号の説明］ １１・・・シミュレーション情報 １２・・・プログラム生成機構 １３・・・大曲処理部　　１４・・・重み関数決定部１
５・・・部分積分実行部 １６・・・境界条件照合・導入部 １７・・・離散化部　　　１８・・・プログラム出力部
１９・・・穴情報Ａ　　　　２０・・・穴情報Ｂ２１・
・・穴情報Ｃ２２・・・穴情報Ｄ２３・・・穴情報Ｅ　
　　　２４・・・メツシュ分割部２５・・・節点要素デ
ータ ２６・・・ＦＯＲＴＲＡＮプログラム ２７・・・翻訳機構　　　２８・・・機械語プログラム
２９・・・計算装置 ３０・・・シミュレーション結果 ３１・・・式処理部 ・、！ 代理人　弁理士　小川勝馬、−１− 笛　／圀 箪２圓 第３面 ＄４田 処た。　　　　ガ仁 ＄ｇ面 Ａγ 菓４面 Ｊ４．ガた ＄７回 第２回 １η　　タ　　Ｅコ ；：：；都ね用榛２１．？髪絆琳・ ノａＢ　　　　　　　　　”””　　　　　　ノブ／／
　　　　　　　７“＄７２図 箋・３面 第　　７４　ｐ弓 ノ３０ノ 番／り固 ツタ１２　　傅喀友予テテ１ 、’ｆ／／ ２７ｆｆｌ ￥２３田 ２１１２ノ 第２り居

Claims (1)

1. 【特許請求の範囲】 １、物理現象を支配する偏微分方程式系、境界条件群、
   上記物理現象の生起する領域形状および離散化の方法に
   関する指定を入力情報として与え、上記物理現象を有限
   要素法を用いて計算機上でシミュレーションするための
   プログラムを自動生成するプログラム生成方法において
   、連立する複数の変数および偏微分方程式を含む上記偏
   微分方程式系に対し、数式的に重み関数を乗じ、数式的
   に部分積分を実行し、数式的に上記境界条件群より境界
   条件の導入処理を行い、数式的に変数の離散化処理を行
   い、係数行列・定数ベクトルへ成分を代入し線形計算に
   よりこれを解くためのプログラムを生成することを特徴
   とするプログラム生成方法。 ２、上記の重み関数を乗ずる処理は、上記偏微分方程式
   系中の各偏微分方程式に対し、最適な重み関数を自動的
   に選択することを特徴とする特許請求の範囲第１項記載
   のプログラム生成方法。 ３、上記の境界条件の導入処理は、上記偏微分方程式系
   中の各偏微分方程式に対し、上記境界条件群からこれに
   適合する境界条件を自動的に選択して導入することを特
   徴とする特許請求の範囲第１項記載のプログラム生成方
   法。 ４、上記変数の離散化処理は、偏微分方程式中の解くべ
   き変数（対象変数）について補間関数による総和と積分
   との順序を交換することにより行い、取り扱える偏微分
   方程式の範囲を拡大したことを特徴とする特許請求の範
   囲第１項記載のプログラム生成方法。 ５、上記変数の離散化処理は、上記入力情報中で指定が
   なされた場合、使用する補間関数を変化させることによ
   り、生成されるプログラム中で構成される係数行列を対
   角化することを特徴とする特許請求の範囲第１項記載の
   プログラム生成方法。 ６、上記プログラムの生成処理は、構成すべき係数行列
   に対して成分を代入する際にインデクスの管理を２重に
   したプログラムを生成することにより、上記偏微分方程
   式中の連立する変数、偏微分方程式に対してプログラム
   の自動生成を行うことを特徴とする特許請求の範囲第１
   項記載のプログラム生成方法。