JP2019124595A - 物体内への物流入量算出方法 - Google Patents
物体内への物流入量算出方法 Download PDFInfo
- Publication number
- JP2019124595A JP2019124595A JP2018005758A JP2018005758A JP2019124595A JP 2019124595 A JP2019124595 A JP 2019124595A JP 2018005758 A JP2018005758 A JP 2018005758A JP 2018005758 A JP2018005758 A JP 2018005758A JP 2019124595 A JP2019124595 A JP 2019124595A
- Authority
- JP
- Japan
- Prior art keywords
- solution
- thing
- simultaneous linear
- linear equations
- 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.)
- Pending
Links
- 238000004364 calculation method Methods 0.000 title claims abstract description 44
- 238000000034 method Methods 0.000 claims abstract description 81
- 238000012937 correction Methods 0.000 claims abstract description 16
- 230000005684 electric field Effects 0.000 abstract description 24
- 239000000463 material Substances 0.000 abstract description 2
- 239000012530 fluid Substances 0.000 abstract 1
- 239000013598 vector Substances 0.000 description 44
- 239000011159 matrix material Substances 0.000 description 13
- 238000004458 analytical method Methods 0.000 description 6
- 230000006698 induction Effects 0.000 description 4
- 238000002939 conjugate gradient method Methods 0.000 description 3
- 238000000354 decomposition reaction Methods 0.000 description 3
- 244000205754 Colocasia esculenta Species 0.000 description 2
- 235000006481 Colocasia esculenta Nutrition 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 230000005672 electromagnetic field Effects 0.000 description 2
- 230000001788 irregular Effects 0.000 description 2
- 238000007781 pre-processing Methods 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 239000000126 substance Substances 0.000 description 2
- 238000004891 communication Methods 0.000 description 1
- 239000004020 conductor Substances 0.000 description 1
- 230000008094 contradictory effect Effects 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 230000002950 deficient Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 238000001727 in vivo Methods 0.000 description 1
- 239000007788 liquid Substances 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
Landscapes
- Investigating Or Analyzing Materials By The Use Of Electric Means (AREA)
- Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)
Abstract
【課題】物質内に流入する電界などの物の量を正確に把握することのできる物流入量算出方法を提供する。【解決手段】スカラーポテンシャルを未知数とする連立１次方程式に、外部環境データを算入して物体表面から物体内への所定の物の流入量を算出する流入量算出工程を具備する方法であって、上記の所定の物は、自由境界を有する流動的に移動可能なものを意味し、上記流入量算出工程において、上記連立1次方程式を補正する工程を具備し、この補正は上記連立1次方程式における各行から所定の値を差し引くことにより行われることを特徴とする。【選択図】なし
Description
本発明は、物体内への物流入量算出方法に関するものである。
近年、中間周波数帯 (300 Hz – 10 MHz)の電磁界を利用した機器の利用拡大に伴い、機器から発生する電磁界による生体への影響に関心が寄せられている。準静的近似が有効な周波数帯における電界ばく露時の人体内誘導電界を計算する手法の一つとして SPFD (Scalar Potential Finite Difference)法が挙げられる（非特許文献１及び２参照）。
T.W. Dawson, J. De Moerloose, and M.A. Stuchly, "Com-parison of magnetically induced ELF fields in humans com-puted by FDTD and scalar potential FD codes," Appl.Comput. Electromagn. Soc. J., vol.11, no.3, pp.63–71, 1996.
T.W. Dawson, K. Caputa, and M.A. Stuchly, "High-resolution organ dosimetry for human exposure to low-frequency electric fields," IEEE Trans. Power Deliv., vol.13,no.2, pp.366–373, 1998.
しかしながら、上述の手法は、スカラーポテンシャルを未知数とする大規模な連立一次方程式を解くことに帰着されるため、連立一次方程式を高速に解ける数値計算ライブラリを用いて解くことで、解析時間の短縮が期待できるものの、SPFD 法の連立一次方程式は、右辺ベクトル b の全要素の和が 0 という条件を満たさない場合、解が存在しないという問題がある。この右辺ベクトルの値は人体表面から流れ込む電流から計算され、物理的に考えるとその総和は 0 になるはずである。しかし、この値を外部電界計算によって求める場合は、数値誤差やモデルの離散化による誤差などが含まれ、電流の総和は 0 にならないことが多い。
そのため、誤差が最も小さくなる妥協解を求めることが必要である。しかし、右辺ベクトルの要素和が 0 でないとき、解法によっては解が収束しないことや大きな誤差が出てしまうことがあり、数値計算ライブラリを適用出来ない場合がある。
従って本発明の目的は、物質内に流入する電界などの物の量を正確に把握することのできる物流入量算出方法を提供することにある。
そのため、誤差が最も小さくなる妥協解を求めることが必要である。しかし、右辺ベクトルの要素和が 0 でないとき、解法によっては解が収束しないことや大きな誤差が出てしまうことがあり、数値計算ライブラリを適用出来ない場合がある。
従って本発明の目的は、物質内に流入する電界などの物の量を正確に把握することのできる物流入量算出方法を提供することにある。
本発明は、以下の各発明を提供するものである。
１．スカラーポテンシャルを未知数とする連立1次方程式に、外部環境データを算入して物体表面から物体内への所定の物の流入量を算出する流入量算出工程を具備する方法であって、
上記の所定の物は、自由境界を有する流動的に移動可能なものを意味し、
上記流入量算出工程において、上記連立1次方程式を補正する工程を具備し、この補正は上記連立1次方程式における各行から所定の値を差し引くことにより行われる
ことを特徴とする物体内への物流入量算出方法。
この方法において上記外部環境データは、実測データでも計算値でもよい。また上記の所定の物の具体例としては、電流、気体又は液体等が挙げられる。
２．１記載の方法であって、上記の所定の物が電流である方法。
１．スカラーポテンシャルを未知数とする連立1次方程式に、外部環境データを算入して物体表面から物体内への所定の物の流入量を算出する流入量算出工程を具備する方法であって、
上記の所定の物は、自由境界を有する流動的に移動可能なものを意味し、
上記流入量算出工程において、上記連立1次方程式を補正する工程を具備し、この補正は上記連立1次方程式における各行から所定の値を差し引くことにより行われる
ことを特徴とする物体内への物流入量算出方法。
この方法において上記外部環境データは、実測データでも計算値でもよい。また上記の所定の物の具体例としては、電流、気体又は液体等が挙げられる。
２．１記載の方法であって、上記の所定の物が電流である方法。
本発明によれば、物質内に流入する電界などの物の量を正確に把握することができる。
以下、本発明についてさらに詳細に説明する。
本研究では、右辺ベクトルの総和が 0 でない場合におけるSPFD 法の連立一次方程式の計算方法について検討を行う。まず、様々な数値解法に適用でき、誤差が最小となる解が得られるような右辺ベクトルの補正方法を提案する。そして、補正方法の妥当性を検討するとともに、実際に補正を用いて誘導電界解析を行い、様々な数値解法への適用可能性を検討する。
2. 1 SPFD 法 [5]
SPFD 法を用いて人体内誘導電界を計算する場合、2 段階の解析を行う。はじめに、電界ばく露時に人体表面に誘起される表面電荷を計算するため、人体外部の電界分布を計算する。次に、得られた表面電荷を境界条件として用い、SPFD 法で体内誘導電界を計算する。SPFD 法では、人体をボクセルでモデル化し、生体内のボクセルに含まれる全ての節点について、スカラーポテンシャルを未知数とする以下の方程式を立てる。
式 (1) は図１に示す節点0とその節点に隣接する6つの接点r(r =1∼6)のスカラーポテンシャルからなる方程式となる。
なお、Srは節点0−r間のアドミタンス、ψr は節点 r におけるスカラーポテンシャルを表している。すなわち、SPFD法は図２のような回路網にモデル化される。また、境界条件は外部電界計算によって求められた表面電荷密度分布ρs を用いて、電流連続の式により式 (2) と表される。
なお、ˆnはモデル表面に垂直な法線単位ベクトルJ は電流密度を意味する。この境界条件は図２の回路網では電流源として表される。
本研究では、右辺ベクトルの総和が 0 でない場合におけるSPFD 法の連立一次方程式の計算方法について検討を行う。まず、様々な数値解法に適用でき、誤差が最小となる解が得られるような右辺ベクトルの補正方法を提案する。そして、補正方法の妥当性を検討するとともに、実際に補正を用いて誘導電界解析を行い、様々な数値解法への適用可能性を検討する。
2. 1 SPFD 法 [5]
SPFD 法を用いて人体内誘導電界を計算する場合、2 段階の解析を行う。はじめに、電界ばく露時に人体表面に誘起される表面電荷を計算するため、人体外部の電界分布を計算する。次に、得られた表面電荷を境界条件として用い、SPFD 法で体内誘導電界を計算する。SPFD 法では、人体をボクセルでモデル化し、生体内のボクセルに含まれる全ての節点について、スカラーポテンシャルを未知数とする以下の方程式を立てる。
式 (1) は図１に示す節点0とその節点に隣接する6つの接点r(r =1∼6)のスカラーポテンシャルからなる方程式となる。
なお、Srは節点0−r間のアドミタンス、ψr は節点 r におけるスカラーポテンシャルを表している。すなわち、SPFD法は図２のような回路網にモデル化される。また、境界条件は外部電界計算によって求められた表面電荷密度分布ρs を用いて、電流連続の式により式 (2) と表される。
なお、ˆnはモデル表面に垂直な法線単位ベクトルJ は電流密度を意味する。この境界条件は図２の回路網では電流源として表される。
2. 2 SPFD 法の連立一次方程式
SPFD 法の連立一次方程式を行列形式Ａｘ＝ｂで表わせば、数値計算ライブラリを用いた計算が可能となる。なお、この時の係数行列Ａはアドミタンス、右辺ベクトル b は表面から流れ込む電流から成り立つ。ここで、 SPFD 法における係数行列 A の行ベクトルの和及び列ベクトルの和は 0 であり、rank-deficientであるため、非正則である。係数行列が非正則である連立一次方程式では、解 x が一意に決まらない。また、この時の連立一次方程式は右辺ベクトルが係数行列 A の列空間に存在するときのみ解を持ち、右辺ベクトルがこの条件を満たさない時、解は存在しない。SPFD 法の連立一次方程式では、右辺ベクトルの全要素の総和が0であれば、右辺ベクトルは列空間に存在し、解は存在するようになる。行ベクトルの和が0であることからも、右辺ベクトルの全要素の和が0でなければ、連立一次方程式Ax =bを満たす解が存在しないことがわかる。また、SPFD 法における右辺ベクトルの全要素の和は表面から流れ込む電流の総和であり、理論的には電流連続の式を満たし、この総和は0になるはずである。しかし、右辺ベクトルの値は外部電界計算によって求めるため、数値誤差や離散化による誤差などが含まれ、多くの場合、この総和は 0 にはならない。
このような連立一次方程式に対して、反復法を用いて解こうとすると解法によっては解が発散してしまうことがある。1 行1 列削除し、係数行列を正則にして解けば解は収束するが、電流の総和が 0 からかけ離れている場合は得られる解に大きな誤差が含まれる。
そのため、このような問題に対して、連立一次方程式 Ax = bを満たす解ではなく、残差ノルム ‖r‖L2 = ‖b Ax‖L2 を最小とするような解を妥協解として得る手法が考えられる。
SPFD 法の連立一次方程式を行列形式Ａｘ＝ｂで表わせば、数値計算ライブラリを用いた計算が可能となる。なお、この時の係数行列Ａはアドミタンス、右辺ベクトル b は表面から流れ込む電流から成り立つ。ここで、 SPFD 法における係数行列 A の行ベクトルの和及び列ベクトルの和は 0 であり、rank-deficientであるため、非正則である。係数行列が非正則である連立一次方程式では、解 x が一意に決まらない。また、この時の連立一次方程式は右辺ベクトルが係数行列 A の列空間に存在するときのみ解を持ち、右辺ベクトルがこの条件を満たさない時、解は存在しない。SPFD 法の連立一次方程式では、右辺ベクトルの全要素の総和が0であれば、右辺ベクトルは列空間に存在し、解は存在するようになる。行ベクトルの和が0であることからも、右辺ベクトルの全要素の和が0でなければ、連立一次方程式Ax =bを満たす解が存在しないことがわかる。また、SPFD 法における右辺ベクトルの全要素の和は表面から流れ込む電流の総和であり、理論的には電流連続の式を満たし、この総和は0になるはずである。しかし、右辺ベクトルの値は外部電界計算によって求めるため、数値誤差や離散化による誤差などが含まれ、多くの場合、この総和は 0 にはならない。
このような連立一次方程式に対して、反復法を用いて解こうとすると解法によっては解が発散してしまうことがある。1 行1 列削除し、係数行列を正則にして解けば解は収束するが、電流の総和が 0 からかけ離れている場合は得られる解に大きな誤差が含まれる。
そのため、このような問題に対して、連立一次方程式 Ax = bを満たす解ではなく、残差ノルム ‖r‖L2 = ‖b Ax‖L2 を最小とするような解を妥協解として得る手法が考えられる。
2. 3 残差ノルム最小の解を得る方法
残差ノルム最小の解を得る手法として、Moore-Penrose の擬似逆行列 [6] や最小化残差法 [7] などがあるが、ここでは、解法にかかわらず右辺ベクトルの補正をすることで残差ノルム最小の解を得る方法を提案する。
連立一次方程式 Ax = b の残差ノルムが最小である時の残差ベクトルを rmin とし、その時の解ベクトルを x∗とすると、これらは式 (3) を満たす。式 (3) を整理すると、式 (4) となる。
式 (4) は係数行列 A、右辺ベクトル b − rmin の連立一次方程式が残差ノルムを最小とする解 x∗を持つことを表しているので、これを解けば残差ノルム最小の解を得ることができる。元の連立一次方程式 Ax = b から補正後の右辺ベクトル˜b を得るのに必要な演算は式 (5) のみなので、様々な解法に適用可能であると予想される。
2. 4 残差ノルムの最小値及び残差ベクトル
SPFD 法における N 元連立一次方程式 Ax = b の残差ノルム ‖r‖L2 の最小値を考える。残差ノルムの定義は式 (6) であり、残差ベクトルの i 番目の要素 ri は右辺ベクトル b、解ベクトル x の i 番目の要素 bi 、 xi と係数行列 A の i 行 j 列の要素aij を用いて、式 (7) のように表される。
式 (6) の最小値を求めるには式 (8) の関数 f (r) の最小値を求めれば良い。
次に右辺ベクトルの全要素の和が式 (9) であるときの SPFD 法の連立一次方程式における残差ベクトル r が受ける制約条件を考える。
残差ノルム最小の解を得る手法として、Moore-Penrose の擬似逆行列 [6] や最小化残差法 [7] などがあるが、ここでは、解法にかかわらず右辺ベクトルの補正をすることで残差ノルム最小の解を得る方法を提案する。
連立一次方程式 Ax = b の残差ノルムが最小である時の残差ベクトルを rmin とし、その時の解ベクトルを x∗とすると、これらは式 (3) を満たす。式 (3) を整理すると、式 (4) となる。
式 (4) は係数行列 A、右辺ベクトル b − rmin の連立一次方程式が残差ノルムを最小とする解 x∗を持つことを表しているので、これを解けば残差ノルム最小の解を得ることができる。元の連立一次方程式 Ax = b から補正後の右辺ベクトル˜b を得るのに必要な演算は式 (5) のみなので、様々な解法に適用可能であると予想される。
2. 4 残差ノルムの最小値及び残差ベクトル
SPFD 法における N 元連立一次方程式 Ax = b の残差ノルム ‖r‖L2 の最小値を考える。残差ノルムの定義は式 (6) であり、残差ベクトルの i 番目の要素 ri は右辺ベクトル b、解ベクトル x の i 番目の要素 bi 、 xi と係数行列 A の i 行 j 列の要素aij を用いて、式 (7) のように表される。
式 (6) の最小値を求めるには式 (8) の関数 f (r) の最小値を求めれば良い。
次に右辺ベクトルの全要素の和が式 (9) であるときの SPFD 法の連立一次方程式における残差ベクトル r が受ける制約条件を考える。
ここで残差ベクトルの全要素の和は式 (10) のように表せる。
第 1 項は右辺ベクトルの要素和であり、今回は式 (9) としたので、Itotal となる。また、SPFD 法の係数行列は任意の行で式(11) を満たす。さらに、対称行列であるので、式 (12) も同時に満たす。
式 (9)、(12) から式 (10) は式 (13) となる。
次に、式 (13) の制約のもとで式 (8) の関数 f (r) の最小値及び、その時の残差ベクトル r をラグランジュの未定乗数法を用いて求める。このときのラグランジュ関数は式 (14) となる。
最小となる r を求めるためには式 (15)、(16) の連立方程式を解けばよい。
これを解くと、ri 、 λ はそれぞれ以下のように求まる。
以上より、残差ノルム最小となる残差ベクトル rmin、残差ノルム ‖rmin‖L2 はそれぞれ式 (19)、(20) となることが示された。
実際に計算を行う際は、右辺ベクトルの総和 Itotal と要素数 Nから残差ノルム最小となる残差ベクトル rmin を求め、右辺ベクトルを式 (5) のように補正すればよい。
第 1 項は右辺ベクトルの要素和であり、今回は式 (9) としたので、Itotal となる。また、SPFD 法の係数行列は任意の行で式(11) を満たす。さらに、対称行列であるので、式 (12) も同時に満たす。
式 (9)、(12) から式 (10) は式 (13) となる。
次に、式 (13) の制約のもとで式 (8) の関数 f (r) の最小値及び、その時の残差ベクトル r をラグランジュの未定乗数法を用いて求める。このときのラグランジュ関数は式 (14) となる。
最小となる r を求めるためには式 (15)、(16) の連立方程式を解けばよい。
これを解くと、ri 、 λ はそれぞれ以下のように求まる。
以上より、残差ノルム最小となる残差ベクトル rmin、残差ノルム ‖rmin‖L2 はそれぞれ式 (19)、(20) となることが示された。
実際に計算を行う際は、右辺ベクトルの総和 Itotal と要素数 Nから残差ノルム最小となる残差ベクトル rmin を求め、右辺ベクトルを式 (5) のように補正すればよい。
2. 5 補正方法の妥当性評価
前節で提案した補正方法を図 3 の回路モデルに用いて実際に計算を行い。その妥当性を評価した。この回路はアドミタンス S1 ∼ S6、電流源 Ia、 Ib、 Ic からなり、未知数を節点電位 V1 ∼ V 4 とする。それぞれのアドミタンスをS1 = 1 S/m、 S2 = 2 S/m、 S3 = 3 S/m、 S4 = 4 S/m、 S5 =5 S/m、 S6 = 6 S/m とし、それぞれの電流値は総和が 0 でなくなるように、Ia = 1 A、 Ib = −2 A、 Ic = −1 A とした。この時の連立一次方程式は式 (21) となる。
前節で提案した補正方法を図 3 の回路モデルに用いて実際に計算を行い。その妥当性を評価した。この回路はアドミタンス S1 ∼ S6、電流源 Ia、 Ib、 Ic からなり、未知数を節点電位 V1 ∼ V 4 とする。それぞれのアドミタンスをS1 = 1 S/m、 S2 = 2 S/m、 S3 = 3 S/m、 S4 = 4 S/m、 S5 =5 S/m、 S6 = 6 S/m とし、それぞれの電流値は総和が 0 でなくなるように、Ia = 1 A、 Ib = −2 A、 Ic = −1 A とした。この時の連立一次方程式は式 (21) となる。
このときの電流の総和は −2.0 であるので、式 (19)、(20) から残差ノルム最小となる残差ベクトルは rmin = −0.5[1 1 1 1]t、 残差ノルムは ‖rmin‖L2 = 1.0 となる。この連立一次方程式の右辺ベクトルに対して補正を行なった上で共役勾配法 (CG: ConjugateGradient Method) を用いて解く。また、 Moore-Penrose の擬似逆行列を用いて同様の連立一次方程式を解く。Moore-Penroseの擬似逆行列 A+と右辺ベクトルとの積 x = A+b によって得られる解 x は残差ノルム最小の解であるという特徴があるため、その時の残差ベクトル及び残差ノルムはそれぞれ、式 (19)、(20)を満たすことが予想される。この条件で計算を行った時、それぞれの手法によって得られる解 V1 ∼ V4、電界 E1 ∼ E6 を比較する。なお、電界値は節点間の電位差から求める。
表 1–3 に各手法で得られたスカラーポテンシャル、電界、残差ベクトル及び残差ノルムを示す。表 2、 3 から補正を用いて得られた電界、 残差ベクトル、 残差ノルムの値は Moore-Penroseの擬似逆行列を用いて得られた値とそれぞれ一致することが確認できた。したがって、補正により残差ノルム最小の解を求められることが確認された。
表 1–3 に各手法で得られたスカラーポテンシャル、電界、残差ベクトル及び残差ノルムを示す。表 2、 3 から補正を用いて得られた電界、 残差ベクトル、 残差ノルムの値は Moore-Penroseの擬似逆行列を用いて得られた値とそれぞれ一致することが確認できた。したがって、補正により残差ノルム最小の解を求められることが確認された。
3. 数値人体モデルを用いた解析
3. 1 ばく露条件
図 4 のように数値人体モデル前方に波源として微小ダイポールを配置した場合の内部誘導電界計算を提案の補正方法とともに行った。波源の周波数は 85 kHz とし、数値人体モデルには、情報通信研究機構によって開発された日本人成人男性モデルTARO モデル [8] を用いた。本モデルは空間分解能 2 mm で51 種類の組織で構成されている。外部電界計算時の TARO モデルの材質は完全導体とし、内部誘導電界計算時の人体組織の電気定数は文献 [9] より引用した。外部電界計算には AET 社の有限積分法に基づく電磁解析シミュレーションソフト CSTSTUDIO SUITE Low Frequency Solver を用いた。このときに得られた電流分布の総和及び最小残差ノルムなどの値を表 4に示す。
3. 2 残差ノルム最小の解の妥当性
補正によって求められた残差ノルム最小の解の妥当性を検討した。外部電界計算によって得られた電流値に対して、式 (22)を行い、電流の総和 Itotal = 0 である電流分布を模擬した。
この電流分布を用いて得られる誘導電界を正しい値であるとし、補正によって求められる誘導電界と比較する。
図 5 に SPFD 法によって求められた内部誘導電界強度分布を示す。また、内部誘導電界の統計量は表 5 のようになっており、補正によって得られた値は概ね一致していた。
3. 1 ばく露条件
図 4 のように数値人体モデル前方に波源として微小ダイポールを配置した場合の内部誘導電界計算を提案の補正方法とともに行った。波源の周波数は 85 kHz とし、数値人体モデルには、情報通信研究機構によって開発された日本人成人男性モデルTARO モデル [8] を用いた。本モデルは空間分解能 2 mm で51 種類の組織で構成されている。外部電界計算時の TARO モデルの材質は完全導体とし、内部誘導電界計算時の人体組織の電気定数は文献 [9] より引用した。外部電界計算には AET 社の有限積分法に基づく電磁解析シミュレーションソフト CSTSTUDIO SUITE Low Frequency Solver を用いた。このときに得られた電流分布の総和及び最小残差ノルムなどの値を表 4に示す。
3. 2 残差ノルム最小の解の妥当性
補正によって求められた残差ノルム最小の解の妥当性を検討した。外部電界計算によって得られた電流値に対して、式 (22)を行い、電流の総和 Itotal = 0 である電流分布を模擬した。
この電流分布を用いて得られる誘導電界を正しい値であるとし、補正によって求められる誘導電界と比較する。
図 5 に SPFD 法によって求められた内部誘導電界強度分布を示す。また、内部誘導電界の統計量は表 5 のようになっており、補正によって得られた値は概ね一致していた。
3. 3 様々な解法を用いた計算
数値計算ライブラリに含まれる連立一次方程式の様々な数値解法を用いて計算を行い、その時の計算時間を測定した。解法として、
• 共役勾配（CG: Conjugate Gradient）法
• 不完全コレスキー分解付き CG（ICCG: Incomplete Cholesky decomposition CG）法
• 不完全 LU 分解付き CG （ILUCG: Incomplete LU CG）法
• 最小残差（MINRES: Minimize Residual）法
• 対称 LQ（SYMMLQ: Symmetric LQ）法
• 代数的多重格子法付き CG （AMG-CG: Alegebraic Multi Grid CG）法
を用いる。CG 法、ICCG 法、ILUCG 法、MINRES 法、SYMMLQ 法は Mathworks 社の MATLAB 2016a の線形代数ライブラリを用いた。AMG-CG 法は代数的多重格子 (AMG)法 [10] を CG 法に適用した解法であり、みずほ情報総研によって開発された数値計算ライブラリを用いた。本ライブラリは、CPU による計算と GPU（Graphics Processor Unit）による並列計算の両方が可能であり、今回は CPU、GPU、両方の計算をそれぞれ行った。なお、CPU の計算では、計算の並列化は行わなかった。これらの計算を表 6 の計算環境で実行した。
収束条件は相対残差ノルム ‖b Ax‖L2/‖b‖L2 < 1.0 × 10−5とし、各解法が収束条件を満たすまでにかかった時間を計測する。なお、計算時間の計測は前処理を含む各数値解法の解を求めるルーチンのみとした。
また、これらの解法に加え比較対象として、従来広く用いられている逐次過緩和（SOR: Successive Over Relaxation）法の計算も行った。SOR 法は、異なる収束基準を用いていることや収束性が緩和係数の値に依存することなどから、単純に上記の解法と計算時間の比較をすることは出来ないが今回は計算の一例として、緩和係数を 1.5 とし、50000 反復にかかった計算時間を比較対象とする。なお、50000 反復目の相対残差ノルムは3.75 × 10−2であったため、‖b Ax‖L2/‖b‖L2 < 1.0 × 10−5上記の収束条件を満たすためにはさらに反復計算を行う必要がある。なお、SOR 法の計算は C 言語で記述し、bash 上で実行した。
それぞれの数値解法ごとの計算時間を図 6 に示す。AMG-CG法は一番計算時間が短く GPU による並列計算を行った場合は3.2 秒、並列化を行わなかった場合は 96.7 秒であった。次に計算時間が短かったのは前処理付き共役勾配法の ILUCG 法、ICCG 法の 2 つで両者とも計算時間は 1200 秒程度であった。これらの解法を SPFD 法の連立一次方程式のとしてを用いた場合、従来用いられてきた CG 法や SOR 法を用いた場合に比べて、より高速な解析が可能であることがわかった。
数値計算ライブラリに含まれる連立一次方程式の様々な数値解法を用いて計算を行い、その時の計算時間を測定した。解法として、
• 共役勾配（CG: Conjugate Gradient）法
• 不完全コレスキー分解付き CG（ICCG: Incomplete Cholesky decomposition CG）法
• 不完全 LU 分解付き CG （ILUCG: Incomplete LU CG）法
• 最小残差（MINRES: Minimize Residual）法
• 対称 LQ（SYMMLQ: Symmetric LQ）法
• 代数的多重格子法付き CG （AMG-CG: Alegebraic Multi Grid CG）法
を用いる。CG 法、ICCG 法、ILUCG 法、MINRES 法、SYMMLQ 法は Mathworks 社の MATLAB 2016a の線形代数ライブラリを用いた。AMG-CG 法は代数的多重格子 (AMG)法 [10] を CG 法に適用した解法であり、みずほ情報総研によって開発された数値計算ライブラリを用いた。本ライブラリは、CPU による計算と GPU（Graphics Processor Unit）による並列計算の両方が可能であり、今回は CPU、GPU、両方の計算をそれぞれ行った。なお、CPU の計算では、計算の並列化は行わなかった。これらの計算を表 6 の計算環境で実行した。
収束条件は相対残差ノルム ‖b Ax‖L2/‖b‖L2 < 1.0 × 10−5とし、各解法が収束条件を満たすまでにかかった時間を計測する。なお、計算時間の計測は前処理を含む各数値解法の解を求めるルーチンのみとした。
また、これらの解法に加え比較対象として、従来広く用いられている逐次過緩和（SOR: Successive Over Relaxation）法の計算も行った。SOR 法は、異なる収束基準を用いていることや収束性が緩和係数の値に依存することなどから、単純に上記の解法と計算時間の比較をすることは出来ないが今回は計算の一例として、緩和係数を 1.5 とし、50000 反復にかかった計算時間を比較対象とする。なお、50000 反復目の相対残差ノルムは3.75 × 10−2であったため、‖b Ax‖L2/‖b‖L2 < 1.0 × 10−5上記の収束条件を満たすためにはさらに反復計算を行う必要がある。なお、SOR 法の計算は C 言語で記述し、bash 上で実行した。
それぞれの数値解法ごとの計算時間を図 6 に示す。AMG-CG法は一番計算時間が短く GPU による並列計算を行った場合は3.2 秒、並列化を行わなかった場合は 96.7 秒であった。次に計算時間が短かったのは前処理付き共役勾配法の ILUCG 法、ICCG 法の 2 つで両者とも計算時間は 1200 秒程度であった。これらの解法を SPFD 法の連立一次方程式のとしてを用いた場合、従来用いられてきた CG 法や SOR 法を用いた場合に比べて、より高速な解析が可能であることがわかった。
4. 結 論
本研究では、右辺ベクトルの全要素の和が 0 でない矛盾したSPFD 法の連立一次方程式に対する補正方法の提案と有効性の検討を行なった。提案した右辺ベクトルの補正方法を用いれば、残差ノルム最小の解を得ることができ、妥当な結果が得られることを確認した。また、右辺ベクトルの要素和が大きい時、解法によっては大きな誤差が出たり、発散する場合があったが、補正を連立一次方程式を解く時の前処理として行うことで残差ノルム最小の解を数値計算ライブラリに含まれる様々な数値解法で得られることが確認できた。また、AMG-CG 法、ILUCG法、ICCG 法などの解法を用いれば、従来手法より高速な解析が可能であることがわかった。
本研究では、右辺ベクトルの全要素の和が 0 でない矛盾したSPFD 法の連立一次方程式に対する補正方法の提案と有効性の検討を行なった。提案した右辺ベクトルの補正方法を用いれば、残差ノルム最小の解を得ることができ、妥当な結果が得られることを確認した。また、右辺ベクトルの要素和が大きい時、解法によっては大きな誤差が出たり、発散する場合があったが、補正を連立一次方程式を解く時の前処理として行うことで残差ノルム最小の解を数値計算ライブラリに含まれる様々な数値解法で得られることが確認できた。また、AMG-CG 法、ILUCG法、ICCG 法などの解法を用いれば、従来手法より高速な解析が可能であることがわかった。
Claims (2)
- スカラーポテンシャルを未知数とする連立1次方程式に、外部環境データを算入して物体表面から物体内への所定の物の流入量を算出する流入量算出工程を具備する方法であって、
上記の所定の物は、自由境界を有する流動的に移動可能なものを意味し、
上記流入量算出工程において、上記連立1次方程式を補正する工程を具備し、この補正は上記連立1次方程式における各行から所定の値を差し引くことにより行われる
ことを特徴とする物体内への物流入量算出方法。
- 請求項１記載の方法であって、
上記の所定の物が電流である方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2018005758A JP2019124595A (ja) | 2018-01-17 | 2018-01-17 | 物体内への物流入量算出方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2018005758A JP2019124595A (ja) | 2018-01-17 | 2018-01-17 | 物体内への物流入量算出方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| JP2019124595A true JP2019124595A (ja) | 2019-07-25 |
Family
ID=67397881
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2018005758A Pending JP2019124595A (ja) | 2018-01-17 | 2018-01-17 | 物体内への物流入量算出方法 |
Country Status (1)
| Country | Link |
|---|---|
| JP (1) | JP2019124595A (ja) |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JPH0542116A (ja) * | 1991-08-13 | 1993-02-23 | Hitachi Ltd | 生体表面電位分布から生体内電流分布を推定する信号処理方法 |
| JPH0749857A (ja) * | 1993-06-03 | 1995-02-21 | Daikin Ind Ltd | 物理量解析方法およびその装置 |
| JPH07191965A (ja) * | 1993-12-27 | 1995-07-28 | Mitsubishi Electric Corp | 対象システムを最適設計または推定する方法およびその装置 |
| JP2008269329A (ja) * | 2007-04-20 | 2008-11-06 | Murata Mfg Co Ltd | 連立一次方程式の解を反復的に決定する方法 |
| US20100290675A1 (en) * | 2006-07-27 | 2010-11-18 | Alvin Wexler | High definition impedance imaging |
| JP2012247973A (ja) * | 2011-05-27 | 2012-12-13 | Denso Corp | 回路連携解析シミュレーション装置および回路連携解析方法 |
| WO2015122369A1 (ja) * | 2014-02-14 | 2015-08-20 | 国立大学法人東京大学 | 脳内電流シミュレーション方法とその装置，及び脳内電流シミュレーション装置を含む経頭蓋磁気刺激システム |
-
2018
- 2018-01-17 JP JP2018005758A patent/JP2019124595A/ja active Pending
Patent Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JPH0542116A (ja) * | 1991-08-13 | 1993-02-23 | Hitachi Ltd | 生体表面電位分布から生体内電流分布を推定する信号処理方法 |
| JPH0749857A (ja) * | 1993-06-03 | 1995-02-21 | Daikin Ind Ltd | 物理量解析方法およびその装置 |
| JPH07191965A (ja) * | 1993-12-27 | 1995-07-28 | Mitsubishi Electric Corp | 対象システムを最適設計または推定する方法およびその装置 |
| US20100290675A1 (en) * | 2006-07-27 | 2010-11-18 | Alvin Wexler | High definition impedance imaging |
| JP2008269329A (ja) * | 2007-04-20 | 2008-11-06 | Murata Mfg Co Ltd | 連立一次方程式の解を反復的に決定する方法 |
| JP2012247973A (ja) * | 2011-05-27 | 2012-12-13 | Denso Corp | 回路連携解析シミュレーション装置および回路連携解析方法 |
| WO2015122369A1 (ja) * | 2014-02-14 | 2015-08-20 | 国立大学法人東京大学 | 脳内電流シミュレーション方法とその装置，及び脳内電流シミュレーション装置を含む経頭蓋磁気刺激システム |
Non-Patent Citations (4)
| Title |
|---|
| "SPFD法における行列構造に関する検討", 電子情報通信学会大会講演論文集(CD-ROM), vol. 2016, JPN6022001928, 6 September 2016 (2016-09-06), ISSN: 0004688863 * |
| "SPFD法による電界ばく露時の人体内誘導電界計算に関する検討", 電子情報通信学会技術研究報告, vol. 117, no. 384, JPN6022001924, 11 January 2018 (2018-01-11), pages 7 - 12, ISSN: 0004827042 * |
| "ボクセルデータ用高速多重極表面電荷法による低周波磁界誘導電界計算", ＩＥＥＪ ＴＲＡＮＳ．ＦＭ, vol. 126, no. 5, JPN6022001926, 2006, ISSN: 0004688862 * |
| "接触電流を模擬した生体モデル内における誘導電界の計測", 宮崎大学工学部紀要, JPN6022001929, 31 July 2014 (2014-07-31), pages 43 - 48, ISSN: 0004688864 * |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US10928433B2 (en) | Method and program for calculating potential, current, and peripheral electromagnetic field in electric circuit | |
| Soleimani | Computational aspects of low frequency electrical and electromagnetic tomography: A review study | |
| Anees et al. | Time domain finite element method for Maxwell’s equations | |
| Vipiana et al. | Optimized numerical evaluation of singular and near-singular potential integrals involving junction basis functions | |
| Rizea et al. | Finite difference approach for the two-dimensional Schrödinger equation with application to scission-neutron emission | |
| Bera et al. | A MATLAB-based boundary data simulator for studying the resistivity reconstruction using neighbouring current pattern | |
| Sobczyk et al. | Spectral density reconstruction with Chebyshev polynomials | |
| Zhang et al. | A numerical computation forward problem model of electrical impedance tomography based on generalized finite element method | |
| Seoane et al. | An accurate and efficient three‐dimensional high‐order finite element methodology for the simulation of magneto‐mechanical coupling in MRI scanners | |
| Meunier et al. | Unstructured–PEEC method for thin electromagnetic media | |
| US20090177456A1 (en) | Mixed Decoupled Electromagnetic Circuit Solver | |
| JP2019124595A (ja) | 物体内への物流入量算出方法 | |
| Vancea | Entanglement entropy in the σ-model with the de Sitter target space | |
| Beddek et al. | Solution of Large Stochastic Finite Element Problems—Application to ECT-NDT | |
| Jin et al. | Numerical estimation of piecewise constant Robin coefficient | |
| Zhao et al. | A novel formulation with coulomb gauge for 3-D magnetostatic problems using edge elements | |
| Akçelik et al. | Shape determination for deformed electromagnetic cavities | |
| Dłotko et al. | Topoprocessor: An efficient computational topology toolbox for h-oriented eddy current formulations | |
| Sorti et al. | Data-driven simulation of transient fields in air–coil magnets for accelerators | |
| Nomura et al. | Analysis of Current Density Inside Human Body Using Geometric Multi-Grid Method | |
| Du et al. | Interval inverse analysis of hyperbolic heat conduction problem | |
| Laakso et al. | Improving the computational speed and reducing the staircasing error for simulations of human exposure to low frequency magnetic fields | |
| Nomura et al. | Contact Current Density Analysis Inside Human Body in Low-Frequency Band Using Geometric Multi-Grid Solver | |
| Faugeras | Tokamak plasma boundary reconstruction using toroidal harmonics and an optimal control method | |
| Aiello et al. | A GMRES iterative solution of FEM‐BEM global systems in skin effect problems |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20210114 |
|
| RD03 | Notification of appointment of power of attorney |
Free format text: JAPANESE INTERMEDIATE CODE: A7423 Effective date: 20210114 |
|
| A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A821 Effective date: 20210114 |
|
| A977 | Report on retrieval |
Free format text: JAPANESE INTERMEDIATE CODE: A971007 Effective date: 20211222 |
|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20220125 |
|
| A02 | Decision of refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A02 Effective date: 20220719 |