— 48 —SENAC Vol. 48, No. 4（2015. 10）［研究成果］

ベクトルスーパーコンピュータ SX-ACE
による暑熱環境下体温上昇の高速解析
西尾 渉
浅野 陽平 平田 晃正
名古屋工業大学 大学院工学研究科

あらまし 熱ストレスの人体へのリスク評価には，体温上昇，発汗の時間変化を把握することが重要となり，それには熱解析と温熱生理反応を融合した定式化を用いる必要がある．しかしながら，これまで名古屋工業大学の研究グループが開発してきた体温上昇解析コードを用いて， 2時間の暑熱環境におけるばく露状況を解析するには，通常のワークステーションを用いた場合 7時間程度必要となる．上記のばく露解析を 20 秒程度で完了することが可能となれば，個体差によるばらつきなど様々な環境を模擬した解析を即時的に処理，普及啓発などへの応用展開が期待できる．本稿では，ベクトルスーパーコンピュータ SX-ACE に上記解析コードを実装し，解析時間の高速化について検討するとともに，解析例を示す．

1. まえがき
暑熱環境などの熱ストレスに対する人体安全性を評価するには，体内温度上昇，発汗の時間変化を解析，検討することが重要となる．名古屋工業大学の研究グループでは，部位ごとの血液温度差を考慮し，かつ温熱生理を考慮した混成熱解析手法を開発してきた[1]．この手法を精巧な数値人体モデルに適用し，暑熱環境に配置した場合の体内温度上昇，発汗量の解析を行い，測定値との比較によりその有効性を確認している[2]．しかしながら，物理計算に加え，システムバイオ
ロジー，具体的には発汗など体温調節反応を加味した複合的な解析を実施する場合，ベクトル化には適しておらず解析負荷も小さくない．例えば，2 時間の暑熱環境におけるばく露状況を解析するには，通常のワークステーション（CPU ：Intel® Xeon® W5590 @3.33GHz , 4 コア×2）において 7 時間程度必要である．たとえば，このばく露状況の解析を 20 秒程度で実行可能となれば，現実的なばく露環境下において，年齢や性別など複数の人体モデル，さらには熱調整機能の相違などの個人差を考慮した，様々なパターンのデータを即時的に処理することができる．得られたデータを総合的に評価することで，熱ストレスに関する準リアルタイムでのリスク評価が可能となる．本稿では，4 コアで構成されるノードを多数接続したベクトルスーパーコンピュータ
SX-ACE[3]へのコードの実装，および並列実行環境による計算時間の加速率について検討した結果を紹介するとともに，解析例を紹介する．

2. 解析手法
2.1 解剖学的数値人体モデル
解 析 に は ， 情 報 通 信 研 究 機 構 (NICT: National Institute of Information and Communication
Technology )開発の成人男性に相当する日本人の解剖学的数値人体モデルを使用した[4]．これらのモデルは皮膚，筋肉，脂肪，骨，骨髄液，心臓，レンズ，血管など 51 種類の組織から構成されており， 2mm の分解能を有している．また，ハイリスク群である 3 歳児相当のモデルについては[4]に基づき開発された[5]を用いた．
本研究では従来の血液温度計算を改良し，部位ごとに血液温度を考慮出来るようにするため，それぞれのモデルを部位ごとに分割して用いることにした．なお，頭部と胴体部における血流量ベクトルスーパーコンピュータ SX-ACE による暑熱環境下体温上昇の高速解析— 49 —

は十分に大きいため，2 つの部位の血液温度は同一なものと仮定して扱う．数値人体モデルの概観と部位分割の定義を図 1 の(a)，(b)にそれぞれ示す．

(a)

(b)

図 1. 数値人体モデルの外観(a)と部位分割(b)の定義図

2.2 温度上昇解析手法
2.2.1 生体熱輸送方程式
生体熱輸送方程式を計算することにより，温度上昇の数値計算を行う．一般に，生体内の温度上昇計算は，組織内および組織間の熱伝導による熱移動，血流による熱冷却，皮膚表面から外気への熱伝達を考慮することで導出され，体内座標を表す位置ベクトルを r とすると，次式で表される[6]．
T (r, t )C(r )   (r )   (K (r )T )  A(r, t )  B(r, t )  (T (r, t )  Tb,m (t ))(1)t
ここで，C は組織の比熱[J/kg·oC]， は密度[kg/m3]，T は組織の温度[oC]，K は熱伝導率[W/m·oC]，
A は代謝熱[W/m3]， B は血流定数[W/m3·oC]， Tb は血液温度[oC]である．また，境界条件はモデルの体表面および肺内部に対して，次式が適用される．
T (r , t ) K (r ) H (r )  (Ts (r , t )  Te (t ))  40.6  ( SW (r , Ts (r , t ))  PI ) / Sn(2)ここで， H は熱伝達率[W/m·oC]， TS は外気に接する組織の温度[oC]， Ta は外気の温度[oC]， SW は汗に関する係数[W/m2]，PI は不感蒸散量，n はその表面の法線方向の座標，S は人体の表面積[m2]である．また，40.6 は変換係数である．各種熱定数は，文献[1]に示したものを用いた．式(1)中の血液温度 Tb,m に関しては文献[7]を拡張し，Tb ,m (t )  Tb 0  

Qbtot (t )  Qb ,n (t )C b   b  Vbdt

(3)

Qbtot (t )   B (t )(Tb (t )  T ( r , t )) dVにより部位間の血液による熱交換を模擬し，部位ごとに血液温度の時間変化を計算することとした．ここで，m(=1,…,5)は分割した部位，n(=1,…,4)は部位間の隣接部分を表している．また，Tb0は血液の初期温度[oC]，Qbtot は当該部位で血液が獲得する熱量，Cb は血液の比熱[W/ oC]，  b は血液の密度[kg/m3]，Vb は血液の体積[m3]である．分割した部位間で血液により交換される熱量は文— 50 —

SENAC Vol. 48, No. 4（2015. 10）

献[8]より，Qb,n (t )   hx (Tb _ trunk (t )  Tb_limb(t ))dV(4)で与えられる．ここで，hx は熱交換係数[W/ oC]である．なお，前述の通り，頭部と胴体は血流量が多いため両部位の血液温度は同一であると仮定した[8]．
2.2.2 発汗のモデル化
汗に関する係数 SW は，文献[8]に基づき，SW ( r , t )   (r )[{11 tanh( 11ΔTS  10 )  10 }ΔTS
 { 21 tanh(  21ΔTH   20 )   20 }ΔTH ](5)で与える．式(4)の  TS ，TH は，それぞれ皮膚と視床下部の平均温度上昇である． および  は，
文献[7][8]で示したように，10 =1.20，11 =0.80，10 =0.19，11 =0.59， 20 =6.30， 21 =5.70， 20 =1.03，
 21 =1.98 とした．また，χ は，部位に依存する発汗係数である．
一般に，皮膚表面の蒸発熱量 E は次式で表される[8]．

E  2 .2  hc wf pcl ( Ps   A PA )(6)

ここで，hc は対流性熱伝達係数，w は皮膚の濡れ率，fpcl は被服の透湿率，PS,PA はそれぞれ皮膚表面および外気の飽和水蒸気圧[mmHg]， A は相対湿度である．皮膚の濡れ率 w は，“皮膚全体が汗で均一に覆われて，すべて蒸発すると仮定した場合の蒸発熱量（w=1 における E）”に対する“実際の人体表面からの蒸発熱量（式(4) に基づく発汗量）”の比で定義されるものである[8]．発汗が皮膚全面を完全に濡らすほどに達した場合は，それ以上の汗は熱冷却には有効に作用しない．このことから本解析では，式 (5) において，w=1 の値をとるときの蒸発熱量を，熱冷却として働く発汗量の上限とし，解析に組み込み，実装した．
高齢者のモデル化において，脚部においてのみ[9][10]（一部の文献では背中[10]も含む），最大発汗率の減少が報告されている．そこで，高齢者における発汗率のモデル化に際し，式(4)に対して，以下の変更を加えることとした． i)脚部の発汗率の減少を考慮するために，係数  r  を導入
する．ii)加齢に伴い，熱感度が低下するとの報告に基づき，発汗反応を誘発する閾値として  TS,decおよび  TH,dec を導入した．これらを加味した上で得られる式は以下の通りとなる．

SW (r , t )   (r )  ( r )[{11 tanh( 10 )  10 }ΔTS
 { 21 tanh(  21ΔTH   20 )   20 }ΔTH ](7)

ここで，TS  TS ,0  ΔTS ,dec0
ΔTS  TS  (TS ,0  ΔTS ,dec ),TS  TS ,0  ΔTS ,dec0
ΔTH  TH  (TH , 0  ΔTH ,dec ),TH  TH , 0  ΔTH ,dec
TH  TH ,0  ΔTH ,dec

(8)

(9)

である．なお，ここでは示さないが，上記の発汗率のモデル化を導入したうえで，暑熱環境下における高齢者の体内温度上昇および発汗量を解析した結果，平均皮膚温度および体内深部温度上昇において文献[10]での実測値とよい一致が得られている．
2.2.3 血流量のモデル化
成人の皮膚における血流定数 B は，

ベクトルスーパーコンピュータ SX-ACE による暑熱環境下体温上昇の高速解析

B(r )  [ B0 (r )  FHB (TH  TH 0 )  FSB  TS ]  2— 51 —(T ( r )T0 ( r ))6(10)で表される[11]．FHB , FSB はそれぞれ 17,500 W/(m3·oC2)， 1100 W/(m3·oC2)という値を有する．また，T (r ) ，T0 (r ) は，皮膚の微小組織における温度と基準温度である．さらに， TS は皮膚の平均温度上昇である．なお，本研究では FDTD 法の時間ステップ Δt および空間差分⊿を，それぞれ 2.0[s]，2.0[mm](8 ヶ月幼児では空間差分 Δ を 1.7[mm])とした．また，皮膚以外の血流定数 B は，

B((r ),T (r ))  B0 (r )T (r )  39o C(11)で表される[6]．式中の T は組織の温度であり，B0 は組織の基準となる血流量であるため，皮膚以外の血流定数 B は一定であると取り扱う．

2.3 解析のフローチャート
上記 2.2 の計算を含む体温上昇解析のフローチャートを図 2 に示す．図より，解析を実行すると，始めに 28 oC 環境下における熱定常状態の初期温度分布を読み込む．次に，人体モデルや発汗モデルの設定，また外気温などの環境設定を入力する．それから，2.2.1 で示した生体熱輸送方
程式による温度上昇を計算，またその計算を行う際，2.2.2，2.2.3 で示した発汗率の変化や血流量の変化など，熱調整機能を組み込んだ計算を行う．このフローチャートに従って，解析を行った．
Start the program

Initialize calculation to obtainthermoneutral state at 28C

Input parameters
Input model data

Loop through Time

Calculate temperature, Eqs. (1), (2)TB,m Eqs. (3), (4)Update thermoregulatory response• Sweating rate, Eqs. (5) – (11)• Blood perfusion, Eq. (12)Thermal conductivity, Specific Heat
Blood perfusion rate

Input sweating model
High, Normal, Low

Setup ambient condition
If Time= Simulation timeNo

Temperature, Humidity

Yes
End simulation (Store results)
temperature distribution, total sweating

図 2. 体温上昇解析のフローチャート
2.4 解析プログラムの最適化
2.4.1 ベクトル化
今回解析に用いたスーパーコンピュータ SX-ACE はベクトル処理を活用すると，高い性能を得ることができる．このことから，ベクトル計算に適したコードを作成する必要がある．ベクトル計算に適したループの例を図 3 に示す．言語は fortran90 にて記載する．図に示すように，最内側

— 52 —SENAC Vol. 48, No. 4（2015. 10）のループ長を最長にすることで，ノード当たりの計算要素数を増加させることが可能となる．
このように，さまざまなループレベル最適化[12]を施すことで，ベクトル化率とノードあたりの性能を向上させた．

MODELX=160
MODELY=320
MODELZ=866

MODELX=866
MODELY=320
MODELZ=160

DO K = 1,MODELZDO J = 1, MODELY
DO I = 1, MODELXTEMP(I,J,K)
END DOEND DOEND DODO I = 1,MODELZDO J = 1, MODELYDO K = 1, MODELXTEMP(K,J,I)
END DOEND DOEND DO

図 3. ベクトル化に適したループ
2.4.2 並列化
コードを並列化する際，MPI(Message Passing Interface)を用いて並列化を行った．MPI とは，分散メモリ間のデータ処理分割やデータ転送を担うメッセージパッシングの標準規格であり，コードから呼び出すサブプログラムのライブラリである．MPI を用いて並列化を行った理由として，並列化手法には，タスク並列用の自動並列化や OpenMP(Open Multi-Processing)，データ並列用のMPI などがある．今回，解析コードを並列化する上で，タスク並列と MPI を併用したハイブリッ
ド MPI，MPI のみを用いたフラット MPI をそれぞれ実装して，実行時間を比較した．その結果，フラット MPI の方が効率よく計算を高速化できたため，フラット MPI を用いた．
具体的に行った並列化処理としては，人体モデルを Z 軸方向各平面に分割し，それぞれの 1 平面を 1 プロセスとして割り当てた．

3. 解析結果異なる並列プロセス数における計算時間を表 1 に示す．表 1 より算出したプロセス数の増加に対する計算時間の加速率を図 4 に示す．図より，並列プロセス数の増加に伴って並列処理効率が顕著に低下することはなく，加速率は上昇している．しかしながら，理論値 64 倍に対して，並列プロセス数が 64 の場合に得られた加速率は，34.6 倍となった．これは主に演算量のインバランスによる並列化率の低下や，ノード間通信時のデータ転送量が増加する際に発生する遅延によるものと考えられる．

表 1. 異なる並列プロセス数における計算時間
並列プロセス数

148163264

計算時間(sec)

2329

721

39121312467

— 53 —ベクトルスーパーコンピュータ SX-ACE による暑熱環境下体温上昇の高速解析

70
Speed Up

60

Theoretical Value

Speed Up

50

Approxiamte Curve

4030201000816243240
Number of Process

485664図 4. 並列プロセス数の増加に対する計算時間の加速率

ベクトル化，並列化を行ったコードを SX-ACE に実装し，暑熱ばく露における体温上昇を解析した例を図 5 に示す．図は 3 時間の暑熱ばく露(37 oC)を行った際の体表面温度を表している．左から 3 歳児，22 歳相当の成人男性，65 歳相当の成人男性，75 歳相当の成人男性，成人女性となっている．
図より，他のモデルに比べて，75 歳相当の成人男性の体表面温度が最も大きく上昇していることが確認できる．これは，発汗機能の相違によるもので，高齢者は若年者よりも発汗量が少なく，発汗の開始が遅れるためである．

図 5. 暑熱ばく露における表面温度上昇の解析結果

— 54 —

SENAC Vol. 48, No. 4（2015. 10）

4. まとめ本研究では，SX-ACE による体温上昇解析の高速化を目標とし，コードのベクトル化，並列化を行い，異なる並列プロセス数における計算時間を算出した．算出した結果，プロセス数における加速率が理論値 64 倍に対して，計算時間の加速率は 34.6 倍となった．これは，逐次処理による並列化率の低下や，データ転送時の遅延によるものと考えられる．また，解析例として暑熱ばく露における年齢や性別など 5 体の人体モデルに対して体表面温度上昇解析を行った．他のモデルに比べて，75 歳相当の成人男性が最も大きく上昇するという結果が得られた．これは，高齢者の発汗反応の遅延および発汗量に起因するものである．
現在，図 5 に示した体表面温度上昇の解析例などの体温上昇解析に必要な計算時間は，ばく露
時間によって差はあるものの，目標である 20 秒に達していない．目標を達成するためには更なる並列化率の向上が今後の課題となる．

謝辞
本研究は，東北大学サイバーサイエンスセンターのスーパーコンピュータ SX-ACE を利用することで実現でき，計算コードの高速化および MPI 並列化にあたっては，同センターの佐々木技術職員ならびに山下技術専門職員に多大なるご協力をいただきました．この場を借りて深謝致します．
参考文献[1] A. Hirata, T. Asano, O. Fujiwara : “FDTD analysis of body-core temperature elevation in children
and adults for whole-body exposure”, Phys. Med. Biol., vol.53. Pp.5223-5238, (2008).
[2] A. Hirata, T. Nomura and I. Laakso : “Computational estimation of decline in sweating in the elderly
from measured body temperatures and sweating for passive heat exposure,” Physiol. Meas., vol. 33,
pp.N51-60, (2012).
[3] R. Egawa, S. Momose, K. Komatsu, Y. Isobe, H. Takizawa, A. Musa, H. Kobayashi : “Early
Evaluation of the SX-ACE Processor”, Proc. Int’l Conf. High Pef. Comput. Storage & Anal., no.14,(2014).
[4] T. Nagaoka, S. Watanabe, K. Sakurai, E. Kunieda, M. Taki, Y. Yamanaka : “Development of realistic
high-resolution whole-body voxel models of Japanese adults male and female of average high andweight, and application of models to radio-frequency electromagnetic-field dosimetry”, Phys. Med.
Biol., vol.49, pp.1-15, (2004).
[5] T. Nagaoka, E. Kunieda, and S. Watanabe : “Proportion-corrected scaled voxel models for Japanesechildren and their application to the numerical dosimetry of specific absorption rate for frequencies
from 30 MHz to 3 GHz”, Phys. Med. Biol., vol.53, no.23, pp.6695-6712, (2008).
[6] H. H. Pennes : “Analysis of tissue and arterial blood temperatures in resting forearm”, J.Appl.Physiol.,vol.1, pp.93-122, (1948).
[7] A. Hirata, T. Asano and O. Fujiwara : “FDTD analysis of human body-core temperature elevation due
to RF far-field energy prescribed in the ICNIRP guidelines,” Phys. Med. Biol., vol.52, pp.5013-23,(2007).
[8] D. Fiala, K. J. Lomas, and M. Stohrer : “Computer prediction of human thermoregulation andtemperature responses to a wide range of environmental conditions,” Int J Biometeorol, vol.45,
pp.143-159, (2001).
[9] A. Dufour and V. Candas : “Ageing and thermal responses during passive heat exposure: sweatingand sensory aspects,” Eur. J. Applied Physiol., vol. 100, pp.19-26, (2007).
[10] Y. Inoue and M. Shibasaki : “Regional differences in age-related decrements of the cutaneous
vascular and sweating responses to passive heating,” Eur. J. Appl. Physiol. & Occup. Physiol., vol.74,pp.78-84, (1996).
[11] P. Bernardi, M. Cavagnaro, S. Pisa and E. Piuzzi : “Specific absorption rate and temperature elevation
in a subject exposed in the far-field of radio-frequency sources operating in the 10-900-MHz range,”
IEEE Trans. Biomed. Eng., vol.50, pp.295-304, (2003).
[12] 王 建青, 岡部 公起, 曽根 秀昭, 藤原 修 : “ベクトル化・並列化に適した FDTD プログラミ
ングの一手法” ,信学技報.,Vol.103, No.223, pp31-35, (2003).