（様式 20）成果報告書

TSUBAME 共同利用 令和 4 年度 学術利用 成果報告書

利用課題名 GPU クラスタを用いたミリ波帯大規模広帯域電波伝搬シミュレーション
英文：Large-Scale Propagation Simulations of Ultra-Wideband Electromagnetic Fields in Millimeters and Terahertz
Wave Frequencies Using GPU Cluster
チャカロタイ ジェドヴィスノプ
Jerdvisanop Chakarothai
国立研究開発法人情報通信研究機構
National Institute of Information and Communications Technology
https://www.nict.go.jp
ミリ波帯・テラヘルツ帯における電波は他の周波数帯と比較して、電波の直進性が大きく、見通し伝搬の形態で利用されることが多い。しかし、ひとたび電波伝搬経路のフレネルゾーン内に吸収体・散乱体・反射体などの障害物が存在すると、受信電力の予測手法の複雑さが増大し、不確かさが大きくなる。特に屋内伝搬環境は、壁、天井、机、棚や椅子等の様々な障害物が存在し、伝搬状況を正しく予測することは困難である。そこで、本研究では、オフィス等の室内環境においてミリ波帯（特に 28 GHz 帯）電波の受信電力を高精度に予測するために、大規模電磁界シミュレーション技術を開発した。特に実際の構造物の離散化及び機能的な電波散乱シート（入射された電波を様々な方向へ散乱させるための構造）を従来のプログラムに組み込むために、改良・拡張した。(1)屋内電波環境の数値解析のために、床・壁のみではなく、より実際のオフィスの構造に近い解析モデルにするために、3 次元構造モデル（テーブル、椅子、パソコンなど）を解析空間の解像度に合わせて離散化して組み込むことを可能にした。 (2)電波散乱シートを解析空間に組み込むために、完全磁気導体及び完全導体を同一のメッシュで取り扱えるようにするために、
FDTD 法の定式化を行い、三次元解析による定式化の妥当性検証を行った。実際に散乱板を解析モデルに組み込み、FDTD 法によって解析した結果とモーメント法によるものと比較した。今後、散乱板を組み込んだ大規模電磁界解析を行う予定である。
Millimeter-wave and terahertz electromagnetic (EM) waves are often used in line-of-sight propagation. However, if
there are absorbers, scatterers, reflectors, or other obstacles within the Fresnel zone of the wave propagation path, theuncertainty of the received power prediction method increases significantly. Especially in indoor propagation
environments, there are various obstacles such as walls, ceilings, desks, shelves, and chairs, making accurate predictionof propagation conditions difficult. Therefore, in this study, we developed large-scale EM simulation techniques for
accurate prediction of the received power of millimeter-wave (especially 28 GHz band) waves in indoor environments.
In particular, we improved and extended conventional programs to discretize actual structures and incorporate
functional structures such as EM scattering sheet. (1) To perform numerical analysis of indoor radio wave environments,
in addition to floors and walls, we discretize 3D structural models (tables, chairs, computers, etc.) to match resolutionof the analysis space. (2) To incorporate EM scattering sheets into the analysis region, we incorporate a surface
impedance into the FDTD method to handle perfect magnetic conductors and perfect conductors on the same surface
plane and demonstrated the validity of the formulation through 3D numerical analysis. We compared the results by theFDTD analysis which includes a scattering sheet with those obtained by the moment method. In the future, we plan toperform large-scale electromagnetic field analysis incorporating scattering sheets.
Keywords: millimeter-wave, FDTD, indoor propagation, large-scale simulation, electromagnetic scattering sheet

背景と目的

内電波環境においては、様々な障害物によって電波が5G 無線通信や Beyond 5G の次世代無線通信技術

散乱され、複雑な電波環境を形成するため、受信電力は主にミリ波帯電波が利用される。しかしながら、ミリ波を正確に予測することが難しい。そこで、本研究では、帯電波の伝搬経路のフレネルゾーン内に吸収体・散乱様々な伝搬環境(特に室内環境)におけるミリ波帯電波体・反射体などの障害物が存在すると、受信電力の予の受信電力を高精度に予測するために、超大規模電測手法の複雑さが増大し、不確かさが大きくなる。また、磁界シミュレーション技術を新たに開発すると共に、機波長が非常に短いため、物体の表面形状等による散能的なの構造物(リフレクトアレーや電波散乱シート等)乱が生じ、従来のレイトレーシング法や光学近似法にに対する数値解析法を提案し、その妥当性を示す．

よる高精度な解析が困難である。特にオフィス等の屋（様式 20）成果報告書有する．

概要
本研究では、オフィス等の室内環境においてミリ波帯

まず，一次元 FDTD 法における媒質境界周辺での電（特に 5G 無線通信に使われる 28 GHz 帯）電波の受信電力を高精度に予測するために、TSUBAME3.0 GPU界および磁界の配置を示す．電界の接線成分𝐸𝐸𝑦𝑦𝑡𝑡 を次クラスタを用いた大規模並列化電磁界シミュレーションによって表すことができる[11]．

式のように複素周波数領域での表面インピーダンス𝑍𝑍𝑠𝑠技術を開発した。解析に用いた計算手法は、電磁界解

𝑍𝑍𝑠𝑠 (𝑠𝑠)� 𝑠𝑠𝐻𝐻𝑧𝑧𝑡𝑡 (𝑠𝑠)𝐸𝐸𝑦𝑦𝑡𝑡 (𝑠𝑠) = �𝑠𝑠析 で よ く 使 わ れ て い る 時 間 領 域 有 限 差 分 (Finitedifference time-domain, FDTD)法である[1]。
本研究に用いたプログラムは文献[2]-[7]にも用いら(1)

ただし，𝐻𝐻𝑧𝑧𝑡𝑡 は境界における磁界の接線成分である．式(1)を時間領域で表すと，れた人体ばく露評価のためのプログラムを改良・拡張し𝐸𝐸𝑦𝑦𝑡𝑡 (𝑡𝑡) = 𝑋𝑋(𝑡𝑡) ∗たものである。昨年度に引き続き、物体形状を忠実に

𝜕𝜕𝐻𝐻𝑧𝑧𝑡𝑡 (𝑡𝑡)
= 𝑋𝑋(𝑡𝑡) ∗ 𝐻𝐻𝑧𝑧′ (𝑡𝑡)𝜕𝜕𝜕𝜕(2)再現できるようにするために、きめ細かくセルの各辺にとなり，X(t) = L-1{X(s)},細かく設定することができるように変更した他、解析モで表される．ここで，Zs(s)は以下のような式で与えられデルのデータをハードディスクに保存せずに、直接プロる．

グラムの中に組み込めるようにした。加えて、新たな機X(s) = Zs(s)/s, 𝐻𝐻𝑧𝑧′ = ∂𝐻𝐻𝑧𝑧𝑡𝑡 (t)/∂t𝑍𝑍𝑠𝑠 (𝑠𝑠) = �

能的な構造物（リフレクトアレー）を解析空間に組み込

𝑠𝑠𝜇𝜇0 𝜇𝜇𝑟𝑟
𝑠𝑠𝜀𝜀0 𝜀𝜀𝑟𝑟 + 𝜎𝜎(3)むために、表面インピーダンスに対して、FDTD 法の定σ, εr, µr, ε0, µ0 はそれぞれ媒質の導電率，比誘電率，比式化を行い、提案する手法の妥当性を検討した。

透磁率，真空の誘電率および透磁率である．次のステ

本研究の解析対象としては、屋内での電波環境を可ップでは，X(s)を FILT によって時間領域波形へ変換し視化することを目的とした。実際にテーブルや椅子等た後，Prony 法を適用すると，z 領域において X(z)を次のオフィス器具を含めた解析モデルは、大きさ 3×2.4式のように表すことができる[2],[3]．
𝑁𝑁𝑙𝑙×2.6 m であり、必要な計算資源は、約 193 億セルの3𝑋𝑋(𝑧𝑧) = �

約 720GB の総使用メモリを要した。大規模 GPU クラス

𝑙𝑙=1𝐴𝐴𝑙𝑙,
1 − 𝑝𝑝𝑙𝑙 𝑧𝑧 −1(4)タにより高速化を行った結果、16 ノード(64 GPUs)によるNl は Prony 法によって得られた実数の極に対応してお計算時間は 1000 ステップ当たり、約 32 分で終えること

り，Al および pl は文献[12],[13]の手順に従って Prony 法ができ、実用的な時間で解析が可能であることを確認

から抽出された係数である．なお，通常 Prony 法を適用した。

すると，複素数の極も得られるが，式(3)で表される表
面インピーダンスは時間領域において単調的に減衰し表面インピーダンスの定式化
本研究において，高伝導性媒質のための表面インピーダンスが適用される．表面インピーダンス法に関しては，これまで様々な定式化手法が提案されてきたがていく関数で求まるため，実数の極のみで表現することができる．式(4)を z 領域での式(2)に代入すると，
𝑁𝑁𝑙𝑙

𝑁𝑁𝑙𝑙

𝑙𝑙=1

𝑙𝑙=1𝑡𝑡𝐸𝐸𝑦𝑦𝑡𝑡 (𝑧𝑧) = � 𝐸𝐸𝑦𝑦,𝑙𝑙(𝑧𝑧) = �𝐴𝐴𝑙𝑙
𝐻𝐻 ′ (𝑧𝑧)1 − 𝑝𝑝𝑙𝑙 𝑧𝑧 −1 𝑧𝑧(5)[8],[9]，本提案手法では，PEC および PMC を同一境界が得られる．式(5)より各実数の極に対応した電界の離面上で定義できることが特徴である．また，従来手法で散的な時間 n ステップ目における更新式は以下のようは，時間領域での厳密解や周波数領域でのベクトルフに求まる．

ィッティング(Vector Fitting, VF)法[10]を用いて定式化が行われたが，本手法は周波数領域での表面インピータ

𝑡𝑡,𝑛𝑛𝑡𝑡,𝑛𝑛−1𝐸𝐸𝑦𝑦,𝑙𝑙= 𝑝𝑝𝑙𝑙 𝐸𝐸𝑦𝑦,𝑙𝑙+11𝐴𝐴𝑙𝑙𝑡𝑡,𝑛𝑛+𝑡𝑡,𝑛𝑛−
�𝐻𝐻𝑧𝑧 2 − 𝐻𝐻𝑧𝑧 2 �∆𝑡𝑡

(6)

Prony 法を用いて時間領域でフィッティングを行うことで，𝑡𝑡,𝑛𝑛の総和で計算できる．
従って，電界の接線成分は𝐸𝐸𝑦𝑦,𝑙𝑙

高精度な解を与えながらも定式化が簡素になる利点を以下のような更新方程式が導出される．

ンスを FILT によって時間領域波形に変換してから一方では，磁界の接線成分は前進差分を適用すると，

（様式 20）成果報告書
1𝑡𝑡,𝑛𝑛+1𝑡𝑡,𝑛𝑛−𝐻𝐻𝑧𝑧 (𝑖𝑖+21) = 𝐻𝐻𝑧𝑧 (𝑖𝑖+21) −222∆𝑡𝑡
�𝐸𝐸 𝑡𝑡,𝑛𝑛 1 − 𝐸𝐸𝑦𝑦𝑛𝑛 (𝑖𝑖) �𝜇𝜇0 ∆𝑥𝑥 𝑦𝑦 (𝑖𝑖+2)(7)600400

最終的に，式(5)の𝐸𝐸𝑦𝑦𝑡𝑡,𝑛𝑛(𝑖𝑖+1) を式(7)の中に代入すれば，
2371

媒質境界面における磁界の更新方程式が得られる．
1𝑡𝑡,𝑛𝑛+1𝑡𝑡,𝑛𝑛−𝐻𝐻𝑧𝑧 (𝑖𝑖+21) = 𝐻𝐻𝑧𝑧 (𝑖𝑖+21) −222∆𝑡𝑡

𝑁𝑁𝑙𝑙𝐴𝐴𝑙𝑙𝜇𝜇0 ∆𝑥𝑥 + 2 ∑𝑙𝑙=1

𝑁𝑁𝑙𝑙𝑡𝑡,𝑛𝑛−1
�𝐸𝐸𝑦𝑦𝑛𝑛 (𝑖𝑖) − � 𝑝𝑝𝑙𝑙 𝐸𝐸𝑦𝑦,𝑙𝑙
�

(8)

𝑙𝑙=1

これによって従来の FDTD 法の磁界の配置において，完全電気導体を表面インピーダンスによって表すことができる．本提案法の計算フローチャートを図 2 に示す．

単位: mm

図 1 三次元モデルのメッシュデータからボクセル化（離散化）して得られた結果
テーブル

なお，本提案法は，高導電性媒質に限らずミリ波帯における生体組織へ適用することが可能であることが示されている[13]．一方，PMC の場合は，単純に境界面における磁界の接線成分をゼロとすることで実現でき椅子る．従って，本提案法を用いて，同一面内で PEC または PMC をモデルすることが可能である．
結果および考察
まず、図 1 のようなテーブルおよび椅子を 3 × 2.4 ×2.6 m の模擬部屋内に配置して解析を行った。扉の材
3

図 2 模擬オフィス部屋内にテーブルおよび椅子を配置したときのミリ波帯(28 GHz)屋内電磁界分布

質は木材とし、比誘電率および導電率は、それぞれ 2.0

次に、機能的な構造を有する電波散乱シートに対し及び 0.001 S/m とした。施工ボードはεr = 2.4, σ = 0.001て解析を行う。図 3 に異なる高さを有する 2 種類の金属

S/m、床面は、εr = 13, σ = 0.01 S/m の厚さ 10 cm のコン板で構成されている EMSS を示す。基準面における反クリートとした。解析モデルの解像度は、1 mm で、8 層射位相は、PMC が 0o で，PEC が-180o となっている。0oの完全整合層 (Perfectly matched layers)を含めた解析の反射位相は、金属板の高さをλ/4 分だけ変えることで領域は、3,024 × 2,424 × 2,624 セルである。解析周実現される。また、図 3 には、実際の EMSS の構造(モ波数は 28 GHz である。また、入射はモーメント法によっーメント法における解析モデル)と FDTD シミュレーショて求めたホーンアンテナの開口分布を波源とした。ホーンで使用される解析モデルを示している。解析周波数ンアンテナを用いた場合、開口面分布は磁流のみで表を 30 GHz と し 、 EMSS の 全 体 サ イ ズ は 20×20現できるが、任意波源の場合は電流及び磁流の両方cm2(20λ×20λ)であり、金属板の幅は 2 cm (2λである．

をモデリングする必要がある。図 2 には、テーブルや椅基準面での反射位相が-180o の金属板は、導電率σ =子を含む模擬部屋の外に波源（K バンドの標準ホーン

58×1010 S/m を有する金属の表面インピーダンスで置きアンテナ）を設置したときの電界分布を示す。図 2 より換えられ、反射位相がゼロの板は完全磁気導体(PMC)

28 GHz 帯電波はスタッドの間を透過することが分かり、で置き換えられている（磁界の接線成分がゼロ）。解析

ホーンアンテナからは鋭いビームが形成されていること

領域のサイズは 241×241×61，セルサイズは∆x = ∆y =を確認できた。結果から、波源が 1.5 m の高さに配置さ∆z = 1 mm、時間ステップ間隔は 1.67 ps とした。入射波れているため、高さ約 37 cm のテーブルおよび椅子か形は連続正弦波とし、20λ×20λ の大きさを有する EMSSらの散乱は少ないことが確認された。

に対する散乱断面積(Bistatic Scattering Cross Section,BSCS)を導出するための計算時間はわずか 23 秒であった。図 6 は垂直入射時(θinc = 0o)の表面インピーダン（様式 20）成果報告書ある。

ExHy

2 cm

謝辞 本研究を進めるあたり，表面インピーダンス法にk関して検討して頂いた法政大学大学院修了生の竹谷和真氏に感謝の意を表す．
20 cm
z
20 cm

PEC

Hyx

Exk

参考文献

yPMC

[1] K. S. Yee, “Numerical solution of initial boundary
value problems involving Maxwell’s equations inisotropic media,” IEEE Trans. Antennas Propag., vol. AP14, no. 3, pp. 302–307, May 1966.
[2] J. Chakarothai, K. Wake, S. Watanabe, “Scalable GPU-

図 3 電波散乱シート(EMSS)の実際の解析モデルおよび表面インピーダンスを用いた場合の FDTD 解析モデル[14]= ∆y = ∆z = 1 mm
= 1.67 ps
FDTD analysis

0oEinc

MoMFDTDθ

parallelized FDTD method for analysis of large-scale
electromagnetic

dosimetryproblems,”

AppliedComputational Electromagnetics (ACES) Journal, vol. 31,
no. 6, pp. 661-668, June 2016.

[3] J. Chakarothai, “Novel FDTD scheme for analysis
of frequency-dependent medium using fast inverse
Laplace transform and Prony’s method,” IEEE Trans.
Antennas Propagat., vol. 67, no. 9, pp. 6076-6089, Sep.

90o
-20 -10 0 10 20 [dBsm]
図 4 垂直入射時の電波散乱シートの散乱特性(実線：MoM, 破線:FDTD)[14]90oスを使った FDTD 解析で得られた BSCS およびモーメント法による電磁界解析で得られた BSCS を示している。
両者は非常によく一致しており、提案手法の妥当性が確認された。

まとめ、今後の課題
本研究では、28 GHz における室内電波伝搬環境を高精度に解析できる時間領域有限差分法を開発した。

2019.
[4] J. Chakarothai, S. Watanabe, K. Wake, “Numerical
dosimetry of electromagnetic pulse exposures using
FDTD method,” IEEE Trans. Antennas Propagat., vol.
66, no. 10, pp. 5397-5408, Oct. 2018.

[5] チャカロタイジェドヴィスノプ, 和氣加奈子,

渡辺聡一, 陳強, 澤谷邦男, “超広帯域電磁界解析のための周波数依存性 FDTD 法,” 電子情報通信学会 和文論文 C, vol. J102-C, no. 5, pp. 102-113,
May 2019.
[6] J. Chakarothai, K. Fujii, “A unified approach fortreatment of frequency-dependent materials in FDTD
method,” ISAP 2019, Xi’an, China, Oct. 30, 2019.

本手法では、これまでのプログラムに対して、(1)より実[7]際のオフィスの構造に近い解析モデルにするために、3“Convergence of a single-frequency FDTD solution in次元構造モデル（テーブル、椅子、パソコンなど）を解析

numerical dosimetry,” IEEE Trans. Microw. Theory

空間の解像度に合わせて離散化して組み込むことを可Techn., vol. 64, no. 3, pp. 707-714, March 2016.

能にした。(2) 機能的な構造体（電波散乱シート等）を解析空間に組み込むために、FDTD 法の定式化を行い、三次元解析による定式化の妥当性検証を行った。モーメント法による解析結果と比較して、同等な解析精度が得られることが示された。
今後は、電波散乱シートやリフレクトアレー等を実際に組み込んだ大規模電波伝搬解析へ拡張する予定で

J.

Chakarothai,

K.

Wake,

S.

Watanabe,[8] J. G. Maloney and G. S. Smith, “The use of surface
impedance concepts in the finite-difference timedomain method,” IEEE Trans. Antennas Propag., vol.
40, no. 1, pp. 38-48, 1992.
[9] J. H. Beggs, R. J. Luebbers, K. S. Yee, and K. S.
Kunz, “Finite-difference time-domain implementationof surface impedance boundary conditions,” IEEETrans. Antennas Propag., vol. 40, pp. 49-56, 1992.

（様式 20）成果報告書[10] B. Gustavsen and A. Semlyen, “Rationalapproximation of frequency domain responses byvector fitting,” IEEE Trans. Power Delivery, vol. 14, no.
3, pp. 1052-1061, July 1999.
[11] M. A. Leontovich, “On the approximate boundary
conditions for electromagnetic fields on the surface ofwellconducting

bodies,”inInvest.

Propag.

Radiowaves, B. A. Vedensky, Ed. Moscow: Academyof Science USSR, pp. 5-20, 1948.
[12] J. Chakarothai, K. Fujii, Y. Suzuki, J. Shibayama,and K. Wake, “Analyses of transient energy depositionin biological bodies exposed to electromagnetic pulses
using parameter extraction method,” IEICE Trans.
Commun., vol. E105-B, no. 6, pp. 694-706, June 2022.
[13] K. Takeya, J. Chakarothai, J. Shibayama, Y.
Suzuki, and K. Fujii, “Novel surface impedance
formulation via FILT and Prony method for FDTD
analyses of lossy media,” IEEE Antennas and WirelessPropagation Letters, vol. 23, no. 5, pp. 1593-1597, May2024.
[14] J. Chakarothai, K. Takeya, J. Shibayama, Y.
Suzuki, K. Fujii, “FDTD Analysis of Electromagnetic
Scattering Sheet Using Surface Impedance Method,”2024 IEEE International Workshop on Antenna
Technology (iWAT2024), April 15-18 2024 (Sendai,Japan).