第 27 回数値流体力学シンポジウム
A10-2

回転するタイヤ周りに発生する空力音の数値解析（Re=100,000）
Computational analysis of aeroacoustic noise generated by rotating tire (Re=100,000)
阿部圭晃，東大院, 相模原市中央区由野台 3-1-1，E-mail：abe@ﬂab.isas.jaxa.jp

野々村 拓，JAXA/ISAS, 相模原市中央区由野台 3-1-1，E-mail：nonomura@ﬂab.isas.jaxa.jp
近藤 勝俊，東理大, 相模原市中央区由野台 3-1-1，E-mail：kondo@ﬂab.isas.jaxa.jp
飯田 大貴，東大院, 相模原市中央区由野台 3-1-1，E-mail：iida@ﬂab.isas.jaxa.jp
渡辺 毅，JAXA/ISAS, 相模原市中央区由野台 3-1-1，E-mail：watanabe@ﬂab.isas.jaxa.jp
池田 俊之，横浜ゴム, 神奈川県平塚市追分 2-1,
小石 正隆，横浜ゴム, 神奈川県平塚市追分 2-1,
山本 誠，東理大, 東京都千代田区九段北 1-14-6, E-mail：yamamoto@rs.kagu.tus.ac.jp
藤井孝藏，JAXA/ISAS, 相模原市中央区由野台 3-1-1，E-mail：fujii@ﬂab.isas.jaxa.jp

Yoshiaki Abe, University of Tokyo, 3-1-1, Yoshinodai, Sagamihara, Japan E-mail：abe@ﬂab.isas.jaxa.jp
Taku Nonomura, JAXA/ISAS, 3-1-1, Yoshinodai, Sagamihara, Japan E-mail：nonomura@ﬂab.isas.jaxa.jp
Katsutoshi Kondo, 3-1-1, Yoshinodai, Sagamihara, Japan E-mail：kondo@ﬂab.isas.jaxa.jp
Daiki Iida, University of Tokyo, 3-1-1, Yoshinodai, Sagamihara, Japan E-mail：iida@ﬂab.isas.jaxa.jp
Takeshi Watanabe, JAXA/ISAS, 3-1-1, Yoshinodai, Sagamihara, Japan E-mail：watanabe@ﬂab.isas.jaxa.jp
Toshiyuki Ikeda, The Yokohama Rubber Corporation, 2-1 Oiwake Hiratsuka Kanagawa, Japan
Masataka Koishi, The Yokohama Rubber Corporation, 2-1 Oiwake Hiratsuka Kanagawa, Japan
Makoto Yamamoto, 3-1-1, Yoshinodai, Sagamihara, Japan E-mail：Kozo Fujii, JAXA/ISAS, 3-1-1, Yoshinodai, Sagamihara, Japan E-mail：fujii@ﬂab.isas.jaxa.jp
The acoustic wave generated from a rotating tire is analyzed by a large-eddy simulation (LES) of a threedimensional compressible Navier-Stokes equations. The Reynolds number is based on the diameter of the tire
and the rotating speed, and set to 1.0 × 105 . The rotating speed is set to Ma= 0.1, and the ground moves withthe corresponding speed. In this paper, the following topics are focused and discussed: 1) Three-dimensional
fine vortex structures appear near the surface of the rotating tire and its contact area on the ground. 2) The
characteristic acoustic wave is extracted in a high-frequency domain (around 1000Hz). 3) POD is conducted todetect the source of the characteristic acoustic wave.

1.

諸言近年，自動車のエンジン等から発生する騒音レベルが減少し，それに伴いタイヤから発生する騒音レベルが相対的に増加した．タイヤの騒音に関しては，1000Hz 付近の音波が支配的であり，その音波はタイヤ表面のパターンから発生する事が知られている．そこで，騒音の少ない新しいタイヤモデルを考える上で，タイヤ表面のパターンから発生する騒音に強い関心が持たれており，騒音の発生メカニズムの理解が望まれている．タイヤのパターンから発生する騒音としては，空力音と振動音が考えられる．実験では，音波の発生要因が流体変動である空力音か，構造振動によるものかを見分ける事は難しく，これらの寄与分は十分に明らかになっていない．本研究ではこのうち流体変動を発生要因とする空力音に着目する．
我々の研究グループでは過去に，レイノルズ数 O(104 )の場合にタイヤパターンから発生する空力音の解析を行った(1) ．解析には高次精度スキーム（6 次精度コンパクト差分）を用い，LES による直接解法を適用した．過去の報
告は，計算コストの兼ね合いから意図的にレイノルズ数を下げた（実際の解析対象とするタイヤ周りの流れはレイノルズ数 O(106 )）解析であった．今回は，レイノルズ数をより実際に近づけた解析（O(105 )）においてパターンの無いタイヤ周りの流体計算を行い，空力音の発生メカニズムを議論する．

図 1: Acoustic waves and ﬁne vortex structures arounda rotating tire (Re=105 )解析対象は無響室での実験を想定し，回転するドラム

1

Copyright c 2013 by JSFM

第 27 回数値流体力学シンポジウム
A10-2

上にタイヤが接し回転している状態とした．その為，実際の走行時にタイヤ周りに発生する，進行方向を風上とする一様流は存在しない．まず，Sec. 2. で本解析の計算条件（計算モデル，手法，格子）を述べる．次に Sec. 3.
で，得られた圧力場と渦構造を議論する．また各周波数帯の音圧レベル（SPL）分布を示し，固有直交分解（POD）による音源推定を行う．最後に，Sec. 4. で結論を述べる．

は 10 点重ねてある．また背景格子は直交の矩形格子であり，進行方向にはタイヤ径の 200 倍，鉛直・タイヤ軸方向には 100 倍の領域を採っている．なお，タイヤ周辺格子では地面以下の領域，背景格子では地面以下及びタイヤ内部の領域が実際には流体計算を行わない内部点として定義される．これらの格子の間で適切に物理量を内挿し合う事により，計算を進める（詳細は Sec.2.4 参照）．

解析モデルと手法

2.

2.1 解析モデル
最終的な解析対象は 122km/h 走行中の自動車のタイヤ（タイヤ径は 60cm）であり，走行速度とタイヤ径で定義されるレイノルズ数は 106 ，マッハ数は 0.1 である．しかし本研究では LES を用いた直接解法によりタイヤ周りの音響場を解析する事を目的としており，上記のレイノルズ数（106 ）で必要とされる計算格子は非常に大きくなる．その為，本論文では計算コストの観点から意図的にレイノルズ数を 105 に下げた状態を解析対象とする．ゆえに実際の走行時の流れ場を厳密には再現しないが，音響場と音源の推定を含めた物理現象の理解という観点からはこのレイノルズ数の流れ場の解析結果も有用な知見となる（過去にレイノルズ数 104 の解析も行っている(1) ）．
図 2 に実際の解析対象を示す．図中に灰色で塗りつぶされた領域（タイヤ表面及び地面上の矩形領域）にマッハ数 0.1 に相当する表面流速を与え，地面移動とタイヤ回転を模擬している．繰り返しになるが，本解析では実際の走行時の状態（図 2 の表面流速に加えてタイヤ進行方向を風上とする一様流が存在）ではなく，回転するドラム上にタイヤが接し回転している状態をモデル化する．

(a) Perspective view(b) Front view(c) Side view(d) Zoom view

図 3: Computational grids are shown. The tire grids
are colored black, and the background grids are colored
red. (a), (b) and (c) show every 10 to 20 grid lines. (d)shows actual grid lines.

zx: Solved on tire grid
: Solved on background grid(a) Side view(b) Perspective view

図 2: Computational models of a rotating tire and moving ground are shown. The grey colored regions indicate
moving surfaces, and the red arrows show the directionin which the tire and ground move.

図 4: The interpolated regions are shown. Red region
is solved on the tire grid and interpolated to the background one. Blue region is solved on the backgroundgrid and interpolated to the tire one.

2.2 計算格子
図 3 に計算格子を示す．本解析では音響場を直接解法により捉える為，コンパクト差分による高解像度スキームを用いた LES を行う．計算には構造格子を採用するが，タイヤと地面を含めた複雑形状に対応する為に重合格子法(2) を用いる．図 3 に示すようにタイヤ周辺格子（クロロ：1500 万点）と背景格子（赤色：5000 万点）に分ける．
タイヤ周辺格子は図 3(b) に示す L-O 型とし，周方向に2.3 解析手法
支配方程式は圧縮性 Navier-Stokes 方程式とし，タイヤ径 D，一様流密度，音速で無次元化を行う．移流項，粘性項，座標変換のメトリックおよびヤコビアンの計算には，空間 6 次精度のコンパクト差分法(3) を用いる．なおタイヤ周辺格子での地面以下の領域，背景格子でのタイヤ内部及び地面以下等の物体内部点を差分時のステン

2Copyright c 2013 by JSFM

第 27 回数値流体力学シンポジウム
A10-2

シルに含めない工夫を施している（境界付近で 2 次精度中心差分に切り替え）．さらに計算の安定化の為，壁面
境界 1 点目と 2 点目には 2 次精度陽的差分を適用する．
また，高周波の数値振動を抑える為，6 次精度 3 重対角フィルター（フィルター係数は 0.40）(3)(4) を用いる．時間方向の離散化は 2 次精度 3 点後退差分で行い，時間積分には ADI-SGS 陰解法(5)（内部反復は 5 回）を用いる．
無次元時間ステップは 5.0 × 10−3 とし，10 万ステップ（0 ≤ tu∞ /D ≤ 500）の計算を行う．またタイヤ周辺格子のタイヤ表面と，背景格子の地面移動領域には表面流速（マッハ数 0.1 に相当）をすべり境界条件として与える．
重合格子法において接地面付近で適切に内挿出来ない場所には，埋め込み境界法（Immersed boundary method:IBM）(6) を用いて物体内部の数点に鏡面的に流速を外挿する．

面近傍を，背景格子では地面近傍の流れ場を高い解像度で捉えられる．しかしタイヤと地面の接触する面（タイヤ接地面と呼ぶ）近傍の扱いは自明では無いので，図 4に示す．図 4 では，灰色はタイヤと地面以下の領域を表している．青色の領域は背景格子，赤色の領域はタイヤ周辺格子での解が高い解像度を持つと考え，各格子で内挿し合う．すなわち，背景格子の赤色の領域にはタイヤ周辺格子の解を内挿し，タイヤ周辺格子の青色の領域には背景格子の解を内挿する．特に接地面付近の内挿領域
ではタイヤ表面と地面から等距離の所に内挿の切り替え面を置き，接地面極近傍まで高い解像度を保った計算を行う．なお先にも述べたように物体境界に近い内部点には IBM に基づいた解が外挿されており，接地面付近でも安定して計算を進める事が出来る．

(a) Perspective view(a) Perspective view(b) Top view(b) Top view

図 6: Instantaneous ﬂow ﬁeld: isosurfaces of second invariant of a velocity gradient tensor are colored grey.

図 5: Instantaneous ﬂow ﬁeld: pressure distribution
around the rotating tire is shown.

流れ場と考察3.

2.4 タイヤ接地面付近の計算手法
構造格子を用いるメリットの１つとして，格子を物体形状に沿わせ物体境界近傍の解像度を高める事があげられる．重合格子法を用いると，複雑形状でも物体境界近傍は質の良い（境界に沿った）格子を用いた計算が可能である．実際，本解析ではタイヤ周辺格子ではタイヤ表以下では得られた計算結果を示す．

3.1

圧力場と渦構造

まず，図 5 に瞬時圧力場（Cp ）を示す．図 5(a) は向かって左下，図 5(b) は左が進行方向である（従って地面移動方向は右上及び右方向）．タイヤ表面，接地面近傍の

3

Copyright c 2013 by JSFM

第 27 回数値流体力学シンポジウム
A10-2

地面上，及びタイヤ後方の地面上の圧力場は複雑に振動している．また，図 6 には速度勾配テンソルの第二不変量の等値面を灰色で示している．接地面近傍には細かい3 次元渦構造が形成され，それらがタイヤ表面やタイヤ後方に流れていく．図 5 に示す複雑な圧力変動は，このような細かい渦構造による流体変動によって生じるものである．接地面前方では，地面の移動によって誘起される地面近傍の流れと，タイヤの回転によって生じたタイヤ表面近傍の流れがぶつかり合う．行き場を失った流れはタイヤの左右側方に回りこみ，細かい 3 次元渦構造を形成しながら後方へ流れていく．タイヤ後方の強い圧力振動が接地面を起点として左右に別れて拡がるのは，このような細かい渦構造の流れに起因する．なお，タイヤ後方で軸方向に長く伸びる圧力振動と渦構造は地面移動領域の境界で生じている．

(a) OASPL

3.2 音響場
次に，タイヤ表面，タイヤ側面から 10cm の鉛直面，地面上の音圧レベルにオクターブバンドフィルタをかけた結果を図 7-9 に示す．解析の対象は 1000Hz 以上の高周波帯である．まず図 7(a) の OASPL（172-1380Hz）に示
すように，タイヤ表面，接地面，タイヤ後方の地面上で特に音圧レベルが強い．これらはいずれも図 6 で指摘した細かい渦構造が生じている領域であり，非定常性の強い流体変動に伴う圧力振動を拾ったと考えられる．実際，これらの領域で音圧レベルが特に強く生じているのは低周波（172-345[Hz]）であり，解析対象とする高周波数帯に近づくほど弱くなる（図 7(b)-(e)）．このような傾向は鉛直面，地面上の音圧レベル分布を可視化した図 8,9 にも現れている．また，上で指摘した音圧レベルの高い領域は細かい渦構造の分布する領域に局在しており，渦構造の存在しない領域までは空間的に拡がっていない．
一方，図 7(d),(e) では，接地面付近を中心として同心円状に拡がる音響波を捉えている．これらの分布は前述の細かい渦構造の空間分布とは一致しておらず，さらに高周波帯（690-1380Hz）でも明確に現れる．また，空間的に広い範囲に拡がるという点で非常に特徴的な音響波である．特に，この同心円の中心は接地面前方にある事から，タイヤ後方よりも前方へ大きく拡がっていると考えられる（図 8(d),(e)，及び図 9(d),(e)）．なお図 9(a),(d)のタイヤ後方の音圧レベル分布で，タイヤ軸方向に軸を持つ振動が見られる．この領域では，タイヤ側方から離れる方向（図上下方向）の格子幅が大きく変化している為，細長い渦構造が定在してしまう事で数値振動が生じたと考えられる．しかし，次節で議論するように今回着目するタイヤ接地面前方付近を中心とする同心円状の音響波には影響しないと考えられる．
このような 1000Hz 以上で強く残る音響波は，接地面付近の格子の歪み等に由来する数値誤差を極力減らした
今回の計算で初めて捉えられたものである．次節では，この同心円状に拡がる特徴的な音響波の音源推定を目的として，周波数領域での固有直交分解（Proper orthogonaldecomposition: POD）を行う．
4(b) St = 3.2 (172 [Hz])(c) St = 6.4 (345 [Hz])(d) St = 12.8 (690 [Hz])(e) St = 25.6 (1380 [Hz])図 7: Sound pressure level in the perspective view (octave band ﬁltered one).
Copyright c 2013 by JSFM

第 27 回数値流体力学シンポジウム
A10-2(a) OASPL(a) OASPL(b) St = 3.2 (172 [Hz])(b) St = 3.2 (172 [Hz])(c) St = 6.4 (345 [Hz])(c) St = 6.4 (345 [Hz])(d) St = 12.8 (690 [Hz])(d) St = 12.8 (690 [Hz])(e) St = 25.6 (1380 [Hz])(e) St = 25.6 (1380 [Hz])図 8: Sound pressure level on the vertical plane at a
distance of 10 cm from the wheel (octave band ﬁlteredone).

5

図 9: Sound pressure level on the ground (octave bandﬁltered one).
Copyright c 2013 by JSFM

第 27 回数値流体力学シンポジウム
A10-2(a) Power ratio for POD modes at St = 2.56(a) Power ratio for POD modes at St = 12.8(b) 1st mode(b) 1st mode(c) 2nd mode

図 10: POD modes at St = 2.56 (138[Hz]) are shown,
where POD is conducted on the ground. Top, middle
and bottom ﬁgures show the real part, amplitude andphase of each mode, respectively. (a) shows the power
ratio for each mode. (b) shows the ﬁrst mode, and (c)shows the second one.

(c) 2nd mode

図 11: POD modes at St = 12.8 (694[Hz]) are shown,
where POD is conducted on the ground. Top, middle
and bottom ﬁgures show the real part, amplitude andphase of each mode, respectively. (a) shows the power
ratio for each mode. (b) shows the ﬁrst mode, and (c)shows the second one.

6Copyright c 2013 by JSFM

第 27 回数値流体力学シンポジウム
A10-2

3.3 周波数領域の POD
ここでは POD を用いて主要な音響波の抽出を試みる．
POD は乱流場から主要なエネルギーモードを抽出する際に使用される手法であり，今回は音響波の抽出に適用する(7)(8)(9) ．POD は一般に，時系列の流れ場データに対して行う事が多い．しかし，乱流を含む流れ場は非常に複雑なパターンを形成する為，元の流れ場を構築する為に非常に多くのモードを必要とし，主要なモードを特定する事が困難となる．これらの問題に対処する為，今回は周波数領域での POD を行う．まず，時系列流れ場（圧力変動量）のデータに短時間フーリエ変換を繰り返し，得られたフーリエ係数を snap shot として複素空間での POD を行う．この時，短時間フーリエ変換をかける範囲を指定する事によって，特定の周波数領域のエネルギーモードのみを抽出出来る（同時に位相の情報も得られる）．今回は，Hanning 窓を用いて着目する周波数（St 数）の周期 4 つ分を切り出し，短時間フーリエ変換を行った．使用する時系列データは同じなので，着目する周波数が小さい程多くの snap shot を得られ，信頼性の高い POD を行える．
以下に St = 2.56, 12.8, 25.6 の 3 つの周波数領域における結果を示す．いずれもタイヤ前方の一辺 O(D) 程度の長方形領域（図 10,11,12 の黒色の線で示す）内を積分領域とし，地面上の圧力変動分布から主要なモードを抽出している．各図の (a) には各モードの比率を示し，(b),(c) には1 番目と 2 番目に主要なモードを表している．図 (b),(c)は積分領域内の主要なモードに追従するモードを全領域で示しており，着目する周波数領域において，積分領域内の主要な圧力変動とその他の領域の圧力変動がどのように関係しているかを表す．
まず，St = 2.56 ではタイヤ接地面前方を中心とする同心円状の第 1 モード（図 10(b)）が卓越しており，位相図で分かるように接地面前方から遠方へ音響波が拡がっていく．第 2 モード（図 10(c)）ではタイヤを中心として左右に双極となる分布が現れており，位相図からタイヤの左右で交互に音圧レベルの高低領域が入れ替わる事が分かる．これらの事から，第 1 モードでは単極子音源，第 2モードでは二重極音源が同定される．単極子音源は流れ場の吸い込み・湧き出しに相当するものであり，接地面前方付近で流れが圧縮・膨張を繰り返している為に音が発生している可能性がある．また，二重極音源は円柱周りから発生するエオルス音等に代表され，物体にかかる流体力から音響波が発生するものである．ここでは着目している周波数が低い事からも乱流による発生音では無い事，大きな流れの構造が発生源となっている事が推測される．特に，第 1,2 モードで大きく特徴の異なるモードが抽出される事は音圧レベルの空間分布からは予想されなかった事である．
次に St = 12.8 でも，第 1 モード（図 11(b)）が接地面前方を中心とする同心円状の分布を持つ一方，第 2 モード（図 11(c)）の分布はタイヤを中心として左右に双極となっている．特に第 2 モードは左右の音響波の伝播が位7(a) Power ratio for POD modes at St = 25.6(b) 1st mode(c) 2nd mode

図 12: POD modes at St = 25.6 (1380[Hz]) are shown,
where POD is conducted on the ground. Top, middle
and bottom ﬁgures show the real part, amplitude andphase of each mode, respectively. (a) shows the power
ratio for each mode. (b) shows the ﬁrst mode, and (c)shows the second one.

Copyright c 2013 by JSFM

第 27 回数値流体力学シンポジウム
A10-2相差を持っている事から，St = 2.56 程顕著では無いものの，二重極のような音源を捉えている．
最後に St = 25.6 の場合には，第 1,2 モード共に良く似た分布となっている（図 12(b)(c)）．まず，タイヤ接地面前方と後方に指向性を持つ強い音響波が 2 つ捉えられている．また，タイヤ接地面側方の左右それぞれから斜め前方と斜め後方に指向性を持つ 4 つの音響波が現れている．ゆえに，St = 25.6 の場合には低周波で見られた単極子，二重極音源よりも複雑な多重極音源が形成されている可能性が高い．また図 12(b) 及び (c) のタイヤ後方で細かく振動するモードは，図 6 におけるタイヤ後方の細かい渦構造に対応する．
最後に，前節の最後に指摘したタイヤ後方の数値振動の音響波への影響を考察する．St = 2.56（図 10(b)(c)）の第 1,2 モードにはタイヤ後方にタイヤ軸方向に伸びる数値振動が現れるが，St = 12.8, 25.6 の高周波領域では見られない．この事は，タイヤ接地面前方（積分領域）のSt = 2.56 付近の低周波領域の音響波にタイヤ後方の数値振動が影響しうる一方で，今回の解析対象である高周波領域の音響波には殆ど影響が無い事を示唆している．

4.

結論

高解像度スキーム（コンパクト差分）を用いた LES に基づく直接解法により，回転するタイヤ周り（タイヤ径，回転速度で定義したレイノルズ数は 105 ）に発生する音響波を解析した．流体計算では IBM と重合格子法を併用し，数値振動を極力抑えつつタイヤ接地面付近の流体場も安定に解き進める事が出来るようになった．計算の結果，タイヤ表面近傍と後方には 3 次元性の強い細かい渦構造が形成される事が分かった．また音圧レベルの空間分布を可視化し，1000Hz 付近の高周波領域でタイヤ接地面前方を中心とする同心円状の特徴的な音響波を初めて捉える事が出来た．さらに周波数領域での POD を行い，低周波では接地面付近に単極子または双極子音源が存在する事を明らかにした．また高周波領域での接地面前方に拡がる音響波は，双極子以上の多重極音源を有する可能性がある事を示した．これらの知見はタイヤ表面のパターン（溝）が無い場合の結果であり，表面にパターンを付けた場合に発生する音響波との違い等を今後解析していきたい．

(2) Fujii, K., “Uniﬁed Zonal Method Based on theFortiﬁed Solution Algorithm,” Journal of Computational Physics, Vol. 118, 1995, pp. 92–108.
(3) Lele, S. K., “Compact Finite Diﬀerence Schemeswith Spectral-like Resolution.” Journal of Computational Physics, Vol. 103, No. 1, 1992, pp. 16–42.
(4) Bogey, C. and Bailly, C., “A family of low dispersive and low dissipative explicit schemes for ﬂow
and noise computations,” Journal of ComputationalPhysics, Vol. 194, 2004, pp. 194–214.
(5) Nishida, H. and Nonomura, T., “ADI-SGS Scheme
on Ideal Magnetohydrodynamics,” Journal of Computational Physics, Vol. 228, 2009, pp. 3182–3188.
(6) Ye, T., Mittal, R., Udaykumar, H. S., and Shyy,W., “An Accurate Cartesian Grid Method for Viscous Incompressible Flows with Complex Immersed
Boundaries,” Journal of Computational Physics,
Vol. 156, No. 9, 1999, pp. 209–240.
(7) Freund, J. B. and Colonius, T., “Turbulence andSound-ﬁeld POD Analysis of a Turbulent Jet,” International Journal of Aeroacoustics, Vol. 8, 2009,
pp. 337–354.
(8) Nonomura, T. and Fujii, K., “Computational Studyof Eﬀects of Near-wall Turbulent Structure on
Aeroacoustic Waves from a Supersonic Jet Impinging on an Inclined Plate,” AIAA 2011-2 , 2011.
(9) Suzuki, T., Bodony, D., Ryu, J., and Lele, S. K.,“Noise Sources of High-Mach-number Jets at
Low Frequencies Studied with a Phased-array
Approach Based on LES Database,” Annual
research briefs 2007, Center for Turbulence Research, NASA Ames and Stanford University, 2007,
http://www.stanford.edu/group/ctr/ResBriefs/ARB07.html.

謝辞
本研究の解析には，JAXA の SX9（JSS-SX9）及び東北大学サイバーサイエンスセンターの SX9 を使用した．
また，本研究の一部は JSPS 科研費 258793 の助成を受けた．ここに記して謝意を表する．
参考文献(1) Kaneda, I., Nonomura, T., Fujii, K., Ikeda, T., andKoishi, M., “Computational Analysis of Aeroacoustic Waves Induced by Rotating Tire,” Proceeding ofEuropean Conference on Computational Fluid Dynamics CFD 2010 , June 2010.
8Copyright c 2013 by JSFM