IP Force 特許公報掲載プロジェクト 2022.1.31 β版

知財求人 - 知財ポータルサイト「IP Force」

▶ 国立大学法人名古屋大学の特許一覧

<>
  • 特開-FDTD解析方法 図1
  • 特開-FDTD解析方法 図2
  • 特開-FDTD解析方法 図3
  • 特開-FDTD解析方法 図4
  • 特開-FDTD解析方法 図5
< >
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】公開特許公報(A)
(11)【公開番号】P2023170345
(43)【公開日】2023-12-01
(54)【発明の名称】FDTD解析方法
(51)【国際特許分類】
   G06F 17/13 20060101AFI20231124BHJP
【FI】
G06F17/13
【審査請求】未請求
【請求項の数】11
【出願形態】OL
(21)【出願番号】P 2022082029
(22)【出願日】2022-05-19
(71)【出願人】
【識別番号】504139662
【氏名又は名称】国立大学法人東海国立大学機構
(74)【代理人】
【識別番号】100165663
【弁理士】
【氏名又は名称】加藤 光宏
(72)【発明者】
【氏名】梅田 隆行
(72)【発明者】
【氏名】関戸 晴宇
【テーマコード(参考)】
5B056
【Fターム(参考)】
5B056BB03
(57)【要約】
【課題】 電磁波および音波の伝搬を解析するFDTD法において精度向上および計算時間短縮の両立を図る。
【解決手段】 電磁波および音波の伝搬を解析するFDTD法において、奇数階微分項を付加した差分式を用いる。6次精度のFDTD法の場合は、4次精度の3階微分項×α、2次精度の5次微分項×βを付加した差分式とすることができる。係数α、βは、クーラン数Cと相関を有しており、分散関係式によって伝搬する波の位相速度を求め、理論値との誤差が最小となるように探索的に設定する
こうして得られた差分式を用いて解析することにより、精度向上を図りつつ、クーラン数を大きくすることができ、演算の所要時間を短縮することができる。
【選択図】 図3
【特許請求の範囲】
【請求項1】
コンピュータによって、FDTD法により電波または音波の伝搬を表す2つの方程式に含まれる2つの状態量を数値計算するFDTD解析方法であって、
前記コンピュータが実行する工程として、
(a)空間内に設定された各格子点における2つの前記状態量その他の初期条件を入力し、前記コンピュータのメモリの所定領域に格納するステップと、
(b)前記方程式の第一式を表す差分式によって、所定の時間刻みΔt後の前記格子点における第1の状態量を計算し、結果を前記コンピュータのメモリの所定領域に格納するステップと、
(c)前記方程式の第二式を表す差分式によって、所定の時間刻みΔt後の前記格子点における第2の状態量を計算し、結果を前記コンピュータのメモリの所定領域に格納するステップと、
(d)所定の終了条件に至るまで前記ステップ(b)、(c)を繰り返し実行するステップとを備え、
前記差分式は、FDTD法の計算式に、奇数階微分に係数を乗じた付加項を1つ以上加算した式であるFDTD解析方法。
【請求項2】
請求項1記載のFDTD解析方法であって、
前記差分式は、FDTD法におけるn(nは偶数)次精度の計算式に、(n-2m)(mはn/2より小さい自然数)次以下の精度の(2m+1)階微分に係数を乗じた付加項を1つ以上加算した式であるFDTD解析方法。
【請求項3】
請求項1記載のFDTD解析方法であって、
前記付加項に乗ぜられるそれぞれの係数は、クーラン数Cと相関を有するFDTD解析方法。
【請求項4】
請求項3記載のFDTD解析方法であって、
前記ステップ(a)において、前記係数またはクーラン数Cのいずれかを入力し、
前記ステップ(b)および前記ステップ(c)の実行前に、前記係数とクーラン数との間の予め設定された相関を利用して、残余の前記係数およびクーラン数を設定するステップを備えるFDTD解析方法。
【請求項5】
請求項1記載のFDTD解析方法であって、
前記付加項に乗ぜられるそれぞれの係数は、前記伝搬における位相速度の誤差を最小とするという評価基準に基づいて設定されているFDTD解析方法。
【請求項6】
請求項1記載のFDTD解析方法であって、
前記差分式は、FDTD法における6次精度の計算式に、4次精度の3階微分に係数を乗じた付加項を加算した式であるFDTD解析方法。
【請求項7】
請求項6記載のFDTD解析方法であって、
前記差分式は、さらに、2次精度の5階微分に係数を乗じた付加項を加算した式であるFDTD解析方法。
【請求項8】
請求項1記載のFDTD解析方法であって、
前記差分式は、FDTD法における4次精度の計算式に、4次精度の3階微分に係数を乗じた付加項を加算した式であるFDTD解析方法。
【請求項9】
請求項1記載のFDTD解析方法であって、
前記差分式は、FDTD法における4次精度の計算式に、2次精度の3階微分に係数を乗じた付加項を加算した式であるFDTD解析方法。
【請求項10】
コンピュータによって、FDTD法により電波または音波の伝搬を表す2つの方程式に含まれる2つの状態量を数値計算するためのコンピュータプログラムであって、
(a)空間内に設定された各格子点における2つの前記状態量その他の初期条件を入力し、前記コンピュータのメモリの所定領域に格納する機能と、
(b)前記方程式の第一式を表す差分式によって、所定の時間刻みΔt後の前記格子点における第1の状態量を計算し、結果を前記コンピュータのメモリの所定領域に格納する機能と、
(c)前記方程式の第二式を表す差分式によって、所定の時間刻みΔt後の前記格子点における第2の状態量を計算し、結果を前記コンピュータのメモリの所定領域に格納する機能と、
(d)所定の終了条件に至るまで前記機能(b)、(c)を繰り返し実行する機能とをコンピュータに実現させることができ、
前記差分式は、FDTD法における計算式に、奇数階微分に係数を乗じた付加項を1つ以上加算した式であるコンピュータプログラム。
【請求項11】
FDTD法により電波または音波の伝搬を表す2つの方程式に含まれる2つの状態量を数値計算するためのFDTD解析装置であって、
空間内に設定された各格子点における2つの前記状態量を記憶する算出結果保持部と、
前記各格子点における前記状態量その他の初期条件を入力し、前記算出結果保持部の所定領域に格納する入力部と、
前記方程式の第一式を表す差分式によって、所定の時間刻みΔt後の前記格子点における第1の状態量を計算し、結果を前記算出結果保持部の所定領域に格納する第1の演算部と、
前記方程式の第二式を表す差分式によって、所定の時間刻みΔt後の前記格子点における第2の状態量を計算し、結果を前記算出結果保持部の所定領域に格納する第2の演算部と、
所定の終了条件に至るまで前記第1の演算部および第2の演算部による計算を繰り返し実行する制御部とを備え、
前記差分式は、FDTD法における計算式に、奇数階微分に係数を乗じた付加項を1つ以上加算した式であるFDTD解析装置。



【発明の詳細な説明】
【技術分野】
【0001】
本発明は、FDTD(Finite-Difference Time-Domain)法による解析方法に関する。
【背景技術】
【0002】
電磁波伝搬及び音波伝搬の数値解析の方法として、FDTD法が広く用いられている。FDTD法について電磁波を例にとって概要を説明する。
【0003】
電磁波の伝搬を記述するMaxwell方程式は次式で与えられる。
【数1】
【0004】
FDTD法は、Maxwell方程式を時間及び空間差分の形で表したものであり、第一式、即ち磁界Bについて2次精度で近似した式は次式となる。
【数2】
【0005】
また第二式、即ち電解Eについて2次精度で表した式は次式となる。
【数3】
【0006】
ここで、E、Bの右上に示したt等の符号は、時刻tにおける電界(E)、磁界(B)であることを意味し、E、Bの右下に示したx、y、zは、それぞれ電界強度及び磁界強度の方向を表している。また、Δx、Δy、Δzは、それぞれx、y、z方向の空間格子の刻み幅を表し、Δtは時間積分の刻み幅を表している。
上記2つの式において、x、y、zをそれぞれy、z、xと置き換えることによりy成分についての関係式を得ることができる。また、x、y、zをそれぞれz、x、yと置き換えることによりz成分についての関係式を得ることができる。
2次精度とは、Δt、Δx、Δy、Δzの2乗の精度で表した式であることを意味する。
【0007】
Maxwell方程式の第一式を4次精度で表した次式も知られている。
【数4】
【0008】
同じく第二式を4次精度で表した式は次式となる。
【数5】
【0009】
FDTD法については、種々の改良が提案されており、例えば、特許文献1は、解析対象となる領域を複数のブロックに分割して演算することにより、解析の所要時間を短縮する方法を提案している。
【先行技術文献】
【非特許文献】
【0010】
【非特許文献1】K. S. Yee, Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media, IEEE Transactions on Antennas and Propagation, AP-14 (1966) 302-307.
【非特許文献2】P. G. Petropoulos, Phase error control for FD-TD methods of second and fourth order accuracy, IEEE Transactions on Antennas and Propagation, 42 (1994) 859-862.
【非特許文献3】E. Ikata and G. Tay, Finite-difference time domain acoustic-wave algorithm, Il Nuovo Cimento D, 20 (1998) 1779-1793.
【特許文献】
【0011】
【特許文献1】特開2006-139723号公報
【発明の概要】
【発明が解決しようとする課題】
【0012】
FDTD法では、数値不安定を生じることなく演算を可能とするための時間刻みΔtあたりの電磁波あるいは音波の伝搬速度の上限値が決まっており、これをクーラン条件とよぶ。空間格子の刻み幅が等しい場合(Δx=Δy=Δz)を仮定すると、2次精度のFDTD法の伝搬速度の上限値を表すクーラン数Cは、3次元では以下の式で与えられる。
【数6】
【0013】
また、4次精度の場合は、クーラン数は次式で与えられる。
【数7】
【0014】
このように、空間差分の精度を上げると、時間刻みの上限値は小さくなるため、電磁波あるいは音波の伝搬を所定の時間まで解析するために必要な時間ステップ数(時間積分の回数)が増える。すなわち、空間差分の精度を上げることによって計算時間が非常に長くなるという課題があった。
本発明は、この課題に鑑み、FDTD法による解析において、計算精度の向上と計算時間の短縮を両立させることを目的とする。
【課題を解決するための手段】
【0015】
本発明は、
コンピュータによって、FDTD法により電波または音波の伝搬を表す2つの方程式に含まれる2つの状態量を数値計算するFDTD解析方法であって、
前記コンピュータが実行する工程として、
(a)空間内に設定された各格子点における2つの前記状態量その他の初期条件を入力し、前記コンピュータのメモリの所定領域に格納するステップと、
(b)前記方程式の第一式を表す差分式によって、所定の時間刻みΔt後の前記格子点における第1の状態量を計算し、結果を前記コンピュータのメモリの所定領域に格納するステップと、
(c)前記方程式の第二式を表す差分式によって、所定の時間刻みΔt後の前記格子点における第2の状態量を計算し、結果を前記コンピュータのメモリの所定領域に格納するステップと、
(d)所定の終了条件に至るまで前記ステップ(b)、(c)を繰り返し実行するステップとを備え、
前記差分式は、FDTD法の計算式に、奇数階微分に係数を乗じた付加項を1つ以上加算した式であるFDTD解析方法。
【0016】
発明者は、FDTD解析方法について研究を重ねた結果、上述の形の式を用いることにより、係数およびクーラン数を調整すれば、精度向上と計算速度の短縮の両立を図ることができることを見いだした。上述の式は、semi-Lagrange法の考え方をFDTD解析に応用することにより着想されたものである。ただし、semi-Lagrange法では偶数階微分項も付加するのに対し、偶数階微分項には波形を減衰させる数値散逸効果があることを考慮し、本発明では奇数階微分項のみを付加している点に特徴を有する。
FDTD法の精度および付加項の精度は種々の組み合わせが可能である。また、付加項は1つだけである必要はなく、2以上の項を付加してもよい。
【0017】
本発明において、
前記差分式は、FDTD法におけるn(nは偶数)次精度の計算式に、(n-2m)(mはn/2より小さい自然数)次以下の精度の(2m+1)階微分に係数を乗じた付加項を1つ以上加算した式としてもよい。
【0018】
研究の結果、差分式に用いる点数を考慮すると、上述した形式が最も効率が良いことが見いだされた。n=4とすれば、差分式は「4次精度のFDTD法+2次精度の3階微分」となる。n=6とすれば、差分式は「6次精度のFDTD法+4次精度の3階微分+2次精度の5階微分」となる。他にも種々の差分式を設定可能である。
なお、本発明の差分式は、上述した形式に限定されるものではなく、例えば、6次精度FDTD法+2次精度3階微分項、2次精度FDTD法+4次精度3階微分項などの形としてもよい。
【0019】
本発明において、
前記付加項に乗ぜられるそれぞれの係数は、クーラン数Cと相関を有することが好ましい。
【0020】
研究の結果、係数とクーラン数とを相関を有する状態で調整することにより、精度と計算速度を両立できることが見いだされた。両者の相関は、関数で表すことは困難であり、具体的な数の組み合わせまたはテーブルの形式で与えられることになる。
本発明による解析方法を利用する場合には、これらの相関に基づいて、演算における空間の格子および時間刻みを、数値不安定が生じない範囲で設定して差分式を特定した上で演算を行うことになる。
【0021】
上述の通り係数およびクーラン数に相関を有する場合、
前記ステップ(a)において、前記係数またはクーラン数Cのいずれかを入力し、
前記ステップ(b)および前記ステップ(c)の実行前に、前記係数とクーラン数との間の予め設定された相関を利用して、残余の前記係数およびクーラン数を設定するステップを備えるものとしてもよい。
【0022】
上記態様によれば、例えば、空間の格子および時間刻みなどの条件を設定すれば、クーラン数が決まるため、これに応じて適切な係数の値が設定され、演算を行うことができる。逆に、係数を設定すれば、これに応じてクーラン数を設定して演算を行うことができる。
予め設定された相関は、関数、テーブルなど種々の形式で用意することができる。
【0023】
本発明において、
前記付加項に乗ぜられるそれぞれの係数は、前記伝搬における位相速度の誤差を最小とするという評価基準に基づいて設定されていることが望ましい。
【0024】
本来であれば、電磁波または音波は全波数において一定の位相速度で波形が伝搬することが理想的であるが、数値シミュレーションにおいては波数によって異なる位相速度で伝搬するという現象が現れる。上記態様によれば、この位相速度の誤差が最小となることを評価基準として係数を設定するため、計算精度を向上させることができる。
なお、「評価基準として」とは、厳密に誤差最小を実現する係数でなくても良いことを意味する。例えば、複数の係数候補について位相速度を評価した上で、その中で誤差が最小のものを選択する方法をとってもよい。
【0025】
先に説明した通り、本発明においては、FDTD法の精度および付加項には種々の組み合わせが可能である。
例えば、
前記差分式は、FDTD法における6次精度の計算式に、4次精度の3階微分に係数を乗じた付加項を加算した式としてもよい。
この場合、
前記差分式は、さらに、2次精度の5階微分に係数を乗じた付加項を加算した式としてもよい。
【0026】
また、前記差分式は、FDTD法における4次精度の計算式に、4次精度の3階微分に係数を乗じた付加項を加算した式としてもよい。
さらに、前記差分式は、FDTD法における4次精度の計算式に、2次精度の3階微分に係数を乗じた付加項を加算した式としてもよい。
【0027】
発明者の研究成果によれば、少なくともこれらの組み合わせにおいては、適切な係数およびクーラン数が設定でき、計算精度および計算速度を両立できることが確認されている。これらの態様における具体的な差分式の例および係数、クーラン数の組み合わせについては後述する。
【0028】
本発明は、上述したFDTD解析方法としての態様に限らず、かかる解析を実現するための解析装置、コンピュータによってかかる解析を実現するためのコンピュータプログラム、およびかかるプログラムを記録したコンピュータ読み取り可能な記録媒体として構成することもできる。ここでコンピュータ読み取り可能な記録媒体としては、フラッシュメモリ、ハードディスク、光ディスクなど種々の媒体を用いることができる。
【図面の簡単な説明】
【0029】
図1】解析装置の構成を示した説明図である。
図2】解析処理のフローチャートである。
図3】3次元における伝搬速度の誤差の計算結果を示す説明図である。
図4】2次元における伝搬速度の誤差の計算結果を示す説明図である。
図5】2次元における電磁波の伝搬の様子を示す説明図である。
【発明を実施するための形態】
【0030】
A.解析装置の構成:
図1は、解析装置の構成を示した説明図である。スーパーコンピュータに、図示する各機能を実現するコンピュータプログラムをインストールすることによって構成した。解析装置は、複数のコンピュータ、サーバ等をネットワークで接続して構成してもよい。
以下、電磁波の伝搬を解析する場合を例にとって説明する。
【0031】
以下、図中の各機能ブロックについて説明する。
入力部11は、運動解析に要する種々の条件を入力する。例えば、空間格子の設定および各格子点における磁界および電界の初期条件、時間刻みなどの条件が挙げられる。
【0032】
算出結果保持部12は、入力部11から入力された種々の初期条件を保持するとともに、電磁波の伝搬の解析を行う過程で算出される種々の計算結果を一時的に保持するバッファである。
磁界演算部14は、与えられた条件下で後述するFDTD法の差分式を用いて磁界の変化を演算する。計算結果は、逐次、算出結果保持部12に格納しながら計算を進める。
電界演算部15は、与えられた条件下で後述するFDTD法の差分式を用いて電界の変化を演算する。計算結果は、逐次、算出結果保持部12に格納しながら計算を進める。
磁界演算部14、電界演算部15は、別の機能部として示したが、両者を統合して構成しても差し支えない。
出力部13は、以上で得られた解析結果を出力する。出力形式は問わない。
【0033】
係数等設定部16は、磁界演算部14および電界演算部15の演算に用いられる係数を設定する。後述する通り、本実施例では、磁界および電界の演算には、FDTD法の差分式に、所定の微分項に係数を乗じた付加項を加算した差分式を用いている。係数等設定部16は、この付加項に乗ぜられる係数を設定するのである。
係数は、後述する通り、FDTD法を演算する際のクーラン数と相関を有している。クーラン数Cは、c(光速)と時間刻みΔt、空間の格子間隔Δx(=Δy=Δzとする)によってC=c×Δt/Δxで算出される。即ち、空間の格子および時間刻みが初期条件として指定されればクーラン数は定まることになる。そして、クーラン数が定まれば、それと相関のある係数も定まることになる。
本実施例では、種々のクーラン数に対して係数を求めた結果を設定テーブル17として予め用意しておいた。従って、係数等設定部16は、クーラン数に応じて、設定テーブル17を参照することにより係数を定めることができるのである。
なお、差分式で用いる係数等を初期条件の一つとしてオペレータが入力する方法をとってもよい。この場合は、係数等設定部16および設定テーブル17は省略することができる。
【0034】
B.解析処理:
図2は、解析処理のフローチャートである。図1に示した解析装置が実行する処理である。
解析処理を開始すると、解析装置は、初期条件を入力する(ステップS1)。入力される解析条件は、入力部11(図1参照)において説明した通りである。
次に、解析装置は、差分式の係数等を算出する(ステップS2)。係数等の算出方法は、係数等設定部16(図1参照)において説明した通りである。係数も入力する場合には、ステップS2の処理は省略してもよい。
【0035】
そして、解析装置は、差分式を用いて磁界を演算し、その結果を格納する(ステップS3)。また、差分式を用いて電界を演算し、その結果を格納する(ステップS4)。上記演算を終了条件が満たされるまで(ステップS5)、繰り返し実行して、結果を出力し(ステップS6)、解析処理を終了する。
終了条件としては、予め設定した時間後の結果が得られた場合、計算時間が所定の制限値を超えた場合、その他エラーが生じた場合などとすることができる。
【0036】
C.差分式:
以下、実施例の解析において用いられる差分式について具体例をあげて説明する。本実施例では、FDTD法における計算式に、奇数階微分に係数を乗じた付加項を1つ以上加算した式を用いている。付加項は、1つでも2以上でも良い。
このように差分式は、精度および付加項の選択によって、多様な組み合わせが可能であるが、以下では、2次元および3次元の場合に、それぞれ係数とクーラン数の相関の検討結果が得られているものを具体例として示す。
【0037】
C1.3次元6次精度+4次精度3階微分項+2次精度5階微分項:
かかるケースの場合、Maxwell方程式の第一式、即ち磁界を表す差分式は、次の通りとなる。
【数8】
【0038】
ここで、T1(4,3)は、4次精度3階微分項を表し、次式で与えられる。次式におけるαは係数である。
【数9】
【0039】
また、T2(2,5)は、2次精度5階微分項を表し、次式で与えられる。次式におけるβは係数である。
【数10】
【0040】
また、Maxwell方程式の第二式、即ち電界を表す差分式は、次の通りとなる。
【数11】
【0041】
ここで、T3(4,3)は、4次精度3階微分項を表し、次式で与えられる。次式におけるαは係数である。
【数12】
【0042】
また、T4(2,5)は、2次精度5階微分項を表し、次式で与えられる。次式におけるβは係数である。
【数13】
【0043】
以上の式では、x成分のみを示しているが、各式において、x、y、zをそれぞれy、z、xと置き換えることによりy成分についての関係式を得ることができる。また、x、y、zをそれぞれz、x、yと置き換えることによりz成分についての関係式を得ることができる。3次元については、以下同様である。
【0044】
付加項の係数α、βは、まず数値不安定を生じないよう設定する必要がある。FDTD法の差分式をフーリエ変換することで、周波数ωとクーラン数、時間刻Δt、格子間隔Δxなどの関数として波動の分散関係式が得られるが、この分散関係式において周波数ωが虚数解となるとき、数値不安定が生じることになる。従って、係数α、βは、分散関係式において周波数ωが虚数解とならないことを条件の一つとして設定する必要がある。
また、電磁波は全波数において一定の位相速度で波形が伝搬することが理想的であるため、分散関係式から得られる位相速度、即ちωの実部を波数κで除した値、と理論値の2乗平均誤差をκについて積分した値が最も小さくなるよう係数を設定することが好ましい。本実施例では、これらの条件に基づいて、クーラン数および係数の組み合わせを探索的に求めることとした。以下で説明する種々の差分式においても同様である。
【0045】
本ケースの場合、分散関係式は次式で与えられる。
【数14】
【0046】
空間格子がx、y、z方向に等しい間隔(Δx=Δy=Δz)の場合、この分散関係式に基づいて得られるクーラン数および係数の組み合わせは、次の表の通りである。Δx、Δy、Δzが等しくない場合は、x方向、y方向、z方向のそれぞれに対して、下表を求めることになる。
【表1】
【0047】
図3は3次元における伝搬速度の誤差の計算結果を示す説明図である。横軸にクーラン数をとり、縦軸に伝搬速度についての演算結果と理論値との誤差を示した。FDTD(2,2)は従来の2次精度のFDTD法による演算結果を示している。クーラン条件を満たす範囲、即ちクーラン数Cが約0.577より小さい範囲でしか解が得られない。FDTD(2,4)は、従来の4次精度のFDTD法による演算結果を示している。同じくクーラン条件を満たす範囲、即ちクーラン数Cが約0.494より小さい範囲でしか解が得られない。
【0048】
Presentは、本実施例における結果である。従来のFDTD法による結果とは異なり、クーラン数Cが1.0以下の範囲で解が得られている。クーラン数が1.0より大きい範囲であっても、それに応じた適切なα、βを探索できれば解を得ることは可能である。このように、本実施例の方法によれば、従来のFDTD法よりも、クーラン数を大きくすることができる。クーラン数は、cΔt/Δxで求まるから、空間格子の間隔Δxが一定であれば時間刻みΔtを大きくすることができることを意味する。
本実施例の方法では、付加項の分、1回の演算量が増大するものの、本実施例によれば、従来のFDTD法よりも時間刻みを大きくすることができ、全体として計算時間を従来の2次精度のFDTD法と同等にすることが可能となる。
【0049】
また、従来のFDTD法で解が得られている範囲で比較すると、本実施例の演算結果は、従来のFDTD法よりも全体に誤差が小さくなっていることがわかる。即ち、本実施例の方法によれば、従来のFDTD法よりも精度を向上させることが可能となる。
【0050】
C2.2次元6次精度+4次精度3階微分項+2次精度5階微分項:
かかるケースの場合、Maxwell方程式の第一式、即ち磁界を表す差分式は、x成分について次の通りとなる。
【数15】
【0051】
また、y成分について次の通りとなる。
【数16】
【0052】
かかるケースの場合、Maxwell方程式の第二式、即ち電界を表す差分式は、次の通りとなる。
【数17】
【0053】
以上の3式において、係数αが乗ぜられた項が4次精度3階微分の付加項であり、係数βが乗ぜられた項が2次精度5階微分の付加項である。
本ケースの場合、係数α、βおよびクーラン数を探索するための分散関係式は次式で与えられる。
【数18】
【0054】
空間格子がx、y方向に等しい間隔(Δx=Δy)の場合、この分散関係式に基づいて得られるクーラン数および係数の組み合わせは、次の表の通りである。Δx、Δyが等しくない場合は、x方向、y方向のそれぞれに対して、下表を求めることになる。
【表2】
【0055】
図4は2次元における伝搬速度の誤差の計算結果を示す説明図である。横軸にクーラン数をとり、縦軸に伝搬速度についての演算結果と理論値との誤差を示した。FDTD(2,2)は従来の2次精度のFDTD法による演算結果、FDTD(2,4)は従来の4次精度のFDTD法による演算結果、Presentは本実施例による演算結果を示している。3次元の場合と同様、本実施例によれば、クーラン数を大きくすることができるため、
全体として計算時間を短縮することが可能となることがわかる。また、従来のFDTD法よりも精度を向上させることも可能となる。
【0056】
図5は、2次元における電磁波の伝搬の様子を示す説明図である。初期条件として原点(x=0,y=0)でピークを有するガウス関数を与えた場合の所定時間経過後の電磁波の伝搬の結果を表している。
図5(a)は、従来の2次精度のFDTD法をクーラン数0.5で実行した結果である。図5(b)は、従来の4次精度のFDTD法をクーラン数0.5で実行した結果である。図5(c)は、本実施例の方法をクーラン数0.5で実行した結果である。図5(d)は本実施例の方法を、クーラン数1.0で実行した結果である。
図示する通り、クーラン数0.5での本実施例の結果(図5(c))は、従来の2次精度(図5(a))および4次精度(図5(b))とほとんど差違のない結果が得られている。クーラン数1.0での本実施例の結果(図5(d))では、x方向およびy方向に、若干の数値振動による縞模様が現れているが、それ以外は、2次精度(図5(a))および4次精度(図5(b))と大きな差違はない。即ち、本実施例の方法によれば、クーラン数を1.0まで大きくしても、看過できないほどの誤差を与えることなく演算結果が得られることがわかる。
【0057】
C3.3次元6次精度+4次精度3階微分項:
かかるケースの場合、Maxwell方程式の第一式、即ち磁界を表す差分式は、次の通りとなる。
【数19】
【0058】
また、Maxwell方程式の第二式、即ち電界を表す差分式は、次の通りとなる。
【数20】
【0059】
以上の2式において、係数αが乗ぜられた項が4次精度3階微分の付加項である。
本ケースの場合、係数αおよびクーラン数を探索するための分散関係式は次式で与えられる。
【数21】
【0060】
空間格子がx、y、z方向に等しい間隔(Δx=Δy=Δz)の場合、この分散関係式に基づいて得られるクーラン数および係数αの組み合わせは、次の表の通りである。Δx、Δy、Δzが等しくない場合は、x、y、z方向のそれぞれに対して、下表を求めることになる。
【表3】
【0061】
既に説明したケースと同様、本ケースにおいても、計算精度を確保または向上しつつ、クーラン数を大きく設定することにより演算の所要時間を短くすることができる。
【0062】
C4.2次元6次精度+4次精度3階微分項:
かかるケースの場合、Maxwell方程式の第一式、即ち磁界を表す差分式は、x成分について次の通りとなる。
【数22】
【0063】
また、y成分について次の通りとなる。
【数23】
【0064】
かかるケースの場合、Maxwell方程式の第二式、即ち電界を表す差分式は、次の通りとなる。
【数24】
【0065】
以上の各式において、係数αが乗ぜられた項が4次精度3階微分の付加項である。
本ケースの場合、係数αおよびクーラン数を探索するための分散関係式は次式で与えられる。
【数25】
【0066】
空間格子がx、y方向に等しい間隔(Δx=Δy)の場合、この分散関係式に基づいて得られるクーラン数および係数の組み合わせは、次の表の通りである。Δx、Δyが等しくない場合は、x方向、y方向のそれぞれに対して、下表を求めることになる。
【表4】
【0067】
既に説明したケースと同様、本ケースにおいても、計算精度を確保または向上しつつ、クーラン数を大きく設定することにより演算の所要時間を短くすることができる。
【0068】
C5.3次元4次精度+4次精度3階微分項:
かかるケースの場合、Maxwell方程式の第一式、即ち磁界を表す差分式は、次の通りとなる。
【数26】
【0069】
また、Maxwell方程式の第二式、即ち電界を表す差分式は、次の通りとなる。
【数27】
【0070】
以上の2式において、係数αが乗ぜられた項が4次精度3階微分の付加項である。
本ケースの場合、係数αおよびクーラン数を探索するための分散関係式は次式で与えられる。
【数28】
【0071】
空間格子がx、y、z方向に等しい間隔(Δx=Δy=Δz)の場合、この分散関係式に基づいて得られるクーラン数および係数αの組み合わせは、次の表の通りである。Δx、Δy、Δzが等しくない場合は、x方向、y方向、z方向のそれぞれに対して、下表を求めることになる。
【表5】
【0072】
既に説明したケースと同様、本ケースにおいても、計算精度を確保または向上しつつ、クーラン数を大きく設定することにより演算の所要時間を短くすることができる。
【0073】
C6.2次元4次精度+4次精度3階微分項:
かかるケースの場合、Maxwell方程式の第一式、即ち磁界を表す差分式は、x成分について次の通りとなる。
【数29】
【0074】
また、y成分について次の通りとなる。
【数30】
【0075】
かかるケースの場合、Maxwell方程式の第二式、即ち電界を表す差分式は、次の通りとなる。
【数31】
【0076】
以上の各式において、係数αが乗ぜられた項が4次精度3階微分の付加項である。
本ケースの場合、係数αおよびクーラン数を探索するための分散関係式は次式で与えられる。
【数32】
【0077】
空間格子がx、y方向に等しい間隔(Δx=Δy)の場合、この分散関係式に基づいて得られるクーラン数および係数の組み合わせは、次の表の通りである。Δx、Δyが等しくない場合は、x方向、y方向のそれぞれに対して、下表を求めることになる。
【表6】
【0078】
既に説明したケースと同様、本ケースにおいても、計算精度を確保または向上しつつ、クーラン数を大きく設定することにより演算の所要時間を短くすることができる。
【0079】
C7.3次元4次精度+2次精度3階微分項:
かかるケースの場合、Maxwell方程式の第一式、即ち磁界を表す差分式は、次の通りとなる。
【数33】
【0080】
また、Maxwell方程式の第二式、即ち電界を表す差分式は、次の通りとなる。
【数34】
【0081】
以上の2式において、係数αが乗ぜられた項が4次精度3階微分の付加項である。
本ケースの場合、係数αおよびクーラン数を探索するための分散関係式は次式で与えられる。
【数35】
【0082】
空間格子がx、y、z方向に等しい間隔(Δx=Δy=Δz)の場合、この分散関係式に基づいて得られるクーラン数および係数αの組み合わせは、次の表の通りである。Δx、Δy、Δzが等しくない場合は、x方向、y方向、z方向のそれぞれに対して、下表を求めることになる。
【表7】
【0083】
既に説明したケースと同様、本ケースにおいても、計算精度を確保または向上しつつ、クーラン数を大きく設定することにより演算の所要時間を短くすることができる。
【0084】
C8.2次元4次精度+2次精度3階微分項:
かかるケースの場合、Maxwell方程式の第一式、即ち磁界を表す差分式は、x成分について次の通りとなる。
【数36】
【0085】
また、y成分について次の通りとなる。
【数37】
【0086】
かかるケースの場合、Maxwell方程式の第二式、即ち電界を表す差分式は、次の通りとなる。
【数38】
【0087】
以上の各式において、係数αが乗ぜられた項が4次精度3階微分の付加項である。
本ケースの場合、係数αおよびクーラン数を探索するための分散関係式は次式で与えられる。
【数39】
【0088】
空間格子がx、y方向に等しい間隔(Δx=Δy)の場合、この分散関係式に基づいて得られるクーラン数および係数の組み合わせは、次の表の通りである。Δx、Δyが等しくない場合は、x方向、y方向のそれぞれに対して、下表を求めることになる。
【表8】
【0089】
既に説明したケースと同様、本ケースにおいても、計算精度を確保または向上しつつ、クーラン数を大きく設定することにより演算の所要時間を短くすることができる。
【0090】
以上で説明した実施例の解析方法によれば、電磁波または音波の解析を精度よく高速に計算することができる。
【0091】
音波の伝搬の解析への適用について説明する。
音波の伝搬を記述する簡略化した流体方程式は次式で与えられる。
【数40】
【0092】
上式は、Maxwell方程式と同じ形式である。従って、実施例で説明した種々の差分式の変数を置き換えることにより、音波の伝搬の解析に適用することができる。
なお、本発明は、本実施例に示した態様に関わらず、その趣旨を変更しない範囲で種々の態様により実現可能である。
【産業上の利用可能性】
【0093】
本発明は、電磁波または音波の伝搬の解析に利用することができる。
【符号の説明】
【0094】
11 :入力部
12 :算出結果保持部
13 :出力部
14 :磁界演算部
15 :電界演算部
16 :係数等設定部
17 :設定テーブル
図1
図2
図3
図4
図5