JP7052552B2 - 風切音解析装置及び風切音解析方法 - Google Patents

Description

本発明は、移動中の構造物に発生する風切音を解析する風切音解析装置及び風切音解析方法に関する。特に、構造物のガラス製の部分における高周波領域の表面音圧変動に起因して発生する風切音を解析する風切音解析装置及び風切音解析方法に関する。

車両を始めとする構造物の移動中に発生する風切音を低減するための技術が研究、開発されている。例えば、特許文献１には、車両モデルを走行させるシミュレーションにより当該車両モデル表面の任意の位置における「圧力変動の振幅」と「平均流速」とを取得し、これらの値に基づいて当該任意の位置における風切音の音源の強度を所望の周波数帯毎に算出する技術が開示されている。なお、風切音とは、構造物周囲の流体の流れによって発生する騒音であり、本明細書では、特に、構造物内部にまで到達する騒音を意味する。

特開２０１７－０６２７２７号公報

風切音の強度は、上述した構造物表面の圧力変動の振幅と平均流速のうち、特に圧力変動の振幅と強く相関している。ここで、移動中の構造物の表面に生じる圧力は、流体流動圧力（convective pressure）と音圧（acoustic pressure）とにより構成されている。流体流動圧力は流体流速で移動しそれほど遠方まで伝播しないのに対し、音圧は音速で移動し遠方まで伝播する。加えて、流体流動圧力は波長が比較的に短く振幅が比較的に大きいのに対し、音圧は波長が比較的に長く振幅が比較的に小さい。具体的には、音圧変動の振幅は、流体流動圧力変動の振幅の約１００分の１から１０００分の１である。このため、通常は、移動中の構造物の表面における圧力の大半は流体流動圧力が占めており、風切音の強度は流体流動圧力変動の振幅と強く相関する。

しかしながら、高周波領域（例えば、２ｋＨｚの周波数帯）では音圧の波長がガラスの振動モードの波長と略等しくなる。このため、構造物の筐体が金属製及び／又は樹脂製の部分だけではなくガラス製の部分を有する場合、高周波領域では当該ガラス製の部分の表面における音圧変動によりガラス製の部分が共鳴してその振動が増大し、風切音の強度が高くなる。この場合、当該風切音の強度は、流体流動圧力変動の振幅ではなく音圧変動の振幅と強く相関する。なお、以下では、構造物表面における、圧力変動の振幅、音圧変動の振幅及び流体流動圧力変動の振幅を、それぞれ「表面圧力変動」、「表面音圧変動」及び「表面流体流動圧力変動」とも称する。

表面圧力変動（即ち、表面音圧変動と表面流体流動圧力変動との和）は、構造物周囲の流れ場（典型的には、流速及び渦度）の影響を受けて変化する。このため、高周波領域における構造物のガラス製の部分の表面音圧変動に起因して発生する風切音を低減するためには、当該ガラス製の部分のうち表面音圧変動が比較的に大きい位置を特定し、構造物周囲の流れ場のうちどの部分（位置）の流れ場が当該特定された位置の表面音圧変動に強く影響を与えているのかを特定し、上記特定された位置の表面音圧変動が低減するように当該部分の流れ場を変更する（即ち、構造物の部品形状を変更する）必要がある。しかしながら、「表面音圧変動に大きく寄与している流れ場の部分を特定する方法」には統一された手法があるわけではなく、実際には技術者が自身のノウハウ及び知見に基づいて当該流れ場の部分を特定している。このため、技術者によって検討結果にばらつきが生じ、試行錯誤が多く発生するという問題があった。

本発明は、上記の問題を解決するためになされたものであり、その目的の一つは、車両を始めとする構造物のガラス製の部分における高周波領域の表面音圧変動に大きく寄与している構造物周囲の流れ場の位置を適切に特定することが可能な技術を提供することにある。

本発明の風切音解析装置は、移動中の構造物のガラス製の部分であるガラス部の表面に発生する風切音を解析する。
この風切音解析装置は、
構造物をモデル化した構造物モデル（２０）を移動させる非定常ＣＦＤシミュレーションを実行して、前記構造物モデル（２０）の周囲の流れ場のうち所定領域（２１）内の流れ場の所定時間における平均流速（Ｕ－）及び平均渦度（Ω－）を、当該所定領域（２１）内の節点である空間節点（ｚ）毎に算出するとともに、風切音の解析対象となる角周波数帯である着目角周波数帯（[θ_ｌ，θ_ｈ]）において、前記所定領域（２１）内における乱れ流速の振幅（ｕ～_θ）に基づく値を、前記空間節点（ｚ）毎に算出する非定常ＣＦＤ計算手段（ステップ８０２、８０４、８０６）と、
前記構造物モデル（２０）のガラス部（２２）表面の風切音の解析対象点（ｘ）における着目角周波数帯（[θ_ｌ，θ_ｈ]）の音圧変動の振幅である表面音圧変動（ｐａ_θ）に対する前記所定領域（２１）内の空間節点（ｚ）における流れ場の寄与度を示す指標である音圧ソース密度（ＡＰＤＳ）を、前記非定常ＣＦＤ計算手段（ステップ８０２、８０４、８０６）によって算出された前記平均流速（Ｕ－）と、前記平均渦度（Ω－）と、前記乱れ流速の振幅（ｕ～_θ）に基づく値と、に基づいて算出する音圧ソース密度算出手段（ステップ８０８）と、
を備える。

本発明の風切音解析装置（以下、「本発明装置」とも称する。）は、構造物モデルを移動させる非定常ＣＦＤシミュレーションにより算出された物理量（平均流速、平均渦度及び乱れ流速）に基づいて着目角周波数帯の音圧ソース密度を算出する。この音圧ソース密度は、「構造物モデルのガラス製の部分（ガラス部）表面の解析対象点における着目角周波数帯の表面音圧変動」に対する「所定領域内の空間節点における流れ場」の寄与度を示す指標である。ここで、所定領域とは、構造物モデルの周囲の流れ場の一部であり、解析対象点における表面音圧変動に影響を与える可能性のある領域として規定された領域である。このため、音圧ソース密度を指標として用いることにより、着目角周波数帯の表面音圧変動に対する流れ場の寄与度を空間節点毎に算出できる。従って、本発明装置によれば、所定領域内のどの位置における流れ場が解析対象点における構造物モデルのガラス部における高周波領域の表面音圧変動に大きく寄与しているかを適切に特定することが可能となる。なお、表面音圧変動は風切音と強く相関しており、表面音圧変動が大きいときは風切音が大きくなる。このため、表面音圧変動を解析することは、広義には、風切音を解析することと同義である。

本発明装置の一側面では、
前記構造物モデル（２０）のガラス部（２２）表面の前記解析対象点（ｘ）についての着目角周波数帯（[θ_ｌ，θ_ｈ]）の前記音圧ソース密度（ＡＰＤＳ）を前記所定領域（２１）で空間積分した値は、前記解析対象点（ｘ）における前記表面音圧変動の関数（ｐａ_θ）とその複素共役関数（ｐａ_θ ^＊）との積を着目角周波数帯（[θ_ｌ，θ_ｈ]）で積分した値の近似値である。

この構成によれば、音圧ソース密度の挙動は、構造物モデルのガラス部における表面音圧変動の挙動と精度良く一致する。即ち、音圧ソース密度は、表面音圧変動に対する流れ場の寄与度を表す指標として高い信頼性を有する。このため、構造物のガラス部における高周波領域の表面音圧変動に大きく寄与している構造物周囲の流れ場の位置を高い精度で特定することができる。

本発明装置の一側面は、
更に、前記音圧ソース密度（ＡＰＤＳ）を構成する複数のパラメータのうち、当該音圧ソース密度（ＡＰＤＳ）への寄与が相対的に大きいパラメータである原因パラメータを特定する原因パラメータ特定手段を備える。

この構成によれば、ＡＰＤＳが高い原因となっているパラメータを容易に把握できるため、構造物モデルの形状をより効率的に検討、変更することができる。

本発明装置の一側面では、
前記複数のパラメータは、前記平均流速（Ｕ－）、前記平均渦度（Ω－）及び前記乱れ流速（ｕ～_θ）である。

本発明装置の一側面は、
更に、前記音圧ソース密度（ＡＰＤＳ）が算出された複数の空間節点（ｚ）の中から、外部から入力される値を有する音圧ソース密度（ＡＰＤＳ）に対応している複数の空間節点（ｚ）を抽出し、当該抽出された複数の空間節点（ｚ）を画像処理して等値面（２８、３０、３２）を作成し、当該等値面（２８、３０、３２）を可視化する画像処理手段（ステップ８２０）を備える。

この構成によれば、ＡＰＤＳの等値面が表示されるため、作業者は、入力値を適宜設定することにより、所望のＡＰＤＳ値を有する等値面を視認することができる。この結果、所望のＡＰＤＳ値に対応している空間節点を効率的に選択することができ、構造物のガラス部における高周波領域の表面音圧変動に大きく寄与している構造物周囲の流れ場の位置を効率的に特定することができる。

本発明装置の一側面は、
更に、前記音圧ソース密度算出手段（ステップ８０８）によって算出された複数の音圧ソース密度（ＡＰＤＳ）の最大値に対応している空間節点（ｚ）を抽出する空間節点抽出手段を備える。

この構成によれば、最大のＡＰＤＳ値に対応している空間節点を作業者が選択する必要がなくなるため、構造物のガラス部における高周波領域の表面音圧変動に大きく寄与している構造物周囲の流れ場の位置をより効率的に特定することができる。なお、「音圧ソース密度の最大値」とは、広義には、「音圧ソース密度がほぼ最大である値」を意味する。

本発明装置の一側面では、
前記音圧ソース密度（ＡＰＤＳ）の計算式は、下記式によって定義される。

加えて、本明細書は、計算装置を用いて移動中の構造物のガラス製の部分であるガラス部の表面に発生する風切音を解析する新規な風切音解析方法を開示する。
この風切音解析方法は、
前記計算装置が備える非定常ＣＦＤ計算手段が、構造物をモデル化した構造物モデル（２０）を移動させる非定常ＣＦＤシミュレーションを実行して、前記構造物モデル（２０）の周囲の流れ場のうち所定領域（２１）内の流れ場の流速の時刻歴データ（ｕ（ｚ，ｔ））及び渦度の時刻歴データ（ω（ｚ，ｔ））を、当該所定領域（２１）内の節点である空間節点（ｚ）毎に所定時間に亘ってそれぞれ算出するステップ（ステップ８０２）と、
前記非定常ＣＦＤ計算手段が、前記算出された流速の時刻歴データ（ｕ（ｚ，ｔ））及び渦度の時刻歴データ（ω（ｚ，ｔ））をそれぞれ平均化処理して、平均流速（Ｕ－（ｚ））及び平均渦度（Ω－（ｚ））を前記空間節点（ｚ）毎に算出するステップ（ステップ８０４）と、
前記非定常ＣＦＤ計算手段が、前記算出された流速の時刻歴データ（ｕ（ｚ，ｔ））を高速フーリエ変換処理して、乱れ流速（ｕ～ _θ （ｚ））の振幅の関数とその複素共役関数との積を風切音の解析対象となる角周波数帯である着目角周波数帯（[θ_ｌ，θ_ｈ]）で積分した値である乱れ流速（ｕ～ _θ （ｚ））の振幅の自己相関関数を前記空間節点（ｚ）毎に算出するステップ（ステップ８０６）と、
前記計算装置が備える音圧ソース密度算出手段が、前記構造物モデル（２０）のガラス部（２２）表面の風切音の解析対象点（ｘ）における着目角周波数帯（[θ_ｌ，θ_ｈ]）の音圧変動の振幅である表面音圧変動（ｐａ_θ）に対する前記所定領域（２１）内の空間節点（ｚ）における流れ場の寄与度を示す指標である音圧ソース密度（ＡＰＤＳ）を、予めＲＯＭに格納されている音圧ソース密度（ＡＰＤＳ）の計算式であって、平均流速（Ｕ－（ｚ））と、平均渦度（Ω－（ｚ））と、乱れ流速（ｕ～ _θ （ｚ））の振幅の自己相関関数と、をパラメータとして含む計算式、に、前記非定常ＣＦＤ計算手段によって前記空間節点（ｚ）毎に算出された前記平均流速（Ｕ－）と、前記平均渦度（Ω－）と、前記乱れ流速の振幅（ｕ～_θ）の自己相関関数と、のうち、対応する空間節点（ｚ）における平均流速（Ｕ－（ｚ））と、平均渦度（Ω－（ｚ））と、乱れ流速（ｕ～ _θ （ｚ））の振幅の自己相関関数と、を代入することにより算出するステップ（ステップ８０８）と、
を備える。

この風切音解析方法によれば、所定領域内のどの位置における流れ場が解析対象点における構造物モデルのガラス部における高周波領域の表面音圧変動に大きく寄与しているかを適切に特定することが可能となる。

なお、上記説明においては、発明の理解を助けるために、実施形態に対応する発明の構成に対して、実施形態で用いた符号を括弧書きで添えているが、発明の各構成要件は前記符号によって規定される実施形態に限定されるものではない。

本発明の実施形態に係る風切音解析装置（以下、「本実施装置」又は「実施装置」と称する。）が有する計算装置を示すブロック図である。 非定常ＣＦＤシミュレーションの対象となる車両モデル及びフロードメインの模式図である。 車両モデルの解析表面（右側フロントサイドガラス）における表面音圧変動の分布図である。 表面音圧変動（風切音）の解析対象となる観測点の位置を示す図である。 着目評価点の選択方法を説明するために用いる図であり、ＡＰＤＳの参照値が基準値である等値面を示す図である。 ＡＰＤＳの参照値が基準値の１６倍である等値面を示す図である。 ＡＰＤＳの参照値が基準値の２５６倍である等値面を示す図である。 原因パラメータ（平均流速）の分布図が作成される平面（着目評価点通過断面）の位置を示す図である。 着目評価点通過断面における右側サイドミラー周囲の平均流速の分布図である。 ＡＰＤＳ計算式の導出方法を説明するために用いる図であり、剛体モデル及びフロードメインの模式図である。 剛体モデル及びフロードメインの模式図である。 風切音解析方法の手順を示すフローチャートである（その１）。 風切音解析方法の手順を示すフローチャートである（その２）。 車両モデルの形状を変更したときの表面音圧変動（風切音）の改善効果を説明するために用いる図であり、形状変更前の車両モデルの正面図の一部を示す。 形状変更前の右側サイドミラーの平面図である。 形状変更後の車両モデルの正面図の一部である。 形状変更後の右側サイドミラーの平面図である。 形状変更前の着目評価点通過断面における右側サイドミラー周囲の平均流速の分布図である。 形状変更後の着目評価点通過断面における右側サイドミラー周囲の平均流速の分布図である。 形状変更後の車両モデルの右側フロントサイドガラスにおける表面音圧変動の分布図である。 指標としてのＡＰＤＳの信頼性及び汎用性を説明するために用いる図であり、車両モデルの模式図を示す。 フォアステップモデルの模式図である。 車両モデルの右側フロントサイドガラスにおける、ＡＰＤＳを用いて算出される予測表面音圧変動の分布図である。 車両モデルの右側フロントサイドガラスにおける、ソフトウエアの非定常ＣＦＤシミュレーションにより算出される表面音圧変動の分布図である。 フォアステップモデルのステップの上面における、ＡＰＤＳを用いて算出される予測表面音圧変動の分布図である。 フォアステップモデルのステップの上面における、ソフトウエアの非定常ＣＦＤシミュレーションにより算出される表面音圧変動の分布図である。 表面音圧変動に対する予測表面音圧変動の誤差を車両モデル及びフォアステップモデルについて示したグラフである。

＜本実施装置の概要＞
まず、実施形態に係る風切音解析装置（以下、「本実施装置」とも称する。）の概要を説明する。本実施装置は、様々な構造物の移動時に発生する風切音を解析するが、本実施形態では、構造物のうち車両を例に挙げて説明する。車両走行時に発生する風切音（車内の乗員に聞こえる騒音）は、車両表面の複数の位置における種々の周波数帯の圧力変動の振幅（即ち、表面圧力変動であり、別言すれば、表面音圧変動と表面流体流動圧力変動との和）に起因した複合的な騒音である。これらの複合的な風切音のうち、高周波領域の風切音は、車両のガラス製の部分（例えば、フロントサイドガラス、リアサイドガラス、フロントガラス及びリアガラス。以下では、「ガラス部」とも称する。）における高周波領域の表面音圧変動と強く相関している。

この表面音圧変動は、車両周囲の流れ場（流速及び渦度）の影響を受けて変化する。本実施装置には、車両周囲の流れ場と表面音圧変動との関係を定量的に示す指標（別言すれば、「車両周囲の任意の位置における流れ場（の状態）」が「車両のガラス部の任意の位置における表面音圧変動」にどれくらい寄与しているかを示す指標）の計算式が予め格納されている。指標の値が大きいほど寄与度は高くなる。この指標は、車両周囲の任意の位置における流れ場の平均流速、乱れ流速及び平均渦度等の物理量をパラメータとして有する。加えて、車両のガラス部の任意の位置における表面音圧変動は、当該指標を所定領域内の流れ場で空間積分した値に基づいて算出（予測）することができる。このため、本実施装置は、非定常ＣＦＤ（数値流体力学）計算により車両モデル周囲のうち所定領域内の流れ場の平均流速、乱れ流速及び平均渦度を所定領域内の位置毎に算出し、これらの値に基づいて、上記指標を所定領域内の位置毎に算出する。そして、当該指標に基づいて、走行中の車両モデルの表面音圧変動を、車両モデルのガラス部表面の位置毎、及び、周波数帯毎、に算出（予測）する。作業者は、当該算出結果に基づいて、表面音圧変動が比較的に大きい位置（別言すれば、大きな風切音の発生原因となっている位置）を選択する。なお、車両モデルは「構造物モデル」の一例に相当する。

作業者は、算出結果に基づいて、上記選択された位置についての指標であって、所定領域内の位置毎に算出されている複数の指標の中から、指標の値が比較的に大きい流れ場の位置（別言すれば、上記選択された位置における表面音圧変動に対する寄与度が比較的に高い位置）を選択し、当該選択された流れ場の位置における指標を構成する各物理量（即ち、平均流速、乱れ流速及び平均渦度）のうち、指標の値が大きい原因となっている物理量を特定する。作業者は、当該特定された物理量の値に基づいて、上記選択された流れ場の位置における指標値が低減するように車両モデルの形状を検討、変更する。この構成によれば、指標値が低減するように車両モデルの形状を変更することにより、「当該指標値を有する位置の流れ場の、上記選択された位置における表面音圧変動への寄与度」を低減することができる。このため、当該選択された位置における表面音圧変動を低減でき、結果として、風切音の低減を実現することができる。

＜本実施装置の具体的な構成＞
以下、図面を参照して本実施装置について具体的に説明する。図１に示すように、本実施装置は計算装置１０を備える。計算装置１０は、入力部１２、演算部１４及び出力部１６を有する。出力部１６は、作業者が視認可能な位置に表示画面１８を有する。入力部１２は、風切音の解析対象となる車両の三次元モデルのデータ、解析範囲を示すフロードメインのデータ、及び、非定常ＣＦＤシミュレーションに必要なデータ等を入力する。加えて、入力部１２は、処理の途中で作業者により選択される座標のデータ、指標の参照値及び原因パラメータ等を入力する（後述）。

図２Ａは、車両の三次元モデル２０（以下、単に「車両モデル２０」と称する。）及びフロードメイン２１を模式的に示す図である。図２Ａに示すように、車両モデル２０の空間座標系はｅ１軸、ｅ２軸及びｅ３軸からなり、当該座標系は、車両モデル２０の前後方向における前方が＋ｅ１方向に一致し、且つ、左右方向における左方向が＋ｅ２方向に一致するように設定される。加えて、原点は、車両モデル２０の基準点（例えば、重心）から所定の方向に所定の距離（ゼロ値を含む）だけ離れた位置に設定される。

フロードメイン２１は、車両モデル２０の周囲の空間のうち、非定常ＣＦＤ計算が実施される範囲を示す領域である。本実施形態では、走行中の車両モデル２０の右側フロントサイドガラス２２における表面音圧変動に起因して発生する風切音の解析が行われる。このため、フロードメイン２１は、「車両モデル２０の周囲の流れ場のうち、右側フロントサイドガラス２２における表面音圧変動に影響を与える可能性のある部分の流れ場を含んだ領域」として設定される。具体的には、フロードメイン２１は、右側フロントサイドガラス２２、右側サイドミラー２３、右側フロントピラー及びフロントガラスの右半分等を含む略直方体形状の空間から車両モデル２０と重複している部分がくりぬかれた形状を有する。フロードメイン２１は、車両モデル２０のうちどの部分を解析対象とするかによって自由に設定可能である。なお、フロードメイン２１は、「所定領域」の一例に相当する。

入力部１２にて入力されたこれらのデータは、演算部１４のＲＡＭ（後述）に格納される。なお、以下では、非定常ＣＦＤ計算を単に「ＣＦＤ計算」とも称する。

図１に戻って説明を続ける。演算部１４は、ＣＰＵ、ＲＯＭ、ＲＡＭ等からなるマイクロコンピュータを主要構成部品として有する。演算部１４は、各種インターフェースを有しており、このインターフェースを介して入力部１２及び出力部１６と信号を入出力可能に接続されている。

演算部１４は、ＲＡＭに格納されている上述したデータを用いて車両モデル２０を走行させる非定常ＣＦＤシミュレーションを実行する。ＣＦＤ計算処理は、時刻歴データの算出処理、平均化処理及び高速フーリエ変換処理を含む。この計算処理は、フロードメイン２１によって区画された範囲で行われる。以下、具体的に説明する。

（ＣＦＤ計算処理）
１．時刻歴データの算出処理
演算部１４は、フロードメイン２１内の各空間節点ｚについて、流れ場の流速の時刻歴データｕ（ｚ，ｔ）及び渦度の時刻歴データω（ｚ，ｔ）を所定時間に亘ってそれぞれ算出する。以下では、空間節点ｚを単に「節点ｚ」と称する。なお、フロードメイン２１内の節点ｚの位置は、節点ｚが車両モデル２０の表面上に存在することがないように（即ち、節点ｚが車両モデル２０の表面から離間するように）予め設定されている。演算部１４は、流速及び渦度の時刻歴データｕ（ｚ，ｔ）及びω（ｚ，ｔ）を、節点ｚの座標及び時刻ｔに関連付けてそのＲＡＭに格納する。

２．平均化処理
演算部１４は、ＲＡＭに格納されている流速及び渦度の時刻歴データｕ（ｚ，ｔ）及びω（ｚ，ｔ）をそれぞれ平均化処理して、フロードメイン２１内の各節点ｚについて、平均流速Ｕ－（ｚ）及び平均渦度Ω－（ｚ）を算出する。なお、Ｕ－（ｚ）及びΩ－（ｚ）の「－」は何れも平均を表す表記であり、後述する式中の文字Ｕ及びΩの上に付されたバーに替わるものである。演算部１４は、平均流速Ｕ－（ｚ）及び平均渦度Ω－（ｚ）を、節点ｚの座標に関連付けてそのＲＡＭに格納する。平均流速Ｕ－（ｚ）及び平均渦度Ω－（ｚ）は、指標を算出する際に用いられる（後述）。

３．高速フーリエ変換処理
演算部１４は、ＲＡＭに格納されている流速の時刻歴データｕ（ｚ，ｔ）を高速フーリエ変換処理して、フロードメイン２１内の各節点ｚについて、周波数（角周波数）毎に乱れ流速（厳密には、乱れ流速の振幅）ｕ～_θ（ｚ）の自己相関関数（後述）を算出する（θ：角周波数。θ＝２πｆ）。なお、ｕ～_θ（ｚ）の「～」は、流速ｕを平均化処理したときの乱れ成分を表す表記であり、後述する式中の文字ｕの上に付されたチルダに替わるものである。なお、乱れ流速の振幅ｕ～_θ（ｚ）の自己相関関数が「乱れ流速の振幅に基づく値」の一例に相当する。

自己相関関数は、風切音の解析対象となる周波数帯において算出される（以下、当該周波数帯を「着目周波数帯」とも称する。）。着目周波数帯は作業者によって選択可能であり、例えば、中心周波数ｆ_ｍが５００Ｈｚ、１ｋＨｚ、２ｋＨｚ又は４ｋＨｚの周波数帯が選択され得る。本実施装置は高周波領域の表面音圧変動に起因した風切音を解析するため、典型的には２ｋＨｚ以上の周波数帯が選択され得る。例えば、着目周波数帯として２ｋＨｚ（下限周波数ｆ_ｌ＝１４２０Ｈｚ、上限周波数ｆ_ｈ＝２８４０Ｈｚ）が選択された場合、乱れ流速ｕ～_θ（ｚ）の自己相関関数は、１４２０Ｈｚから２８４０Ｈｚの範囲内で、フロードメイン２１内の各節点ｚについて算出される。

なお、任意の関数Ｆ_θの自己相関関数は、下記式（１）によって定義される。

ここで、上記式（１）中のθ_ｌ及びθ_ｈは、着目角周波数帯の下限角周波数及び上限角周波数をそれぞれ表し、Ｆ_θ ^＊は、関数Ｆ_θの複素共役を表す。即ち、本明細書では、自己相関関数は、「任意の関数Ｆ_θとその複素共役関数Ｆ_θ ^＊との積（後述）を着目角周波数帯において積分した値」として定義される。

演算部１４は、乱れ流速ｕ～_θ（ｚ）の自己相関関数を、節点ｚの座標に関連付けてそのＲＡＭに格納する。乱れ流速ｕ～_θ（ｚ）の自己相関関数は、指標を算出する際に用いられる（後述）。演算部１４は、以上の処理をＣＦＤ計算として実行する。

（ＡＰＤＳの算出）
車両モデル２０の右側フロントサイドガラス２２の表面を「解析表面」と規定すると、演算部１４のＲＯＭには、フロードメイン２１内の流れ場と、解析表面における着目周波数帯の表面音圧変動との関係を定量的に示す指標の計算式が予め格納されている。この指標は、別言すれば、「フロードメイン２１内の任意の節点ｚにおける流れ場（の状態）」が「解析表面上の任意の表面節点ｘ（以下、単に「節点ｘ」と称する。）における着目周波数帯の表面音圧変動ｐａ_θ（ｘ）」にどれくらい寄与しているかを示す指標ということもできる。解析表面上の任意の基準節点から数えてｉ番目の節点ｘをｘ_ｉ（ｉ：１～ｍ）と規定し、フロードメイン２１内の任意の基準節点から数えてｊ番目の節点ｚをｚ_ｊ（ｊ：１～ｎ）と規定すると、ｋ番目の節点ｘ_ｋ（ｋ：1～ｍ）についての指標は、フロードメイン２１内の全節点ｚ_ｊにおいてそれぞれ算出される。

ここで、節点ｘ_ｋについての着目周波数帯における指標をフロードメイン２１内の全節点ｚ_ｊで空間積分した値は、「節点ｘ_ｋにおける着目周波数帯の表面音圧変動ｐａ_θ（ｘ_ｋ）の自己相関関数」の近似値である（後述）。このため、各節点ｚ_ｊにおける着目周波数帯の指標の値は、フロードメイン２１の単位体積当たりの「節点ｘ_ｋにおける着目周波数帯の表面音圧変動ｐａ_θ（ｘ_ｋ）の自己相関関数」の近似値を表していると解釈することができる。従って、以下では、各節点ｚ_ｊにおける着目周波数帯の指標を「ＡＰＤＳ（Acoustic Pressure Density Source（音圧ソース密度））」とも称する。着目周波数帯のＡＰＤＳ（ｘ_ｋ，ｚ_ｊ）は、「フロードメイン２１内の節点ｚ_ｊにおける流れ場」の「解析表面上の節点ｘ_ｋにおける着目周波数帯の表面音圧変動ｐａ_θ（ｘ_ｋ）」への寄与度を示す指標であり、その値が大きいほど寄与度が高いことを示す。ＡＰＤＳの計算式の表記及び導出方法については後述する。なお、以下では、ある節点ｘ_ｋを他の節点ｘ_ｉと区別する必要がない場合は、「ｋ」及び「ｉ」の表記を省略する。節点ｚ_ｉについても同様である。加えて、以下では、着目周波数帯の表面音圧変動を単に「表面音圧変動」とも称し、着目周波数帯のＡＰＤＳを単に「ＡＰＤＳ」とも称する。

演算部１４は、ＣＦＤ計算を終了すると、ＲＯＭに格納されているＡＰＤＳの計算式に基づいて、節点ｘ_ｋについてのＡＰＤＳをフロードメイン２１内の各節点ｚ_ｊにおいて算出する処理を、解析表面上の全節点ｘ_ｉについて行う。即ち、演算部１４は、節点ｘ_ｋについてのＡＰＤＳをｎ個（節点ｚの個数）算出するとともに、全部でｍｎ個のＡＰＤＳを算出する。ここで、ある節点ｘ_ｋについてのＡＰＤＳ（ｘ_ｋ，ｚ_ｊ）の値は節点ｚ_ｊ毎に評価されるため、以下では、節点ｚを「評価点ｚ」とも称する。演算部１４は、ＡＰＤＳ（ｘ，ｚ）の値を節点ｘの座標及び評価点ｚの座標に関連付けてそのＲＡＭに格納する。ＡＰＤＳは、着目周波数帯の中心周波数ｆ_ｍにおける予測表面音圧変動ｐａ_ｐθ（ｘ）（後述）の分布図を作成する際に用いられるとともに、ＡＰＤＳの等値面を作成する際に用いられる（後述）。

（予測表面音圧変動の算出及び分布図表示指令）
上述したように、「任意の節点ｘについてのＡＰＤＳをフロードメイン２１で空間積分した値」は、「当該節点ｘにおける表面音圧変動ｐａ_θ（ｘ）の自己相関関数」の近似値である。このため、演算部１４は、ＲＡＭに格納されているｍｎ個のＡＰＤＳ（ｘ，ｚ）のうち、節点ｘ_ｋについてのｎ個のＡＰＤＳ（ｘ_ｋ，ｚ_ｊ）を抽出し、抽出されたＡＰＤＳ（ｘ_ｋ，ｚ_ｊ）を足し合わせる（即ち、ＡＰＤＳ（ｘ_ｋ，ｚ_ｊ）をフロードメイン２１で空間積分する）ことにより、節点ｘ_ｋにおける表面音圧変動ｐａ_θ（ｘ_ｋ）の自己相関関数の近似値を算出する。そして、この近似値に基づいて、節点ｘ_ｋにおける表面音圧変動ｐａ_θ（ｘ_ｋ）を算出（予測）する。以下では、このようにして算出（予測）された表面音圧変動を、「予測表面音圧変動ｐａ_ｐθ」と称する。演算部１４は、この処理を、解析表面上の全節点ｘ_ｉについて行う。即ち、演算部１４は、ｍ個（節点ｘの個数）の予測表面音圧変動ｐａ_ｐθ（ｘ）を算出する。演算部１４は、予測表面音圧変動ｐａ_ｐθ（ｘ）の値を節点ｘの座標に関連付けてそのＲＡＭに格納する。

演算部１４は、ＲＡＭに格納されている予測表面音圧変動ｐａ_ｐθ（ｘ）にその大きさに応じた色データを対応づけたデータを各節点ｘについて作成し、作成したデータを表示画面１８に画像表示させる表示指令を出力部１６に送信する。予測表面音圧変動ｐａ_ｐθ（ｘ）が大きいほど濃色の色データが対応づけられる。出力部１６は、当該表示指令を受信すると、その表示画面１８に当該データを画像表示する。

図３Ａは、当該表示指令に基づいて表示画面１８に画像表示された解析表面（右側フロントサイドガラス２２）における予測表面音圧変動ｐａ_ｐθ（ｘ）の分布図の一例である。着目周波数帯は２ｋＨｚ（即ち、ｆ_ｍ＝２ｋＨｚ）に設定されている。作業者は、この分布図に基づいて予測表面音圧変動ｐａ_ｐθ（ｘ）が比較的に大きい節点ｘを選択する。図３Ａによると、着目周波数帯が２ｋＨｚのときは領域２４及び領域２６において予測表面音圧変動ｐａ_ｐθ（ｘ）が顕著であるため、作業者は、領域２４、２６内の節点ｘの中から任意の節点を選択する。この例では、領域２４内の節点ｘが選択される。作業者は、以後、当該節点における表面音圧変動の低減を目的として車両モデル２０の形状を変更して当該節点における予測表面音圧変動の挙動を観測する。このため、以下では、上記選択された節点を「観測点ｘ」とも称する（図３Ｂ参照）。観測点ｘは、「解析表面のうち、着目周波数帯における大きな風切音の発生原因となっている位置」を表す。作業者により観測点ｘが選択されると、入力部１２は、当該観測点ｘの座標データを入力する。入力された座標データは演算部１４のＲＡＭに格納される。なお、観測点ｘは、「解析対象点」の一例に相当する。

（ＡＰＤＳ等値面表示指令）
演算部１４は、観測点ｘが選択されると、ＲＡＭに格納されているｍｎ個のＡＰＤＳ（ｘ，ｚ）の中から、観測点ｘについてのｎ個のＡＰＤＳ（ｘ，ｚ）を抽出する。そして、抽出されたＡＰＤＳ（ｘ，ｚ）の中から、作業者により設定された任意の値（以下、「参照値」とも称する。後述。）を有するＡＰＤＳ（ｘ，ｚ）に関連付けられた評価点ｚの座標を更に抽出する。そして、抽出された複数の評価点ｚ同士を接続した等値面を表す等値面データを作成し、作成した等値面データを表示画面１８に画像表示させる表示指令を出力部１６に送信する。出力部１６は、当該表示指令を受信すると、その表示画面１８に「参照値を有するＡＰＤＳの等値面データ」を画像表示する。即ち、演算部１４は、観測点ｘが選択されると、ＡＰＤＳの参照値の入力を促すメッセージを表示画面１８に表示させる。これにより、作業者により任意の参照値が入力される。演算部１４は、上記表示指令の送信処理を、作業者により参照値が入力される毎に繰り返し行う。

図４Ａ～図４Ｃは、参照値が互いに異なるＡＰＤＳの等値面を示す。図４Ａは、参照値が基準値（本例では６５．５）である等値面２８を示し、図４Ｂは、参照値が基準値の１６倍（１０４８．６）である等値面３０を示し、図４Ｃは、参照値が基準値の２５６倍（６７１０８．９）である等値面３２を示す。着目周波数帯は何れも２ｋＨｚに設定されている。図４Ａ～図４Ｃに示すように、等値面２８、３０、３２は何れも閉曲面である。加えて、各等値面２８、３０、３２は、複数の独立した等値面を有している。ＡＰＤＳ（ｘ，ｚ）は、観測点ｘ（図３Ｂ参照）から評価点ｚまでの距離ｒが減少するほど増加する性質を有する。即ち、評価点ｚが観測点ｘに接近するほどＡＰＤＳの値が増加する。このため、図４Ａ～図４Ｃに示すように、参照値が増加するほど等値面２８、３０、３２の表面積が減少している。加えて、等値面３０は等値面２８の内部に含まれており、等値面３２は等値面３０の内部に含まれている。即ち、ある参照値を有する等値面の内部には、当該参照値より大きい参照値を有する等値面が包含されている。このため、参照値が徐々に増加すると等値面の表面積が徐々に減少していき、最終的には数個の等値面が互いに離間して存在することになる。

図４Ｃにその一例を示す。図４Ｃに示すように、等値面３２は３個の等値面３２ａ、３２ｂ及び３２ｃを有している（厳密には、右側フロントサイドガラス２２上に複数の線上の等値面が形成されているが、本例ではこれらの等値面は解析対象から除外している。）。図３Ｂを参照すると、これら３個の等値面３２ａ～３２ｃの中では、等値面３２ｂが観測点ｘに最も近く（即ち、距離ｒが最も短い）、等値面３２ａが観測点ｘに最も遠い（即ち、距離ｒが最も長い）。このため、等値面３２ｂのＡＰＤＳは、等値面３２ａ又は３２ｃのＡＰＤＳよりも距離ｒの影響を強く受けている（即ち、距離ｒのＡＰＤＳ値への寄与度が高い）と言える。ここで、上述したように、ＡＰＤＳ（ｘ，ｚ）は、パラメータとして以下の物理量、即ち、平均流速Ｕ－（ｚ）、乱れ流速ｕ～_θ（ｚ）、及び、平均渦度Ω－（ｚ）を有する。ＡＰＤＳの値は、これらの物理量が大きいほど大きくなる。このため、「等値面３２ｂよりも距離ｒが長い等値面３２ａ又は３２ｃ」が等値面３２ｂと同一のＡＰＤＳ値を有しているのは、等値面３２ａ又は３２ｃのＡＰＤＳが等値面３２ｂのＡＰＤＳよりも上記物理量（Ｕ－（ｚ）、ｕ～_θ（ｚ）、Ω－（ｚ））の少なくとも１つの影響を強く受けている（即ち、上記物理量のＡＰＤＳ値への寄与度が高い）ためであると解釈することができる。加えて、距離ｒが最も長い等値面（本例では等値面３２ａ）が上記物理量の少なくとも１つの影響を最も強く受けている（即ち、上記物理量のＡＰＤＳ値への寄与度が最も高い）と解釈することができる。

従って、作業者は、距離ｒが最も長い等値面３２ａ上の任意の評価点ｚを、ＡＰＤＳ値が比較的に大きい（ほぼ最大である）評価点ｚとして選択する。以下では、選択された評価点ｚを「着目評価点ｚ」とも称する（図５Ａ参照）。着目評価点ｚは、「フロードメイン２１内の流れ場のうち、観測点ｘにおける予測表面音圧変動に対する寄与度が比較的に高い（ほぼ最大である）位置」を表す。作業者により着目評価点ｚが選択されると、入力部１２は、当該着目評価点ｚの座標データを入力する。入力された座標データは演算部１４のＲＡＭに格納される。

（ＡＰＤＳパラメータ表示指令）
演算部１４は、ＲＡＭに格納されている着目評価点ｚにおけるＡＰＤＳ（ｘ，ｚ）を構成している複数のパラメータ（物理量）の中から平均流速Ｕ－（ｚ）、乱れ流速ｕ～_θ（ｚ）、及び、平均渦度Ω－（ｚ）を選択し、これらの数値データを表示画面１８に表示させる表示指令を出力部１６に送信する。出力部１６は、当該表示指令を受信すると、その表示画面１８に選択されたパラメータの数値データを表示する。作業者は、表示された数値データに基づいて、これらのパラメータの中でＡＰＤＳ値への寄与が相対的に大きいパラメータ（即ち、ＡＰＤＳ値が大きい原因となっているパラメータ）を特定する。このとき、必要に応じて他の評価点ｚにおけるＡＰＤＳを構成するパラメータの数値データが表示されてもよい。作業者は、これらの数値データを比較しながらパラメータを特定してもよい。作業者によりパラメータが特定されると、入力部１２は、当該パラメータを入力する。入力されたパラメータは演算部１４のＲＡＭに格納される。以下では、このパラメータを「原因パラメータ」とも称する。

（物理量分布図表示指令）
演算部１４は、ＲＡＭに格納されている評価点ｚ群の座標データから、「着目評価点ｚを通りｅ１ｅ２平面（図２Ａ参照）に平行な平面（図５Ａの破線Ｌ参照）」上に位置する評価点ｚ群を抽出する。そして、抽出された評価点ｚ群に関連付けられた原因パラメータの分布図を表す分布図データを作成し、作成した分布図データを表示画面１８に画像表示させる表示指令を出力部１６に送信する。出力部１６は、当該表示指令を受信すると、その表示画面１８に当該データを画像表示する。なお、以下では、「着目評価点ｚを通りｅ１ｅ２平面に平行な平面」を「着目評価点通過断面」とも称する。

図５Ｂは、原因パラメータとして平均流速Ｕ－（ｚ）が特定された場合における平均流速Ｕ－（ｚ） [m/sec]の分布図を示す。平均流速Ｕ－（ｚ）が大きいほど濃色で表示されている。図５Ａに示すように、この例では着目評価点ｚが右側フロントサイドガラス２２の前方と右側サイドミラー２３との間に位置しているため、分布図の範囲は、右側フロントサイドガラス２２の前方及び右側サイドミラー２３を含む範囲とされている。なお、右側フロントサイドガラス２２及び右側サイドミラー２３の位置には、車両モデル２０を平面視したときの右側フロントサイドガラス２２及び右側サイドミラー２３の画像がはめ込み表示されている。着目評価点ｚは、図５Ｂの領域３４に含まれている。

作業者は、分布図に基づいて着目評価点ｚにおいて原因パラメータ（本例では平均流速Ｕ－（ｚ））が大きくなっている原因を分析する。図５Ｂの例では、右側サイドミラー２３の車体側近傍に位置している領域３４の平均流速Ｕ－（ｚ）が顕著である。一般に、流れ場の流速は空気抵抗が小さいほど大きい。このため、領域３４では車両走行時における空気抵抗が小さいと考えられる。従って、作業者は、領域３４の空気抵抗を増大できるような右側サイドミラー２３の形状を検討してその形状を変更する。形状変更後の予測表面音圧変動ｐａ_ｐθ（ｘ）の具体的な改善効果については後述する。

出力部１６は、演算部１４からの表示指令を受信して表示画面１８に各種データを表示する。詳細は上述したとおりである。以上が、本実施装置の具体的な構成についての
説明である。

＜ＡＰＤＳ計算式の導出＞
続いて、上述したＡＰＤＳの計算式の導出方法を剛体モデルを用いて説明する。ＡＰＤＳは、公知のパウエル（Powell）式に基づいて算出される「剛体モデルにおける表面圧力変動（即ち、表面音圧変動と表面流体流動圧力変動との和）とその周囲の流れ場との関係式」を時間ドメインから周波数ドメインに変換した式に基づいて定義される。但し、上記パウエル式は、ＡＰＤＳ計算式の導出に適した形式で表現される。このため、以下では、まず、図６を参照して剛体モデル及びフロードメインについて説明し、続いて、ＡＰＤＳ計算式の導出に適した形式のパウエル式の導出について説明する。次いで、時間ドメインにおける表面圧力変動と流れ場との関係式の導出について説明し、その後、ＡＰＤＳ計算式の導出について説明する。

（剛体モデル及びフロードメイン）
図６は、剛体モデル３６及びフロードメイン３８の模式図を示す。剛体モデル３６（以下、単に「剛体３６」とも称する。）は、その表面Ｓ上に観測点ｘを有する。フロードメイン３８は、剛体３６の周囲の流れ場のうち、観測点ｘにおける表面圧力変動に影響を与える可能性のある部分の流れ場である。フロードメイン３８は、その内部に複数の評価点ｚを有する。ＡＰＤＳ（ｘ，ｚ）は、観測点ｘにおける表面音圧変動に対する評価点ｚにおける流れ場の寄与度を表す指標である。なお、図６中のその他の文字については後述する。

（パウエル式の導出）
下記式（２）は、体積力がゼロの場合のナビエ－ストークス方程式（Navier-Stokes equation）を示す。

ここで、ρは流れ場の密度であり、ｕは流れ場の流速（ベクトル量）であり、Ｄ／Ｄｔは物質微分であり、ｐは剛体３６の表面における圧力（圧力変動、表面圧力変動）であり、ηは流れ場の粘性係数である。

上記式（２）の物質微分を偏微分を用いて書き直し、その後、以下の関係式；

を用いて式変形すると、下記式（３）が得られる。なお、上記の関係式中のＡは任意のベクトルを表す。

ここで、νは流れ場の動粘性係数であり、ν＝η／ρである。

続いて、上記式（３）を以下の関係式；

を用いて式変形すると、下記式（４）が得られる。なお、ωは流れ場の渦度（ベクトル量）であり、Ｂはエンタルピーである。

本実施形態の流れ場は非圧縮性流体又は微圧縮性流体であるが、この場合、以下の関係式；

が成立するため、上記式（４）は下記式（５）によって近似できる。

上記式（５）の両辺に流れ場の密度ρを乗算して発散をとると、下記式（６）が得られる。

なお、上記式（６）中の「×」は乗算記号を表す。本明細書では基本的に乗算記号は省略して記載するが、明示したほうが理解し易い場合は例外的に記載する。

ここで、流れ場の流速ｕと流れ場の密度ρとの間には、下記式（７）に示す関係式が成立する。

「流れ場の密度ρ」と「流れ場の流速ｕの時間についての偏微分」との積の発散（即ち、下記式（８）の左辺）を上記式（７）を用いて式変形すると、下記式（８）の右辺のようになる。

ここで、ｃは流れ場の音速である。表面圧力変動ｐ、流れ場の密度ρ、及び、流れ場の音速ｃとの間には以下の関係が成立する。

ところで、本実施形態の流れ場は高レイノルズ数流れ（即ち、流速が比較的に大きい流れであり、例えば、１００km/h～１４０km/h。）である。高レイノルズ数流れでは、以下の関係式；

が成立する。なお、レイノルズ数Ｒｅは、Ｒｅ＝ｕＬ／νとして定義され、Ｌは流れ場の代表長さを表す。上記式（８）の右辺の括弧内の項を以下に示すように式変形して上記関係式を適用すると、式（８）の右辺の括弧内の項は下記式（９）によって近似される。

上記式（６）を上記式（８）及び式（９）を用いて式変形（近似）すると、下記式（１０）が得られる。

ここで、本実施形態の流れ場は低マッハ数流れ（即ち、音速に対する流体の流速の比が比較的に低い流れ場）である。低マッハ数流れでは以下の関係式；

が成立するため、上記式（１０）をこれらの関係式を用いて式変形すると、下記式（１１）に示すパウエル式が得られる。なお、上記の関係式のｃ_０は静止状態のときの流れ場の音速であり、ρ_０は静止状態のときの流れ場の密度である。

（時間ドメインにおける表面圧力変動と流れ場との関係式の導出）
グリーン（Green）関数Ｇは、下記式（１２）の関係式を満たす。

ここで、観測点ｘ及び評価点ｚは上述したとおりであり（図６参照）、τは過去の時刻である。なお、δ（ｘ）はディラックデルタ関数である。加えて、上記式（１２）中のグリーン関数Ｇは、自由空間においては下記式（１３）のように表される。

物体境界条件を考慮すると、上記式（１１）に示すパウエル式は、下記式（１４）のように表される。

Ｈ（ｚ）はヘビサイド関数（Heaviside Unit Function）であり、下記式（１５）によって定義される。

ここで、図６に示すように、Ｖはフロードメイン３８を表すかたまりであり、Ｓ＋は剛体３６の表面Ｓを取り囲む仮想表面である。上記式（１５）中の「z in V」とは、評価点ｚがフロードメイン３８の内部に位置している場合を示し、「z inside S+」とは、評価点ｚが仮想表面Ｓ＋の内部に位置している場合を示す。

このヘビサイド関数Ｈ（ｚ）は、下記式（１６）の積分方程式を満たす。

ここで、式中の（・）は任意の値との内積を表す。加えて、図６に示すように、剛体３６は、その表面Ｓのうち観測点ｘとは異なる位置に評価点ｙを有する。式中のｎは、当該評価点ｙの法線ベクトルｎ（ｙ）を表し、ｄＳは、評価点ｙ近傍の表面面積ｄＳ（ｙ）を表す。

上記式（１４）の左辺の１項目、２項目及び右辺は、それぞれ下記式(14-1)、(14-2)及び(14-3)のように変形できる。

上記式(14-1)～(14-3)を上記式（１４）に代入すると、下記式（１７）が得られる。

上記式（１７）の右辺の２項目（－符号は省略）は、下記式(17-1)のように式変形できる。

上記式(17-1)を上記式（１７）に代入すると、下記式（１８）が得られる。

上記式（１３）を用いると、上記式（１８）の積分表示は下記式（１９）のようになる。

上記式（１９）の右辺の１項目を上記式（１８）を用いて展開すると、下記式(19-1)が得られる。なお、以下では、表記の都合上、上記式（１９）の右辺の各項の式変形を記載する際は、時間積分を省略して記載する。

上記式（１９）の右辺の２項目及び３項目を上記式（１８）を用いて展開すると、下記式(19-2)が得られる。

上記式（１９）の右辺の４項目は、下記式(19-3)のように式変形できる。

上記式(19-1)～(19-3)を上記式（１９）に代入すると、下記式（２０）が得られる。

仮想表面Ｓ＋を剛体３６の表面Ｓに収束させると、上記式（２０）は下記式（２１）のように書き直すことができる。

ここで、高レイノルズ数流れでは、上記式（２１）の右辺の３項目は４項目に対して無視できるほど小さい。加えて、剛体３６が静止している場合は流速がないため、当該右辺の２項目はゼロに等しい。従って、上記式（２１）は、下記式（２２）のように簡略化できる。

静止物体表面上では流速がないため、エンタルピーＢは下記式のように表される。

上記式の右辺ｐ／ρをｐ／ρ＝ｐ_＊と定義してエンタルピーＢを上記式（２２）に代入すると、下記式（２３）の関係式が得られる。この関係式が、「時間ドメインにおける表面圧力変動と流れ場との関係式」である。

（ＡＰＤＳ計算式の導出）
上記式（２３）の関係式は、表面圧力変動ｐに基づく値（関数）ｐ_＊（＝ｐ／ρ）をその両辺に含むため、観測点ｘにおける表面圧力変動ｐを流速ｕ及び渦度ωから直接的に計算することができない。このため、以下では、式変形により右辺からｐ_＊を含む項を消去し、「観測点ｘにおける表面圧力変動ｐ」と「流れ場の流速ｕ及び渦度ω」との関係式を導出する。

図７に示すように、剛体３６の表面Ｓ上の観測点ｘの近傍領域を表面領域Ｓｘと定義し、表面領域Ｓｘの法線ベクトルをｎ（ｘ）とする。表面領域Ｓｘは十分に大きい平面である。加えて、法線ベクトルｎ（ｘ）上に、観測点ｘからの距離が等しい２個の点ｘ＋、ｘ－を定義する。点ｘ＋は剛体３６の外部に位置しており、点ｘ－は剛体３６の内部に位置している。なお、その他の文字については後述する。点ｘ＋、ｘ－をそれぞれ観測点として上記式（２２）に代入すると、下記式(22-1)、(22-2)が得られる。

式(22-1)と式(22-2)の両辺をそれぞれ足し合わせてｘ－とｘ＋についてｘに近づくときの極限をとると、下記式（２４）が得られる。なお、以下の式変形では、Ｈ（ｘ＋）＝１及びＨ（ｘ－）＝０を利用している（式（１５）参照）。

上記式（２４）の右辺の１項目及び２項目は、それぞれ下記式(24-1)及び(24-2)のように式変形できる。なお、式中のＳ＼Ｓｘは、剛体３６の表面Ｓのうち表面領域Ｓｘを含まない部分を表す（図７参照）。

上記式(24-1)、(24-2)を上記式（２４）に代入すると、下記式（２５）が得られる。

フーリエ変換は下記式によって定義される。

ここで、式中のＡ_θ（ｘ）は角周波数θにおける関数Ａ（ｘ，ｔ）の大きさ（振幅）を表す。Ａ_θ（ｘ）は複素数であり、その複素共役Ａ_θ ^＊（ｘ）は下記（２６）の関係式を満たす。

上記式（２４）にフーリエ変換を作用し、畳み込み積分のフーリエ変換結果を用いると、下記式（２７）が得られる。

ここで、グリーン関数Ｇをフーリエ変換したＧ_θ（ｘ，ｚ）は下記式（２８）のように表され、下記式（２９）の関係式（いわゆるヘルムホルツ(Helmholtz)方程式）を満たす。

ここで、ｋ_０は角周波数θ成分の波数（いわゆる音響波数）であり、ｋ_０＝θ／ｃ_０である。ｒは観測点ｘと評価点ｚとの距離である。

Ｇ_θ（ｘ，ｚ）の勾配は、上記式（２８）を用いて下記式（３０）のように表すことができる。なお、式中のベクトルｒは、観測点ｘから評価点ｚへ向かうベクトルである。

ここで、剛体３６の表面に生じる圧力は流体流動圧力変動と音圧とにより構成されている。流体流動圧力変動は流体流速で移動しそれほど遠方まで伝播しないため、音響コンパクト条件（ｋ_０ｒ＜＜１）を満たす。従って、上記式（３０）のＧ_θ（ｘ，ｚ）の勾配のうち表面流体流動圧力変動に寄与する部分は、下記式(30-1)で表されるように角周波数θに依存しておらず、距離ｒの関数のみ近似できる。

一方、音圧は音速で移動し遠方まで伝播するため、音響コンパクト条件は成立しない（即ち、遠方点では（ｋ_０ｒ＞＞１）となる。）。このため、上記式（３０）のＧ_θ（ｘ，ｚ）の勾配のうち表面音圧変動に寄与する部分は下記式(30-2)のように表される。

上記式（２７）の左辺である「観測点ｘにおける表面圧力変動ｐ_θ（ｘ）に基づく関数ｐ_＊θ（ｘ）」の自己相関関数（上記式（１）参照）は、下記式（３１）のように表される。

上記式（３１）は、Ｇ_θ（ｘ，ｚ）の勾配（式（３０）参照）の音響コンパクト条件を用いて下記式（３２）のように表される。

ここで、着目角周波数帯[θ_ｌ，θ_ｈ]における任意の観測点ｘの物理量の相関距離ｌ_ｘ,ｃｒ、相関面積Ｓ_ｘ,ｃｒ（図７参照）及び相関体積Ｖ_ｘ,ｃｒをそれぞれ下記式(33-1)、(33-2)及び(33-3)によって定義する。なお、式中のθ_ｍは着目角周波数帯の中心角周波数である（θ_ｍ＝２πｆ_ｍ）。

即ち、相関面積Ｓ_ｘ,ｃｒは平均流速の二乗に比例しており、相関体積Ｖ_ｘ,ｃｒは平均流速の三乗に比例している。

本実施形態の着目角周波数帯は比較的に大きいため、相関距離ｌ_ｘ,ｃｒは比較的に小さい。従って、相関面積Ｓ_ｘ,ｃｒ及び相関体積Ｖ_ｘ,ｃｒも比較的に小さい。このため、図７に示すように、相関面積Ｓ_ｘ,ｃｒは表面領域Ｓｘ内に含まれている（Ｓ_ｘ,ｃｒ⊂Ｓｘ）。ここで、上記式（３２）の右辺の１項目の被積分項には剛体３６の表面Ｓ上の評価点ｙが含まれており、表面Ｓのうち表面領域Ｓｘを含まない領域（Ｓ＼Ｓｘ）で積分されている。このため、観測点ｘの物理量は評価点ｙの物理量と無相関となり、この結果、上記式（３２）の右辺の１項目はゼロに等しくなる。

加えて、以下では、任意の関数Ｆ_θとその複素共役関数Ｆ_θ ^＊との積を「関数Ｆ_θのノルム」と称し、下記式（３４）のように定義する。

即ち、自己相関関数（式（１）参照）は、任意の関数Ｆ_θのノルムを着目角周波数帯において積分した値である。

これにより、上記式（３２）は、下記式（３５）のように表すことができる。

上記式（３５）によると、着目角周波数帯における観測点ｘの表面圧力変動ｐ_θ（ｘ）は、フロードメイン３８内の流れ場の渦度ωと流速ｕの外積に依存している。本願発明者はこの点に着目し、上記式（３５）の右辺の空間積分の自己相関関数の値を低減すれば、観測点ｘにおける表面圧力変動ｐ_θ（ｘ）を低減でき、結果として表面音圧変動ｐａ_θ（ｘ）も低減できると考えた。この自己相関関数は、下記式（３６）に示すように式変形できる。

ここで、式中のｚ１、ｚ２は、フロードメイン３８内の任意の評価点であり、ｚ１≠ｚ２である。着目角周波数帯における評価点ｚ１の物理量と評価点ｚ２の物理量との相関は、評価点ｚ１が評価点ｚ２の相関体積Ｖ_{ｚ２,ｃｒ}に含まれていない場合は無視することができる。このため、上記式（３６）は、更に下記式（３７）のように変形できる。

上記式（３７）の右辺の「着目角周波数帯について積分される部分」のうち表面音圧変動に寄与する部分は、上記式(30-2)を用いて下記式（３８）のように近似できる。

上記式（３８）の流速ｕは、下記式(39-1)、(39-2)に示すように、平均成分としての平均流速Ｕ－と、乱れ成分としての乱れ流速ｕ～に分解でき、渦度ωは、平均成分としての平均渦度Ω－と、乱れ成分としての乱れ渦度ω～に分解できる。

加えて、 剛体３６の空間座標系において、ベクトルｒを下記式（４０）によって規定する。

ここで、ｒ１、ｒ２及びｒ３は、観測点ｘを原点としたときのベクトルｒのｅ１成分、ｅ２成分及びｅ３成分である。

これらの式(39-1)、(39-2)及び（４０）を上記式（３８）の右辺に代入すると、当該右辺は、下記式（４１）のように式変形できる。なお、乱れ成分同士の積は無視できるほど小さいため、下記式（４１）では、その記載を省略している。

ここで、２つの任意な関数Ｅ_θと関数Ｆ_θの相関関数を着目角周波数帯において下記式（４２）のように定義する。

上記式（１）及び式（４２）を用いると、渦一様等方性仮説より、下記式(43-1)～(43-6)に示す関係式が得られる。

ここで、ｋは流速ｕ_θと渦度ω_θとの関係を規定した定数である。

表面流体流動圧力変動への寄与度を算出する際は、空間微分と時間微分との間に以下の関係式；

が成立するため、乱れ渦度ω～の自己相関関数は下記式(43-7)のように近似できる。

一方、表面音圧変動への寄与度を算出する際は、音速で移動する音圧の空間微分と時間微分との間に以下の関係式；

が成立するため、乱れ渦度ω～の自己相関関数は下記式(43-8)のように近似できる。

上記式(43-1)～(43-8)を用いて上記式（４１）の右辺を整理すると、上記式（３８）の左辺は、最終的に下記式（４４）のように表される。

上記式（４４）の右辺は、上記式（３７）の右辺の「着目角周波数帯について積分される部分」のうち表面音圧変動に寄与する部分に相当する。このため、式（４４）の右辺を式（３７）の該当部分に代入すると、下記式（４５）が得られる。

上記式（４５）の右辺の被積分項に「流れ場の密度の二乗ρ^２」を乗算し、評価点ｚ２を一般の表記ｚに書き換えると、下記式（４６）に定義されるＡＰＤＳの計算式が得られる。本明細書では、このようにして導出された式（４６）を、「剛体３６の観測点ｘにおける表面音圧変動ｐａ_θ（ｘ）」に対する「フロードメイン３８内の評価点ｚにおける流れ場」の寄与度を表す指標として定義している。なお、ρ^２を乗算するのは、ＡＰＤＳ計算式の導出過程でｐ^＊＝ｐ／ρを用いたことによるものである。

なお、観測点ｘにおける表面音圧変動ｐａ_θ（ｘ）とＡＰＤＳ（ｘ，ｚ）との間には、下記式（４７）に示す近似式が成立する（式（３５）、（４５）及び（４６）参照）。

即ち、「観測点ｘにおける表面音圧変動ｐａ_θ（ｘ）のノルム」の着目角周波数帯における積分値（別言すれば、観測点ｘにおける表面音圧変動ｐａ_θ（ｘ）の自己相関関数）は、「観測点ｘについてのＡＰＤＳをフロードメインで空間積分した値」として近似することができる。このため、ＡＰＤＳは、観測点ｘにおける表面音圧変動ｐａ_θ（ｘ）に対する評価点ｚにおける流れ場の寄与度を表す指標として高い精度で機能する。以上が、ＡＰＤＳ計算式の導出についての説明である。

＜風切音解析方法＞
続いて、走行中の車両モデル２０の風切音解析方法について具体的に説明する。計算装置１０の演算部１４のＣＰＵは、図８にフローチャートにより示したルーチンを実行するようになっている。ＣＰＵは、図８のステップ８００から処理を開始して以下に述べるステップ８０２～ステップ８１２の処理を順に行う。

ステップ８０２：ＣＰＵは、車両モデル２０を走行させる非定常ＣＦＤシミュレーションを開始し、解析範囲（フロードメイン２１内）において、流速の時刻歴データｕ（ｚ，ｔ）及び渦度の時刻歴データω（ｚ，ｔ）を所定時間に亘ってそれぞれ算出する。これらのデータは節点ｚ毎に算出される。
ステップ８０４：ＣＰＵは、ステップ８０２で算出された流速及び渦度の時刻歴データｕ（ｚ，ｔ）及びω（ｚ，ｔ）をそれぞれ平均化処理して、節点ｚ毎に、平均流速Ｕ－（ｚ）及び平均渦度Ω－（ｚ）を算出する。
ステップ８０６：ＣＰＵは、ステップ８０２で算出された流速の時刻歴データｕ（ｚ，ｔ）を高速フーリエ変換処理して、着目角周波数帯における乱れ流速ｕ～_θ（ｚ）の自己相関関数を算出する。この自己相関関数は節点ｚ毎に算出される。ステップ８０２～ステップ８０６までの処理が非定常ＣＦＤ計算処理に相当する。

ステップ８０８：ＣＰＵは、演算部１４のＲＯＭに格納されているＡＰＤＳ計算式（式（４６）参照）を用いて、解析表面上の節点ｘについての着目角周波数帯のＡＰＤＳを節点ｚ毎に算出する。ＣＰＵは、この処理を解析表面上の全ての節点ｘについて行う。即ち、ＡＰＤＳは、全部でｍｎ個算出される（ｍ：節点ｘの個数、ｎ：節点ｚの個数）。この計算式には、ステップ８０４で算出された平均流速Ｕ－（ｚ）及び平均渦度Ω－（ｚ）並びにステップ８０６で算出された乱れ流速ｕ～_θ（ｚ）の自己相関関数等が代入される。
ステップ８１０：ＣＰＵは、ステップ８０８で算出されたｍｎ個のＡＰＤＳから、任意の節点ｘについてのｎ個のＡＰＤＳを抽出し、抽出されたＡＰＤＳをフロードメイン２１で空間積分することにより当該節点ｘにおける着目角周波数帯の表面音圧変動ｐａ_θ（ｘ）の自己相関関数の近似値を算出する。そして、当該近似値に基づいて当該節点ｘにおける予測表面音圧変動ｐａ_ｐθ（ｘ）を算出（予測）する。ＣＰＵは、この処理を解析表面上の全ての節点ｘについて行う。即ち、予測表面音圧変動ｐａ_ｐθ（ｘ）は全部でｍ個算出される。

ステップ８１２：ＣＰＵは、ステップ８１０で算出された着目角周波数帯の予測表面音圧変動ｐａ_ｐθ（ｘ）に基づいて、「着目角周波数帯の中心角周波数θ_ｍにおける予測表面音圧変動ｐａ_ｐθ（ｘ）の分布図を表すデータ」を作成し、当該データの表示指令を出力部１６に送信する。これにより、予測表面音圧変動ｐａ_ｐθ（ｘ）の分布図が表示画面１８に画像表示される（図３Ａ参照）。その後、ＣＰＵは、ステップ８１４に進む。

ステップ８１４では、ＣＰＵは、作業者により観測点ｘが選択されたか否かを判定する。この観測点ｘは、予測表面音圧変動ｐａ_ｐθ（ｘ）の分布図（ステップ８１２参照）に基づいて作業者により選択される。観測点ｘが選択されたと判定した場合、ＣＰＵは、ステップ８１４にて「Ｙｅｓ」と判定し、後述するステップ８１６の処理を行う。これに対し、観測点ｘが選択されていないと判定した場合、ＣＰＵは、ステップ８１４にて「Ｎｏ」と判定し、再度ステップ８１４の判定を行う。ＣＰＵは、観測点ｘが選択されたと判定するまでこの処理を所定の演算時間が経過する毎に繰り返す。

ステップ８１６：ＣＰＵは、ＡＰＤＳの参照値の入力を促すメッセージの表示指令を出力部１６に送信する。これにより、当該メッセージが表示画面１８に表示され、作業者はＡＰＤＳの参照値を入力できる。その後、ＣＰＵは、ステップ８１８に進む。

ステップ８１８では、ＣＰＵは、作業者によりＡＰＤＳの参照値が入力されたか否かを判定する。参照値が入力されたと判定した場合、ＣＰＵは、ステップ８１８にて「Ｙｅｓ」と判定し、以下のステップ８２０の処理を行う。

ステップ８２０：ＣＰＵは、ステップ８０８で算出されたｍｎ個のＡＰＤＳの中から、ステップ８１４で選択された観測点ｘについてのｎ個のＡＰＤＳを抽出する。そして、抽出されたＡＰＤＳの中から、ステップ８１８で入力された参照値を有するＡＰＤＳの等値面を表すデータを作成し、当該データの表示指令を出力部１６に送信する。これにより、作業者が入力した参照値を有するＡＰＤＳの等値面が表示画面１８に画像表示される（図４Ａ～図４Ｃ参照）。その後、ＣＰＵは、ステップ８２２に進む。

ステップ８２２では、ＣＰＵは、作業者により着目評価点ｚが選択されたか否かを判定する。この着目評価点ｚは、ＡＰＤＳの等値面（ステップ８２０参照）に基づいて作業者により選択される。着目評価点ｚが選択されたと判定した場合、ＣＰＵは、ステップ８２２にて「Ｙｅｓ」と判定し、後述するステップ８２４（図９）の処理を行う。

これに対し、着目評価点ｚが選択されていないと判定した場合（具体的には、作業者が現時点におけるＡＰＤＳ等値面からは着目評価点ｚを選択できない場合、又は、作業者が着目評価点ｚを選択中の場合）、ＣＰＵは、ステップ８２２にて「Ｎｏ」と判定し、再度ステップ８１８の判定を行う。

そして、ステップ８１８にて参照値が入力されたと判定した場合（即ち、作業者によりＡＰＤＳの参照値が再入力された場合）、ＣＰＵは、ステップ８１８にて「Ｙｅｓ」と判定してステップ８２０の処理を行い、ステップ８２２にて着目評価点ｚが選択されたか否かを再度判定する。一方、ステップ８１８にて参照値が入力されていないと判定した場合、ＣＰＵは、ステップ８１８にて「Ｎｏ」と判定してステップ８２２に進み、ステップ８２２にて着目評価点ｚが選択されたか否かを再度判定する。ＣＰＵは、着目評価点ｚが選択されたと判定するまでステップ８１８～ステップ８２２の処理を繰り返す。

これに対し、ステップ８１６から直接ステップ８１８に進んだ場合において参照値が入力されていないと判定した場合（即ち、初回の参照値が入力されていない場合）、ＣＰＵは、ステップ８１８にて「Ｎｏ」と判定してステップ８２２に進む。上述したように、着目評価点ｚはＡＰＤＳの等値面に基づいて選択される（即ち、ステップ８２０を経て選択される）ため、初回の参照値が入力されていない場合（即ち、ステップ８２０を経ていない場合）、ＣＰＵはステップ８２２にて「Ｎｏ」と判定する。

ステップ８２４（図９）：ＣＰＵは、ステップ８１８で入力された参照値を有するＡＰＤＳ（ｘ，ｚ）の中から、ステップ８２２で選択された着目評価点ｚにおけるＡＰＤＳ（ｘ，ｚ）を抽出し、抽出されたＡＰＤＳ（ｘ，ｚ）を構成するパラメータの数値データの表示指令を出力部１６に送信する。これらのパラメータは、具体的には、着目評価点ｚにおける平均流速Ｕ－（ｚ）、乱れ流速ｕ～_θ（ｚ）、及び、平均渦度Ω－（ｚ）である（ステップ８０４、８０６参照）。これにより、ＡＰＤＳ（ｘ，ｚ）を構成するパラメータの数値データが表示画面１８に表示される。その後、ＣＰＵは、ステップ８２６に進む。

ステップ８２６では、ＣＰＵは、原因パラメータが特定されたか否かを判定する。この原因パラメータは、ＡＰＤＳ（ｘ，ｚ）を構成するパラメータの数値データ（ステップ８２４参照）に基づいて、作業者により特定される。原因パラメータが特定されたと判定した場合、ＣＰＵは、ステップ８２６にて「Ｙｅｓ」と判定し、後述するステップ８２８の処理を行う。これに対し、原因パラメータが特定されていないと判定した場合、ＣＰＵは、ステップ８２６にて「Ｎｏ」と判定し、再度ステップ８２６の判定を行う。ＣＰＵは、原因パラメータが特定されたと判定するまでこの処理を所定の演算時間が経過する毎に繰り返す。

ステップ８２８：ＣＰＵは、ステップ８２２で選択された着目評価点ｚを通る所定の平面（典型的には、ｅ１ｅ２平面に平行な平面。即ち、着目評価点通過断面。）上における原因パラメータの分布図を表すデータを作成する。ＣＰＵは、当該データの表示指令を出力部１６に送信する。これにより、原因パラメータの分布図が表示画面１８に表示される（図５Ｂ参照）。その後、ＣＰＵは、ステップ８３０に進み、本ルーチンを終了する。

＜本実施装置による表面音圧変動（風切音）の改善効果＞
次に、本実施装置を用いて車両モデル２０の解析範囲の形状を変更したときの着目周波数帯（本例では２ｋＨｚ）の表面音圧変動ｐａ_θ（ｘ）の改善効果について図１０Ａ～図１３を参照して具体的に説明する。図５Ａ及び図５Ｂの例では、領域３４における空気抵抗を増大できるような右側サイドミラー２３の形状が検討される。図１０Ａは、車両モデル２０の正面図の－ｅ２方向側の一部を示し、図１０Ｂは右側サイドミラー２３の平面図の概形を示す。図１０Ａ及び図１０Ｂに示すように、右側サイドミラー２３の車体側（ｅ２方向側）の面２３ａは、車両モデル２０の前進方向（ｅ１方向）と略平行となっている。本願発明者は、これにより面２３ａ近傍の空気抵抗が小さくなっていると考え、面２３ａ近傍の空気抵抗を増大できる形状を検討した。図１１Ａは、この検討に基づいて形状が変更された車両モデル１２０の正面図の－ｅ２方向側の一部を示し、図１１Ｂは車両モデル１２０の右側サイドミラー１２３の平面図の概形を示す。図１１Ａ及び図１１Ｂに示すように、右側サイドミラー１２３の車体側（ｅ２方向側）の面１２３ａは、車両モデル１２０の前進方向（ｅ１方向）に対して交差する向きとなっており、これにより、面１２３ａ近傍の空気抵抗を増大させている。

図１２Ａは、形状変更前（車両モデル２０）の着目評価点通過断面（図１０Ａの破線Ｌ参照）における右側サイドミラー２３の周囲の平均流速Ｕ－（ｚ）の分布図を示す。図１２Ｂは、形状変更後（車両モデル１２０）の着目評価点通過断面（図１１Ａの破線Ｌ参照）における右側サイドミラー１２３の周囲の平均流速Ｕ－（ｚ）の分布図を示す。図１２Ａ及び図１２Ｂに示すように、右側サイドミラーの形状を変更することにより右側サイドミラー１２３の車体側近傍の領域の平均流速Ｕ－（ｚ）は領域３４と比べて大幅に低減している。

図１３は、形状変更後の解析表面（右側フロントサイドガラス１２２）における予測表面音圧変動ｐａ_ｐθ（ｘ）の分布図を示す。図３Ａ（形状変更前の解析表面における予測表面音圧変動ｐａ_ｐθ（ｘ）の分布図）と比較すると、図１３では、観測点ｘにおける予測表面音圧変動ｐａ_ｐθ（ｘ）が大幅に低減している。以上より、本実施装置によれば、ＡＰＤＳを利用して車両モデルの表面音圧変動に大きく寄与している流れ場の位置を特定することにより、車両モデルの高周波領域における表面音圧変動（風切音）を効率的に且つ確実に低減できることが示された。

＜指標としてのＡＰＤＳの信頼性及び汎用性＞
上述したように、「節点ｘについての着目周波数帯のＡＰＤＳをフロードメインで空間積分した値」は、「当該節点ｘにおける着目周波数帯の表面音圧変動ｐａ_θ（ｘ）の自己相関関数」の近似値である（式（４７）参照）。本実施装置は、これを利用してＡＰＤＳに基づいて節点ｘにおける予測表面音圧変動ｐａ_ｐθ（ｘ）を算出している。以下では、予測表面音圧変動ｐａ_ｐθの近似精度を検証してＡＰＤＳの指標としての信頼性について検討する。具体的には、予測表面音圧変動ｐａ_ｐθを、非定常ＣＦＤ計算用のソフトウエアによって算出される表面音圧変動ｐａ_θと比較することにより、予測表面音圧変動ｐａ_ｐθの近似精度を検証する。本実施形態では上記ソフトウエアとしてＰｏｗｅｒＦＬＯＷが使用されるが、ソフトウエアの種類はこれに限られない。なお、本実施装置が表面音圧変動を非定常ＣＦＤ計算により算出するのではなくＡＰＤＳに基づいて算出（予測）する（ステップ８１０参照）のは、非定常ＣＦＤ計算による算出方法では多大な時間がかかるからである。例えば、車両モデル２０の右側フロントサイドガラス２２の表面音圧変動を節点ｘ毎に算出する場合、非定常ＣＦＤ計算では３ヶ月かかるのに対し、ＡＰＤＳに基づく計算では２日で済む。

更に、以下では、車両モデルだけではなく、車両モデルとは異なる形状を有するフォアステップモデルについても上記近似精度を検証することにより、ＡＰＤＳの指標としての汎用性についても検討する。

図１４Ａは、上記の説明で用いられた車両モデル２０である。図１４Ｂは、フォアステップモデル４０の模式図である。フォアステップモデル４０は、基準面４２と、基準面４２上に配置された直方体形状を有するステップ４４と、基準面４２上にステップ４４から離間して配置されたノズル４６と、を有する。ステップ４４の上面４４ａは基準面４２と平行である。ステップ４４の側面４４ｂ（ノズル４６側の面）と基準面４２との間には段差が形成されている。ノズル４６は、ステップ４４の側面４４ｂと平行に対向している送風口４６ａを有する。送風口４６ａは、ステップ４４に向かって風を送る。送風方向と直交する方向を「幅方向」と規定すると、幅方向におけるステップ４４の長さａは、幅方向におけるノズル４６の長さｂよりも長い。ステップ４４の高さはノズル４６の高さよりも低い。ノズル４６の幅方向における中心は、ステップ４４の幅方向における中心と、幅方向において一致している。本検証では、車両モデル２０の右側フロントサイドガラス２２とフォアステップモデル４０のステップ４４の上面４４ａが解析対象とされる。

図１５Ａ及び図１５Ｂは、それぞれ、車両モデル２０の右側フロントサイドガラス２２における予測表面音圧変動ｐａ_ｐθ及び表面音圧変動ｐａ_θの分布図である。着目周波数帯は２ｋＨｚに設定されている。図１６Ａ及び図１６Ｂは、それぞれ、フォアステップモデル４０のステップ４４の上面４４ａにおける予測表面音圧変動ｐａ_ｐθ及び表面音圧変動ｐａ_θの分布図である。着目周波数帯は２ｋＨｚに設定されている。予測表面音圧変動ｐａ_ｐθ及び表面音圧変動ｐａ_θが大きいほど濃色で表示されている。図１５Ａと図１５Ｂを比較すると、予測表面音圧変動ｐａ_ｐθの大きさ及び位置は、特に値が大きい領域において、表面音圧変動ｐａ_θの大きさ及び位置と精度良く一致している。同様に、図１６Ａと図１６Ｂを比較すると、予測表面音圧変動ｐａ_ｐθの大きさ及び位置は、特に値が大きい領域において、表面音圧変動ｐａ_θの大きさ及び位置と精度良く一致している。以上より、ＡＰＤＳに基づいて算出される予測表面音圧変動ｐａ_ｐθは、ソフトウエアのＣＦＤ計算により算出される表面音圧変動ｐａ_θに対して高い近似精度を有することが示された。これは、ＡＰＤＳが指標として高い信頼性を有することを示している。

図１７は、「表面音圧変動ｐａ_θに対する予測表面音圧変動ｐａ_ｐθの誤差」を車両モデル２０及びフォアステップモデル４０についてそれぞれ示したグラフである。車両モデル２０についての誤差は、図１５Ａ及び図１６Ｂに示す２つの分布図の対応する各節点ｘにおける「表面音圧変動ｐａ_θに対する予測表面音圧変動ｐａ_ｐθの差」の平均を音圧変動レベル（dB）に変換した値である。フォアステップモデル４０についての誤差は、図１６Ａ及び図１６Ｂに示す２つの分布図の対応する各節点ｘにおける「表面音圧変動ｐａ_θに対する予測表面音圧変動ｐａ_ｐθの差」の平均を音圧変動レベル（dB）に変換した値である。車両モデル２０の誤差とフォアステップモデル４０の誤差を比較すると、何れも約３ｄＢでありほぼ一致している。以上より、表面音圧変動ｐａ_θに対する予測表面音圧変動ｐａ_ｐθの誤差は、解析対象モデルの形状に関係なくほぼ一定であることが示された。これは、ＡＰＤＳが指標として高い汎用性を有することを示している。

本実施装置の作用効果について説明する。本実施装置では、車両モデルを走行させる非定常ＣＦＤシミュレーションにより算出された物理量（平均流速Ｕ－、平均渦度Ω－及び乱れ流速ｕ～_θ）に基づいて着目角周波数帯のＡＰＤＳを算出する。このＡＰＤＳ（ｘ，ｚ）は、「車両モデル２０の解析表面（右側フロントサイドガラス２２）上の観測点ｘにおける着目角周波数帯の表面音圧変動ｐａ_θ（ｘ）」に対する「解析範囲（フロードメイン２１）内の評価点ｚにおける流れ場」の寄与度を示す指標である。このため、ＡＰＤＳを指標として用いることにより、着目角周波数帯の表面音圧変動ｐａ_θ（ｘ）に対する流れ場の寄与度を評価点ｚ毎に算出できる。従って、フロードメイン２１内のどの位置における流れ場が観測点ｘにおける車両モデル２０のガラス製の部分（本実施形態では右側フロントサイドガラス２２）における高周波領域の表面音圧変動ｐａ_θ（ｘ）に大きく寄与しているかを適切に特定することが可能となる。

加えて、本実施装置では、観測点ｘについてのＡＰＤＳをフロードメイン２１で空間積分した値は、当該観測点ｘにおける表面音圧変動ｐａ_θ（ｘ）の自己相関関数の近似値である。このため、ＡＰＤＳの挙動は、表面音圧変動ｐａ_θ（ｘ）の挙動と精度良く一致する。即ち、ＡＰＤＳは、表面音圧変動ｐａ_θ（ｘ）に対する流れ場の寄与度を表す指標として高い信頼性を有する（図１４Ａ～図１７参照）。このため、車両モデルのガラス製の部分における高周波領域の表面音圧変動ｐａ_θ（ｘ）に大きく寄与している流れ場の位置を高い精度で特定することができる。

具体的には、従来は、上記関係を定量的に評価することができなかったため、技術者によって流れ場の解釈が異なることがあった。加えて、解析は、流れ場の断面における物理量の分布に基づいて行われるため、三次元空間における流れ場の情報が欠落し、表面音圧変動が大きい原因を見逃してしまう可能性があった。これに対し、本実施装置によれば、ＡＰＤＳを用いることにより表面音圧変動ｐａ_θ（ｘ）と流れ場との関係を定量的に評価できるため、技術者間で流れ場の解釈が異なる可能性を大幅に低減でき、風切音の解析精度を向上することができる。加えて、本実施装置では、表面音圧変動への寄与度が高い流れ場の位置を特定してから当該位置を通る断面における物理量の分布を評価するため、流れ場の情報が欠落する可能性を大幅に低減でき、結果として、表面音圧変動が大きい原因を確実に発見することができる。

更に、本実施装置では、表示画面１８にＡＰＤＳの等値面が表示されるため、作業者は、参照値を適宜設定することにより、所望のＡＰＤＳ値を有する等値面を視認することができる。この結果、所望のＡＰＤＳ値に対応している評価点ｚを効率的に選択することができ、車両モデル２０のガラス製の部分における高周波領域の表面音圧変動ｐａ_θ（ｘ）に大きく寄与している流れ場の位置を効率的に特定することができる。

更に、本実施装置では、観測点ｘは、ＡＰＤＳに基づいて算出（予測）された予測表面音圧変動ｐａ_ｐθの大きさに基づいて選択される。予測表面音圧変動ｐａ_ｐθは風切音と強い相関を有している。このため、この構成によれば、風切音をより適切に解析することができる。

（変形例１）
次に、変形例１に係る風切音解析装置（以下、「第１変形装置」とも称する。）について説明する。第１変形装置は、原因パラメータの特定が作業者ではなく装置により行われる点で上記実施装置と異なっている。具体的には、第１変形装置は、図９のステップ８２４及び８２６の処理に代えて、以下の処理を実行する。即ち、第１変形装置の計算装置の演算部は、ステップ８２２で選択された着目評価点ｚにおけるＡＰＤＳを構成する３つのパラメータ（平均流速Ｕ－、平均渦度Ω－及び乱れ流速ｕ～_θ）の中からＡＰＤＳ値への寄与が相対的に大きいパラメータを原因パラメータとして特定する。具体的には、演算部は、着目評価点ｚにおける３つのパラメータの数値を、「着目評価点ｚにおけるＡＰＤＳ値よりも低いＡＰＤＳ値に対応している評価点ｚ」における３つのパラメータの数値と比較することにより原因パラメータを特定する。この構成によれば、原因パラメータを作業者が特定する必要がなくなるため、ＡＰＤＳが高い原因となっているパラメータを容易に把握でき、車両モデルの形状をより効率的に検討、変更することができる。

（変形例２）
続いて、変形例２に係る風切音解析装置（以下、「第２変形装置」とも称する。）について説明する。第２変形装置は、着目評価点ｚの選択（抽出）が作業者ではなく装置により行われる点で上記実施装置と異なっている。具体的には、第２変形装置は、図８のステップ８１６～ステップ８２２の処理に代えて、以下の処理を実行する。即ち、第２変形装置の計算装置の演算部は、ステップ８１４にて「Ｙｅｓ」と判定した場合、任意の値を有するＡＰＤＳの等値面データを作成し、当該データにより示される等値面が数個になるまでＡＰＤＳの値を変更しながら等値面データを作成する（厳密には、数個になってからも各等値面の表面積ができるだけ小さくなるまでＡＰＤＳの値を増加する。）。そして、観測点ｘまでの距離が最も長い等値面上の任意の点を着目評価点ｚとして抽出する。この構成によれば、着目評価点ｚを作業者が選択する必要がなくなるため、車両モデルのガラス製の部分における高周波領域の表面音圧変動に大きく寄与している流れ場の位置をより効率的に特定することができる。なお、第２変形装置では最大（ほぼ最大）のＡＰＤＳ値に対応している評価点ｚが着目評価点ｚとして抽出されるが、この構成に限られず、予め設定された条件を満たす着目評価点ｚが抽出される構成であっても良い。

以上、実施形態及び変形例の風切音解析装置及び風切音解析方法について説明したが、本発明は上記実施形態及び変形例に限定されるものではなく、本発明の目的を逸脱しない限りにおいて種々の変更が可能である。

例えば、上記実施形態及び変形例では解析対象となる構造物モデルとして車両モデルが用いられたが、この構成に限られず、例えば、航空機及び船舶等の構造物モデル等が用いられてもよい。

加えて、上記実施形態及び変形例では、表面音圧変動が比較的に大きい観測点が解析対象とされたが、この構成に限られず、例えば、所定範囲における表面音圧変動の平均値を有する節点ｘを観測点として解析してもよい。

更に、上記実施形態及び変形例では、等値面上の節点ｚが着目評価点ｚとして選択可能であったが、この構成に限られず、例えば、等値面内部の節点ｚが着目評価点ｚとして選択可能であってもよい。

更に、上記実施形態及び変形例では、観測点ｘは作業者により選択されたが、この構成に限られず、例えば、表面音圧変動の値に基づいて風切音解析装置により選択されてもよい。

更に、上記実施形態及び変形例では、観測点ｘは表面音圧変動の分布図に基づいて選択されたが、この構成に限られず、例えば、表面音圧変動以外の物理量に基づいて選択されてもよいし、任意の節点ｘが作業者により観測点ｘとして選択されてもよい。

１０：計算装置、１２：入力部、１４：演算部、１６：出力部、１８：表示画面、２０：車両モデル、２１：フロードメイン、２２：右側フロントサイドガラス、２３：右側サイドミラー、２４，２６：領域、２８，３０，３２（３２ａ，３２ｂ，３２ｃ）：等値面、３４：領域、３６：剛体モデル、３８：フロードメイン、４０：フォアステップモデル、４２：基準面、４４：ステップ、４６：ノズル

Claims (8)

1. 移動中の構造物のガラス製の部分であるガラス部の表面に発生する風切音を解析する風切音解析装置であって、
   構造物をモデル化した構造物モデルを移動させる非定常ＣＦＤシミュレーションを実行して、前記構造物モデルの周囲の流れ場のうち所定領域内の流れ場の所定時間における平均流速及び平均渦度を、当該所定領域内の節点である空間節点毎に算出するとともに、風切音の解析対象となる角周波数帯である着目角周波数帯において、前記所定領域内における乱れ流速の振幅に基づく値を、前記空間節点毎に算出する非定常ＣＦＤ計算手段と、
   前記構造物モデルのガラス部表面の風切音の解析対象点における着目角周波数帯の音圧変動の振幅である表面音圧変動に対する前記所定領域内の空間節点における流れ場の寄与度を示す指標である音圧ソース密度を、前記非定常ＣＦＤ計算手段によって算出された前記平均流速と、前記平均渦度と、前記乱れ流速の振幅に基づく値と、に基づいて算出する音圧ソース密度算出手段と、
   を備える風切音解析装置。
2. 請求項１に記載の風切音解析装置であって、
   前記構造物モデルのガラス部表面の前記解析対象点についての着目角周波数帯の前記音圧ソース密度を前記所定領域で空間積分した値は、前記解析対象点における前記表面音圧変動の関数とその複素共役関数との積を着目角周波数帯で積分した値の近似値である、
   風切音解析装置。
3. 請求項１又は請求項２に記載の風切音解析装置であって、
   更に、前記音圧ソース密度を構成する複数のパラメータのうち、当該音圧ソース密度への寄与が相対的に大きいパラメータである原因パラメータを特定する原因パラメータ特定手段を備える、
   風切音解析装置。
4. 請求項３に記載の風切音解析装置であって、
   前記複数のパラメータは、前記平均流速、前記平均渦度及び前記乱れ流速である、
   風切音解析装置。
5. 請求項１～４の何れか一項に記載の風切音解析装置であって、
   更に、前記音圧ソース密度が算出された複数の空間節点の中から、外部から入力される値を有する音圧ソース密度に対応している複数の空間節点を抽出し、当該抽出された複数の空間節点を画像処理して等値面を作成し、当該等値面を可視化する画像処理手段を備える、
   風切音解析装置。
6. 請求項１～４の何れか一項に記載の風切音解析装置であって、
   更に、前記音圧ソース密度算出手段によって算出された複数の音圧ソース密度の最大値に対応している空間節点を抽出する空間節点抽出手段を備える、
   風切音解析装置。
7. 請求項１～６の何れか一項に記載の風切音解析装置であって、
   前記音圧ソース密度の計算式は、下記数１により表される式によって定義される、
   風切音解析装置。
   

   
8. 計算装置を用いて移動中の構造物のガラス製の部分であるガラス部の表面に発生する風切音を解析する風切音解析方法であって、
   前記計算装置が備える非定常ＣＦＤ計算手段が、構造物をモデル化した構造物モデルを移動させる非定常ＣＦＤシミュレーションを実行して、前記構造物モデルの周囲の流れ場のうち所定領域内の流れ場の流速の時刻歴データ及び渦度の時刻歴データを、当該所定領域内の節点である空間節点毎に所定時間に亘ってそれぞれ算出するステップと、
   前記非定常ＣＦＤ計算手段が、前記算出された流速の時刻歴データ及び渦度の時刻歴データをそれぞれ平均化処理して、平均流速及び平均渦度を前記空間節点毎に算出するステップと、
   前記非定常ＣＦＤ計算手段が、前記算出された流速の時刻歴データを高速フーリエ変換処理して、乱れ流速の振幅の関数とその複素共役関数との積を風切音の解析対象となる角周波数帯である着目角周波数帯で積分した値である乱れ流速の振幅の自己相関関数を前記空間節点毎に算出するステップと、
   前記計算装置が備える音圧ソース密度算出手段が、前記構造物モデルのガラス部表面の風切音の解析対象点における着目角周波数帯の音圧変動の振幅である表面音圧変動に対する前記所定領域内の空間節点における流れ場の寄与度を示す指標である音圧ソース密度を、予めＲＯＭに格納されている音圧ソース密度の計算式であって、平均流速と、平均渦度と、乱れ流速の振幅の自己相関関数と、をパラメータとして含む計算式、に、前記非定常ＣＦＤ計算手段によって前記空間節点毎に算出された前記平均流速と、前記平均渦度と、前記乱れ流速の振幅の自己相関関数と、のうち、対応する空間節点における平均流速と、平均渦度と、乱れ流速の振幅の自己相関関数と、を代入することにより算出するステップと、
   を備える風切音解析方法。

Images