宇宙航空研究開発機構研究開発報告 JAXA Research and Development Report

宇宙航空研究開発機構研究開発報告

JAXA Research and Development Report

JAXA極超音速風洞・超音速風洞シュリーレン画像 へのCTの適用

藤井　啓介，平林　則明，小山　忠勇，津田　尚一，中川　宗敬，

板橋　幸広，中村　晃祥，岡林　希依，高間　良樹

2015年3月

宇宙航空研究開発機構

Japan Aerospace Exploration Agency

ISSN 1349-1113 JAXA-RR-14-011

本印刷物は、グリーン購入法に基づく基本方針の判断基準を満たす紙を使用しています。

宇宙航空研究開発機構研究開発報告

JAXA-RR-14-011

This document is provided by JAXA.

目　次

1

2

2.1 Computerized tomography ... 2

2.1.1 FBP

2.1.2 ART

2.2

3

3.1 0.5 m/1.27 m

3.2 1 m × 1 m

3.3

4

4.1 HB-2

4.2 FY11HOPE

4.3 AGARD-B

4.4 HST

4.5 HRV

5

/

This document is provided by JAXA.

JAXA 極超音速風洞・超音速風洞シュリーレン画像への CT の適用

藤井啓介、平林則明、小山忠勇、津田尚一、中川宗敬、板橋幸広、中村晃祥、

岡林希依、高間良樹

Application of Computerized Tomography onto set of regular schlieren images taken at JAXA Hypersonic / Supersonic Tunnels

Keisuke Fujii, Noriaki Hirabayashi, Tadao Koyoma, Shoichi Tsuda, Muneyoshi Nakagawa, Yukihiro Itabashi, Akiyoshi Nakamura, Kie Okabayashi and Yoshiki Takama

ABSTRACT

Although the schlieren visualization has been indeed one of so commonly utilized techniques to observe flow field at a high speed regime, it gives only integrated information related to the density gradient along the ray. In order to reconstruct three dimensional structure from a set of images of integrated information, the computerized tomographic (CT) algorithm is applied to them. Evaluation tests are conducted both at JAXA 1.27 m/0.5 m hypersonic wind tunnel and at JAXA 1 m × 1 m supersonic wind tunnel, which result in tomographic images where three dimensional structure of shock waves and expansion waves are easily recognized. It should be noted that three dimensional shock structures are reconstructed from ordinary equipments for industrial wind tunnels without any major modifications.

概要

シュリーレン可視化法は高速気流観測に非常に広く使われている密度勾配に関した計測手法の一つである が、密度勾配の観測線上にわたる積分値を与えるという性質上

3

computerized tomography(CT)

1.27m

1m × 1m

膨張波などの三次元構造の把握に有用な定性的再構成データの取得できることを確認した。

目 次

1

2

2

3

2.1 Computerized tomography . . . . 3

2.1.1 FBP

. . . . 3

2.1.2 ART

. . . . 5

2.2

CT

. . . . 7

JAXA 極超音速風洞・超音速風洞シュリーレン画像への CT の適用 藤井啓介^*1、平林則明^*1、小山忠勇^*1、津田尚一^*1、中川宗敬^*2、板橋幸広^*2、中村晃祥^*2、 岡林希依^*3、高間良樹^*3

Application of Computerized Tomography onto set of regular schlieren images taken at JAXA Hypersonic / Supersonic Tunnels Keisuke Fujii

, Noriaki Hirabayashi

, Tadao Koyoma

, Shoichi Tsuda

, Muneyoshi Nakagawa

, Yukihiro Itabashi

, Akiyoshi Nakamura

, Kie Okabayashi

and Yoshiki Takama

ABSTRACT Key Word:

technology 概要 JAXA 極超音速風洞・超音速風洞シュリーレン画像への CT の適用

Application of Computerized Tomography onto set of regular schlieren images taken at JAXA Hypersonic / Supersonic Tunnels Keisuke Fujii, Noriaki Hirabayashi, Tadao Koyoma, Shoichi Tsuda, Muneyoshi Nakagawa, Yukihiro Itabashi, Akiyoshi Nakamura, Kie Okabayashi and Yoshiki Takama ABSTRACT Although the schlieren visualization has been indeed one of so commonly utilized techniques to observe flow field at a high speed regime, it gives only integrated information related to the density gradient along the ray. In order to reconstruct three dimensional structure from a set of images of integrated information, the computerized tomographic (CT) algorithm is applied to them. Evaluation tests are conducted both at JAXA 1.27 m/0.5 m hypersonic wind tunnel and at JAXA 1 m × 1 m supersonic wind tunnel, which result in tomographic images where three dimensional structure of shock waves and expansion waves are easily recognized. It should be noted that three dimensional shock structures are reconstructed from ordinary equipments for industrial wind tunnels without any major modifications.

3

computerized tomography(CT)

1.27m

1m × 1m

1

2 2

3 2.1 Computerized tomography . . . . 3

2.1.1 FBP

. . . . 3

2.1.2 ART

. . . . 5

2.2

CT

. . . . 7

*

26

12

19

*1 航空本部　風洞技術開発センター

(Wind Tunnel Technology Center, Institute of Aeronautical Technology)

*2 一般財団法人　航空宇宙技術振興財団

(Japan AeroSpace Technology Foundation)

*3 航空本部　空力技術研グループ

(Aerodynamics Research Group, Institute of Aeronautical Technology)

This document is provided by JAXA.

3

9

3.1 0.5 m/1.27 m

. . . . 9

3.2 1 m × 1 m

. . . . 10

3.3

. . . . 11

4

14 4.1 HB-2

HWT2

. . . . 16

4.2 FY11HOPE

HWT2

. . . . 21

4.3 AGARD-B

SWT1

. . . . 25

4.4 HST

HWT1

. . . . 28

4.5 HRV

HWT2RCS

. . . . 29

5

31

1

遷音速以上の速度を持つ圧縮性流体における流れ場の把握は、学術的必要性のみでなく航空宇宙機の設 計情報取得といった工学的観点からも重要な要素である。その様な観点から流れ場の可視化技術としてシュ リーレン法やシャドーグラフ法などが古くからまた極めて幅広く活用されてきており、JAXA 0.5 m/1.27 m 極超音速風洞や

JAXA 1 m × 1 m

3

CFD

気流密度勾配の光路にわたった積分値に従った濃淡としての情報のみを与え、断面内の密度勾配分布を直接 与えないためである。そのような複雑形状における流れ場再構成のための手法の開発はこれまでに多くの 研究者らにより進められており、それらの例としてフォーカシングシュリーレンや、PIV

^[ 1 ]

^[ 2 ]

PLIF ^[ 3 ^]

Computerized Tomography (CT)

X

X

CT

MRI(Magnetic Resonance Imaging)

Positron Emission Tomography (PET)、超音波利用など様々な適

CT

CT

[4]、開

CT

X

JAXA 0.5 m/1.27 m

(HWT1/2)

JAXA 1 m × 1 m

(SWT1)

この報告では、それら開発試験を想定した模型周り可視化を行った際の、試験・計測手法としての成果を 取りまとめ、設備背景、計測手法、データ処理、結果及び課題をそれぞれ記述、整理することで今後の発 展、改良あるいは他風洞への適用に際しての情報提供も目的とする。そのため、本報告では今回使用した

CT

2

2.1 Computerized tomography

風洞シュリーレン画像へ

CT

Filtered BackProjection(FBP)

Algebraic Reconstruction Techniques(ART)

[5]

2.1.1 FBP

数値的ではない

‘traditional tomography’

1

C

f C

X

t

X t

A

A t

A

A

A

C

Backprojection

2

θ

K

l ^�

L

(r, φ)

L

(l ^� , θ)

l ^� = r cos (φ − θ)

であり、ここで空間分布関数

f (r, φ)、その投影データ（f

p(l, θ)

[ R f ] (l, θ) =

∞

−∞

f

l ² + z ² , θ + tan ⁻¹ z l

dz (l � = 0)

[ R f ] (0, θ) =

∞

−∞

f

z, θ + π 2

dz

(1)

図

1: Traditional Linear Tomography (Herman [5] p.29

2.1　Computerized tomography

　風洞シュリーレン画像へ CT 手法を適用するに当たり、当初は計算時間が比較的短くてすみ且つ解析的 バックグランドを持つ代表的な手法である Filtered BackProjection(FBP) 法を用いた。その後に、全周 囲計測をしない場合、模型背景を活用しようとする場合、屈折方向を積極的に活用する場合などを将来検 討するに当たり、より自由度の大きな

Algebraic Reconstruction Techniques(ART) 法を採用することと

1　序

2　背景

宇宙航空研究開発機構研究開発報告　JAXA-RR-14-011

2

This document is provided by JAXA.

展、改良あるいは他風洞への適用に際しての情報提供も目的とする。そのため、本報告では今回使用した

CT

2

2.1 Computerized tomography

風洞シュリーレン画像へ

CT

Filtered BackProjection(FBP)

Algebraic Reconstruction Techniques(ART)

[5]

2.1.1 FBP

数値的ではない

‘traditional tomography’

1

C

f C

X

t

X t

A

A t

A

A

A

C

Backprojection

2

θ

K

l ^�

L

(r, φ)

L

(l ^� , θ)

l ^� = r cos (φ − θ)

であり、ここで空間分布関数

f (r, φ)、その投影データ（f

p(l, θ)

[ R f ] (l, θ) =

∞

−∞

f

l ² + z ² , θ + tan ^{− 1} z l

dz (l � = 0)

[ R f ] (0, θ) =

∞

−∞

f z, θ + π

2

dz

(1)

図

1: Traditional Linear Tomography (Herman [5] p.29

2.1.1　FBP 法

図

2: The relationship between the (r, φ) space and the(l, θ) space (Herman [5] p.104

(Radon

Backprojection

(r, φ)

(p)

π

0

p(l, θ)dθ (2)

という操作に相当する。この

Backprojection

Radon

(式 (1))

Backprojection

Filtered Back Projection (FBP)

p(l, θ)

f (r, φ)

Radon

R ⁻¹

f (r, φ) = [ R ⁻¹ p] (r, φ)

= 1 2π ²

π

0

E

− E

1 r cos (θ − φ) − l

∂p (l, θ)

∂l dldθ (3)

となることが知られている。この

Radon

R ^{− 1}

q(l, θ) = [ D ^Y p] (l, θ) = ∂

∂l p (l, θ) (4)

t(l ^� , θ) = [ H ^Y q] (l ^� , θ) = − 1 π

∞

−∞

q (l, θ)

l ^� − l dl (5)

[ B t] (r, φ) =

π

0

t (r cos (θ − φ) , θ) dθ (6)

R ⁻¹ = − 1

2π BH ^Y D ^Y (7)

ここで

(4)

Hilbert

(6)

Backprojection

(式 (2))

バックプロジェクション

上述の

Backprojection

θ = mΔ, l = nd（但し m, n

p

[ B t] (r, φ) =

π

0

t (r cos (θ − φ) , θ) dθ

≈ Δ

M − 1

m=0

t (r cos (mΔ − φ) , mΔ)

≈ Δ

M−1

m=0

(n + 1) d − r cos (mΔ − φ)

d t (nd, mΔ) + r cos (mΔ − φ) − nd

d t ((n + 1) d, mΔ)

(8)

3

This document is provided by JAXA.

図

2: The relationship between the (r, φ) space and the(l, θ) space (Herman [5] p.104

(Radon

Backprojection

(r, φ)

(p)

π

0

p(l, θ)dθ (2)

という操作に相当する。この

Backprojection

Radon

(式 (1))

Backprojection

Filtered Back Projection (FBP)

p(l, θ)

f (r, φ)

Radon

R ^{− 1}

f (r, φ) = [ R ^{− 1} p] (r, φ)

= 1 2π ²

π

0

E

− E

1 r cos (θ − φ) − l

∂p (l, θ)

∂l dldθ (3)

となることが知られている。この

Radon

R ^{− 1}

q(l, θ) = [ D ^Y p] (l, θ) = ∂

∂l p (l, θ) (4)

t(l ^� , θ) = [ H ^Y q] (l ^� , θ) = − 1 π

∞

−∞

q (l, θ)

l ^� − l dl (5)

[ B t] (r, φ) =

π

0

t (r cos (θ − φ) , θ) dθ (6)

R ^{− 1} = − 1

2π BH ^Y D ^Y (7)

ここで

(4)

Hilbert

(6)

Backprojection

(式 (2))

バックプロジェクション

上述の

Backprojection

θ = mΔ, l = nd（但し m, n

p

[ B t] (r, φ) =

π

0

t (r cos (θ − φ) , θ) dθ

≈ Δ

M −1

m=0

t (r cos (mΔ − φ) , mΔ)

≈ Δ

M − 1

m=0

(n + 1) d − r cos (mΔ − φ)

d t (nd, mΔ) + r cos (mΔ − φ) − nd

d t ((n + 1) d, mΔ)

(8)

とかけ、分布関数を得る。

コンボリューション

上記

BackProjection

(7)

Hilbert

Hilbert

H Φ](l ^� ) = − π ¹

Φ(l)

l

− l dl）は畳み込み積分を使って以下のようにもかける:

[ H Φ] = Φ ∗ ρ ρ(l) = − 1

πl

が、ここで上記の

Hilbert

l ^� = l

‘improper integral’

‘regularization’

ρ A

A→∞ lim Φ ^∗ ρ A = H Φ

これを満たす

regularization family

ρ A (u) = − 2

A/2

0

F A (U ) sin (2πU u) dU (9)

ここで

F A (U )

[5]Table 8.1

A (U ) ≡ α + (1 − α) cos ^2πU _A

∂

∂l ^� p ∗ ρ A

(l ^� ) =

∞

−∞

p ^� _θ (l)ρ A (l ^� − l)dl

= −

∞

−∞

p θ (l)ρ ^� _A (l ^� − l)dl

ρ A

(9)

ρ ^� _A (u) = − 4π

∞

−∞

U F A (U ) cos(2πU u)dU

q(u)

q(u) ≡ − 1

2π ρ ^� _A (u) = 2

∞

−∞

U F A (U ) cos(2πU u)dU (10)

と定義することで、Radon逆変換は、

R ⁻¹ p

(r, φ) ≈ f ^∗ (r, φ) = B ([p θ ∗ q] (l))

= B

_∞

−∞

p(l, θ)q(l ^� − l)dl

(11)

と表されることが分かる。そのため、式

(8)

p(l, θ)

f (r, φ)

Radon

(10)

q

2.1.2 ART

上述の

FBP

180

Algebraic Reconstruction Techniques (ART)

宇宙航空研究開発機構研究開発報告　JAXA-RR-14-011

4

This document is provided by JAXA.

とかけ、分布関数を得る。

コンボリューション

上記

BackProjection

(7)

Hilbert

Hilbert

H Φ](l ^� ) = − π ¹

Φ(l)

l

− l dl）は畳み込み積分を使って以下のようにもかける:

[ H Φ] = Φ ∗ ρ ρ(l) = − 1

πl

が、ここで上記の

Hilbert

l ^� = l

‘improper integral’

‘regularization’

ρ A

A→∞ lim Φ ^∗ ρ A = H Φ

これを満たす

regularization family

ρ A (u) = − 2

A/2

0

F A (U ) sin (2πU u) dU (9)

ここで

F A (U )

[5]Table 8.1

A (U ) ≡ α + (1 − α) cos ^2πU _A

∂

∂l ^� p ∗ ρ A

(l ^� ) =

∞

−∞

p ^� _θ (l)ρ A (l ^� − l)dl

= −

∞

−∞

p θ (l)ρ ^� _A (l ^� − l)dl

ρ A

(9)

ρ ^� _A (u) = − 4π

∞

−∞

U F A (U ) cos(2πU u)dU

q(u)

q(u) ≡ − 1

2π ρ ^� _A (u) = 2

∞

−∞

U F A (U ) cos(2πU u)dU (10)

と定義することで、Radon逆変換は、

R ^{− 1} p

(r, φ) ≈ f ^∗ (r, φ) = B ([p θ ∗ q] (l))

= B

_∞

−∞

p(l, θ)q(l ^� − l)dl

(11)

と表されることが分かる。そのため、式

(8)

p(l, θ)

f (r, φ)

Radon

(10)

q

2.1.2 ART

上述の

FBP

180

Algebraic Reconstruction Techniques (ART)

例えば、空間を

J

J-dim

x x x

I

I-dim

yyy

× J

R (projection matrix)

yyy = R x x x (12)

で表される。計測結果である

yyy

x x x

x x x ^(o)

x x x ^(k)

k + 1

x x x ^(k+1)

yyy

i

y i

R x x x ^(k+1)

i

(rrr i , x x x ^(k+1) )

x x x ^(k)

rrr i

R ^T

i

x x x ^(k)

x x x ^(k+1)

i

x x x ^(k+1)

x x

x ^(k+1) = x x x ^(k) + y i − (rrr i , x x x ^(k) ) (rrr i , rrr i ) rrr i

で得られる。このステップでは

x x x

yyy

x x x

しかし

R

yyy

x x x

そこで基礎式

(12)

eee ∈ Range( R ) ^⊥

yyy = eee + R x x x (13)

ここで、u

u u = reee(r

≡ x x x − μ x x x

x x x

yyy = 1

r u u u + R (zzz + μ x x x ) r (yyy − R μ x x x ) = u u u + r R zzz ( I , r R )

u u u zzz

= r (yyy − R μ x x x ) (14)

となるので、(

I , r R )

u u u zzz

を一つの列ベクトルと考えると、上記手法をそのまま適用する ことで、

u u u ^(k+1) zzz ^(k+1)

− u u u ^(k)

zzz ^(k)

=

r { y i − (rrr i , μ x x x ) } −

ttt, u u u ^(k)

zzz ^(k)

(ttt, ttt) ttt

= r { y i − (rrr i , x x x) } − u ^(k) _i 1 + r ² | rrr i | ²

iii i

rrrr i

(15)

となる。ここで、

I

i

I

i

≡ iii i

rrrr r

とする。よって、u

u u

x x

x = zzz + μ x x x

u u u ^(k+1) − u u u ^(k) = c ^(k) iii i (16)

x x

x ^(k+1) − x x x ^(k) = rc ^(k) rrr i (17)

(18) 2.1.2　ART 法

JAXA 極超音速風洞 ・ 超音速風洞シュリーレン画像への CT の適用

5

This document is provided by JAXA.

であったり、投影画像が単純な分布関数の線積分で表すことができない場合にも適用可能であることなどか ら、その有効性が認められている手法である。

例えば、空間を

J

J-dim

x x x

I

I-dim

yyy

× J

R (projection matrix)

yyy = R x x x (12)

で表される。計測結果である

yyy

x x x

x x x ^(o)

x x x ^(k)

k + 1

x x x ^(k+1)

yyy

i

y i

R x x x ^(k+1)

i

(rrr i , x x x ^(k+1) )

x x x ^(k)

rrr i

R ^T

i

x x x ^(k)

x x x ^(k+1)

rrr i

x x x ^(k+1)

x x

x ^(k+1) = x x x ^(k) + y i − (rrr i , x x x ^(k) ) (rrr i , rrr i ) rrr i

で得られる。このステップでは

x x x

yyy

x x x

しかし

R

yyy

x x x

そこで基礎式

(12)

eee ∈ Range( R ) ^⊥

yyy = eee + R x x x (13)

ここで、u

u u = reee(r

≡ x x x − μ x x x

x x x

yyy = 1

r u u u + R (zzz + μ x x x ) r (yyy − R μ x x x ) = u u u + r R zzz ( I , r R )

u u u zzz

= r (yyy − R μ x x x ) (14)

となるので、(

I , r R )

u u u zzz

を一つの列ベクトルと考えると、上記手法をそのまま適用する ことで、

u u u ^(k+1) zzz ^(k+1)

− u u u ^(k)

zzz ^(k)

=

r { y i − (rrr i , μ x x x ) } −

ttt, u u u ^(k)

zzz ^(k)

(ttt, ttt) ttt

= r { y i − (rrr i , x x x) } − u ^(k) _i 1 + r ² | rrr i | ²

iii i

rrrr i

(15)

となる。ここで、

I

i

I

i

≡ iii i

rrrr r

とする。よって、u

u u

x x x = zzz + μ x x x

u

u u ^(k+1) − u u u ^(k) = c ^(k) iii i (16)

x x x ^(k+1) − x x x ^(k) = rc ^(k) rrr i (17)

(18)

c ^(k) ≡ λ ^(k) r y i −

rrr i , x x x ^(k)

− u ^(k) _i 1 + r ² | rrr i | ²

である。λ

^(k)

< λ ^(k) ≤ 1

λ ^(k) = 0.05

重み係数

R

ART

(12)

projection matrix R

f

y =

_∞

−∞

f (x x x)dl

y i =

j

x j dl ij (19)

であるため、この重み係数はそれぞれの点（領域）によって区切られる積分経路の長さとすればよい

R ij ≡ dl ij

[5]

blob

3

2N × 2N

1 ∼ m 0

Line-i

(j = 2N n + m)

R ^ij

Line-i

(m 0 , n)

R ^{i,2N n+m}

= (m(n) − m 0 )/ sin θ、

R ^{i,2N n+m} = 1/ sin θ, m 1 < m < m 0

R ^{i,2N n+m}

= (m(n + 1) − m 1 )/ sin θ

n

0

Line-i

0

R

0

x x

(17)

0

rrr i

i

3

4N ²

× 2N

2.2

CT

シュリーレン可視化画像に

CT

CT

X

また、シュリーレン可視化により得られる画像は、X線

CT

y =

dε =

∞

−∞

dε dl dl

u

ε u = l

−∞

1 n

∂n

∂u ds ≈ 1 n 0

l

−∞

∂n

∂u ds n = Gρ + 1

G

Gradstone-Dale

(s

(x

(u

_∂u ^∂ρ = ^∂ρ _∂x

宇宙航空研究開発機構研究開発報告　JAXA-RR-14-011

6

This document is provided by JAXA.

ここで、

c ^(k) ≡ λ ^(k) r y i −

rrr i , x x x ^(k)

− u ^(k) _i 1 + r ² | rrr i | ²

である。λ

^(k)

< λ ^(k) ≤ 1

λ ^(k) = 0.05

重み係数

R

ART

(12)

projection matrix R

f

y =

∞

−∞

f (x x x)dl

y i =

j

x j dl ij (19)

であるため、この重み係数はそれぞれの点（領域）によって区切られる積分経路の長さとすればよい

R ij ≡ dl ij

[5]

blob

3

2N × 2N

1 ∼ m 0

Line-i

(j = 2N n + m)

R ^ij

Line-i

(m 0 , n)

R ^{i,2N n+m}

= (m(n) − m 0 )/ sin θ、

R ^{i,2N n+m} = 1/ sin θ, m 1 < m < m 0

R ^{i,2N n+m}

= (m(n + 1) − m 1 )/ sin θ

n

0

Line-i

0

R

0

x x

(17)

0

rrr i

i

3

4N ²

× 2N

2.2

CT

シュリーレン可視化画像に

CT

CT

X

また、シュリーレン可視化により得られる画像は、X線

CT

y =

dε =

∞

−∞

dε dl dl

u

ε u = l

−∞

1 n

∂n

∂u ds ≈ 1 n 0

l

−∞

∂n

∂u ds n = Gρ + 1

G

Gradstone-Dale

(s

(x

(u

^∂ρ _∂u = ^∂ρ _∂x

0

2N

0 n n + 1 2N

m

m

m(n)

m(n + 1) Line-i

θ

図

3:

の関数となるが、投影面内の

t

^∂ρ _∂t = cos θ ^∂ρ _∂y + sin θ ^∂ρ _∂z

θ

_∂x ^∂ρ

_∂x ^∂ρ

^∂ρ _∂t

s − t

y − z

θ

I t

の

t

dI t

dt ∝

_∞

−∞

∂ ² ρ

∂t ² ds =

_∞

−∞

∂ ² ρ

∂t ² + ∂ ² ρ

∂s ²

ds

=

∞

−∞ ∇ ² ρds

という関係から、

∇ ² ρ = ^∂ _∂y

^ρ

+ ^∂ _∂z

^ρ

通常のシュリーレン光学系では

^∂ρ _∂x

CT

FBP

シュリーレン画像に

CT

CT

‘metal artifact’

[6]。シュリーレン画像を基にし

2.2　シュリーレン可視化画像への CT 適用

JAXA 極超音速風洞 ・ 超音速風洞シュリーレン画像への CT の適用

7

This document is provided by JAXA.