(57)【要約】      （修正有）【課題】構造物の表面圧力変動に大きく寄与している構造物周囲の流れ場の位置を適切に特定する技術を提供する。【解決手段】風切音解析装置は、平均流速及び平均渦度を、当該所定領域内の節点である空間節点毎に算出するとともに、着目角周波数帯において、所定領域内における乱れ流速の振幅に基づく値を、空間節点毎に算出する非定常ＣＦＤ計算手段と、構造物モデル表面の風切音の解析対象点における圧力変動の振幅である表面圧力変動に対する、所定領域内の空間節点における流れ場の寄与度を示す指標である圧力ソース密度を、非定常ＣＦＤ計算手段によって算出された平均流速と、平均渦度と、乱れ流速の振幅に基づく値と、に基づいて算出する圧力ソース密度算出手段とを備える。【選択図】図５Ｂ
2019191110.tif 000002
【特許請求の範囲】
【請求項１】
移動中の構造物の表面に発生する風切音を解析する風切音解析装置であって、
構造物をモデル化した構造物モデルを移動させる非定常ＣＦＤシミュレーションを実行して、前記構造物モデルの周囲の流れ場のうち所定領域内の流れ場の所定時間における平均流速及び平均渦度を、当該所定領域内の節点である空間節点毎に算出するとともに、風切音の解析対象となる角周波数帯である着目角周波数帯において、前記所定領域内における乱れ流速の振幅に基づく値を、前記空間節点毎に算出する非定常ＣＦＤ計算手段と、
前記構造物モデル表面の風切音の解析対象点における圧力変動の振幅である表面圧力変動に対する、前記所定領域内の空間節点における流れ場の寄与度を示す指標である圧力ソース密度を、前記非定常ＣＦＤ計算手段によって算出された前記平均流速と、前記平均渦度と、前記乱れ流速の振幅に基づく値と、に基づいて算出する圧力ソース密度算出手段と、
を備える風切音解析装置。
【請求項２】
請求項１に記載の風切音解析装置であって、
前記構造物モデル表面の前記解析対象点についての前記圧力ソース密度を前記所定領域で空間積分した値は、前記解析対象点における前記表面圧力変動の関数とその複素共役関数との積を着目角周波数帯で積分した値の近似値である、
風切音解析装置。
【請求項３】
請求項１又は請求項２に記載の風切音解析装置であって、
更に、前記圧力ソース密度を構成する複数のパラメータのうち、当該圧力ソース密度への寄与が相対的に大きいパラメータである原因パラメータを特定する原因パラメータ特定手段を備える、
風切音解析装置。
【請求項４】
請求項３に記載の風切音解析装置であって、
前記複数のパラメータは、前記平均流速、前記平均渦度及び前記乱れ流速である、
風切音解析装置。
【請求項５】
請求項１〜４の何れか一項に記載の風切音解析装置であって、
更に、前記圧力ソース密度が算出された複数の空間節点の中から、外部から入力される値を有する圧力ソース密度に対応している複数の空間節点を抽出し、当該抽出された複数の空間節点を画像処理して等値面を作成し、当該等値面を可視化する画像処理手段を備える、
風切音解析装置。
【請求項６】
請求項１〜４の何れか一項に記載の風切音解析装置であって、
更に、前記圧力ソース密度算出手段によって算出された複数の圧力ソース密度の最大値に対応している空間節点を抽出する空間節点抽出手段を備える、
風切音解析装置。
【請求項７】
請求項１〜６の何れか一項に記載の風切音解析装置であって、
前記非定常ＣＦＤ計算手段は、更に、前記着目角周波数帯において、前記構造物モデルの表面における前記表面圧力変動を、当該表面上の節点である表面節点毎に算出する、
風切音解析装置。
【請求項８】
請求項１〜７の何れか一項に記載の風切音解析装置であって、
前記圧力ソース密度の計算式は、下記式によって定義される、
風切音解析装置。
【数１】
2019191110.tif 000066
【請求項９】
移動中の構造物の表面に発生する風切音を解析する風切音解析方法であって、
構造物をモデル化した構造物モデルを移動させる非定常ＣＦＤシミュレーションを実行して、前記構造物モデルの周囲の流れ場のうち所定領域内の流れ場の所定時間における平均流速及び平均渦度を、当該所定領域内の節点である空間節点毎に算出するとともに、風切音の解析対象となる角周波数帯である着目角周波数帯において、前記所定領域内における乱れ流速の振幅に基づく値を、前記空間節点毎に算出する非定常ＣＦＤ計算工程と、
前記構造物モデル表面の風切音の解析対象点における圧力変動の振幅である表面圧力変動に対する、前記所定領域内の空間節点における流れ場の寄与度を示す指標である圧力ソース密度を、前記非定常ＣＦＤ計算工程によって算出された前記平均流速と、前記平均渦度と、前記乱れ流速の振幅に基づく値と、に基づいて算出する圧力ソース密度算出工程と、
を備える風切音解析方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、風切音解析装置及び風切音解析方法に関する。
【背景技術】
【０００２】
車両走行時の風切音を低減するための技術が研究、開発されている。例えば、特許文献１には、車両モデルを走行させるシミュレーションにより当該車両モデル表面の任意の位置における「圧力変動の振幅」と「平均流速」とを取得し、これらの値に基づいて当該任意の位置における風切音の音源の強度を所望の周波数帯毎に算出する技術が開示されている。なお、以下では、車両モデル表面における圧力変動の振幅を、単に「表面圧力変動」とも称する。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特開２０１７−０６２７２７号公報
【発明の概要】
【０００４】
ここで、風切音の音源の強度は、表面圧力変動と平均表面流速のうち、特に表面圧力変動と強く相関している。表面圧力変動は、車両周囲の流れ場（典型的には、流速及び渦度）の影響を受けて変化する。このため、車両表面のある特定の位置の風切音を低減するためには、車両周囲の流れ場のうちどの部分（位置）の流れ場が当該特定の位置の表面圧力変動に強く影響を与えているのかを特定し、当該特定の位置の表面圧力変動が低減するように当該部分の流れ場を変更する（即ち、車両の部品形状を変更する）必要がある。しかしながら、「表面圧力変動に大きく寄与している流れ場の部分を特定する方法」には統一された手法があるわけではなく、実際には技術者が自身のノウハウ及び知見に基づいて当該流れ場の部分を特定している。このため、技術者によって検討結果にばらつきが生じ、試行錯誤が多く発生するという問題があった。
【０００５】
本発明は、上記の問題を解決するためになされたものであり、その目的の一つは、車両を始めとする構造物の表面圧力変動に大きく寄与している構造物周囲の流れ場の位置を適切に特定することが可能な技術を提供することにある。
【０００６】
本発明の風切音解析装置は、移動中の構造物の表面に発生する風切音を解析する。
この風切音解析装置は、
構造物をモデル化した構造物モデル（２０）を移動させる非定常ＣＦＤシミュレーションを実行して、前記構造物モデル（２０）の周囲の流れ場のうち所定領域（２１）内の流れ場の所定時間における平均流速（Ｕ−）及び平均渦度（Ω−）を、当該所定領域（２１）内の節点である空間節点（ｚ）毎に算出するとともに、風切音の解析対象となる角周波数帯である着目角周波数帯（[θｌ，θｈ]）において、前記所定領域（２１）内における乱れ流速の振幅（ｕ〜θ）に基づく値を、前記空間節点（ｚ）毎に算出する非定常ＣＦＤ計算手段（ステップ８０２、８０６、８０８）と、
前記構造物モデル（２０）表面の風切音の解析対象点（ｘ）における圧力変動の振幅である表面圧力変動（ｐθ）に対する、前記所定領域（２１）内の空間節点（ｚ）における流れ場の寄与度を示す指標である圧力ソース密度（ＰＤＳ）を、前記非定常ＣＦＤ計算手段（ステップ８０２、８０６、８０８）によって算出された前記平均流速（Ｕ−）と、前記平均渦度（Ω−）と、前記乱れ流速の振幅（ｕ〜θ）に基づく値と、に基づいて算出する圧力ソース密度算出手段（ステップ８１４）と、
を備える。
【０００７】
本発明の風切音解析装置（以下、「本発明装置」とも称する。）は、構造物モデルを移動させる非定常ＣＦＤシミュレーションにより算出された物理量（平均流速、平均渦度及び乱れ流速）に基づいて圧力ソース密度を算出する。この圧力ソース密度は、「構造物モデル表面の解析対象点における表面圧力変動」に対する「所定領域内の空間節点における流れ場」の寄与度を示す指標である。ここで、所定領域とは、構造物モデルの周囲の流れ場の一部であり、解析対象点における表面圧力変動に影響を与える可能性のある領域として規定された領域である。このため、圧力ソース密度を指標として用いることにより、表面圧力変動に対する流れ場の寄与度を空間節点毎に算出できる。従って、本発明装置によれば、所定領域内のどの位置における流れ場が解析対象点における構造物モデルの表面圧力変動に大きく寄与しているかを適切に特定することが可能となる。なお、表面圧力変動は風切音と強く相関しており、表面圧力変動が大きいときは風切音が大きくなる。このため、表面圧力変動を解析することは、広義には、風切音を解析することと同義である。加えて、圧力変動は、流体流動圧力変動と音圧から構成されている。本明細書における表面圧力変動は、これらのうち流体流動圧力変動と同義である。
【０００８】
本発明装置の一側面では、
前記構造物モデル（２０）表面の前記解析対象点（ｘ）についての前記圧力ソース密度（ＰＤＳ）を前記所定領域（２１）で空間積分した値は、前記解析対象点（ｘ）における前記表面圧力変動の関数（ｐθ）とその複素共役関数（ｐθ＊）との積を着目角周波数帯（[θｌ，θｈ]）で積分した値の近似値である。
【０００９】
この構成によれば、圧力ソース密度の挙動は、構造物モデルの表面圧力変動の挙動と精度良く一致する。即ち、圧力ソース密度は、表面圧力変動に対する流れ場の寄与度を表す指標として高い信頼性を有する。このため、構造物の表面圧力変動に大きく寄与している構造物周囲の流れ場の位置を精度良く特定することができる。
【００１０】
本発明装置の一側面は、
更に、前記圧力ソース密度（ＰＤＳ）を構成する複数のパラメータのうち、当該圧力ソース密度（ＰＤＳ）への寄与が相対的に大きいパラメータである原因パラメータを特定する原因パラメータ特定手段を備える。
【００１１】
この構成によれば、ＰＤＳが高い原因となっているパラメータを容易に把握できるため、構造物モデルの形状をより効率的に検討、変更することができる。
【００１２】
本発明装置の一側面では、
前記複数のパラメータは、前記平均流速（Ｕ−）、前記平均渦度（Ω−）及び前記乱れ流速（ｕ〜θ）である。
【００１３】
本発明装置の一側面は、
更に、前記圧力ソース密度（ＰＤＳ）が算出された複数の空間節点（ｚ）の中から、外部から入力される値を有する圧力ソース密度（ＰＤＳ）に対応している複数の空間節点（ｚ）を抽出し、当該抽出された複数の空間節点（ｚ）を画像処理して等値面（２８、３０、３２）を作成し、当該等値面（２８、３０、３２）を可視化する画像処理手段（ステップ８２０）を備える。
【００１４】
この構成によれば、ＰＤＳの等値面が表示されるため、作業者は、入力値を適宜設定することにより、所望のＰＤＳ値を有する等値面を視認することができる。この結果、所望のＰＤＳ値に対応している空間節点を効率的に選択することができ、構造物の表面圧力変動に大きく寄与している構造物周囲の流れ場の位置を効率的に特定することができる。
【００１５】
本発明装置の一側面は、
更に、前記圧力ソース密度算出手段（ステップ８１４）によって算出された複数の圧力ソース密度（ＰＤＳ）の最大値に対応している空間節点（ｚ）を抽出する空間節点抽出手段を備える。
【００１６】
この構成によれば、最大のＰＤＳ値に対応している空間節点を作業者が選択する必要がなくなるため、構造物の表面圧力変動に大きく寄与している構造物周囲の流れ場の位置をより効率的に特定することができる。なお、「圧力ソース密度の最大値」とは、広義には、「圧力ソース密度がほぼ最大である値」を意味する。
【００１７】
本発明装置の一側面では、
前記非定常ＣＦＤ計算手段（ステップ８０２、８０６、８０８）は、更に、前記着目角周波数帯において、前記構造物モデル（２０）の表面における前記表面圧力変動（ｐθ）を、当該表面上の節点である表面節点（ｘ）毎に算出する。
【００１８】
この構成によれば、表面圧力変動に基づいて解析対象点を選択することが可能になる。表面圧力変動は風切音と強い相関を有している。このため、風切音をより適切に解析することができる。
【００１９】
本発明装置の一側面では、
前記圧力ソース密度（ＰＤＳ）の計算式は、下記式によって定義される。
【数１】
2019191110.tif 000003
【００２０】
加えて、本明細書は、移動中の構造物の表面に発生する風切音を解析する新規な風切音解析方法を開示する。
この風切音解析方法は、
構造物をモデル化した構造物モデル（２０）を移動させる非定常ＣＦＤシミュレーションを実行して、前記構造物モデル（２０）の周囲の流れ場のうち所定領域（２１）内の流れ場の所定時間における平均流速（Ｕ−）及び平均渦度（Ω−）を、当該所定領域（２１）内の節点である空間節点（ｚ）毎に算出するとともに、風切音の解析対象となる角周波数帯である着目角周波数帯（[θｌ，θｈ]）において、前記所定領域（２１）内における乱れ流速の振幅（ｕ〜θ）に基づく値を、前記空間節点（ｚ）毎に算出する非定常ＣＦＤ計算工程（ステップ８０２、８０６、８０８）と、
前記構造物モデル（２０）表面の風切音の解析対象点（ｘ）における圧力変動の振幅である表面圧力変動（ｐθ）に対する、前記所定領域（２１）内の空間節点（ｚ）における流れ場の寄与度を示す指標である圧力ソース密度（ＰＤＳ）を、前記非定常ＣＦＤ計算工程（ステップ８０２、８０６、８０８）によって算出された前記平均流速（Ｕ−）と、前記平均渦度（Ω−）と、前記乱れ流速の振幅（ｕ〜θ）に基づく値と、に基づいて算出する圧力ソース密度算出工程（ステップ８１４）と、
を備える。
【００２１】
この風切音解析方法によれば、所定領域内のどの位置における流れ場が解析対象点における構造物モデルの表面圧力変動に大きく寄与しているかを適切に特定することが可能となる。
【００２２】
なお、上記説明においては、発明の理解を助けるために、実施形態に対応する発明の構成に対して、実施形態で用いた符号を括弧書きで添えているが、発明の各構成要件は前記符号によって規定される実施形態に限定されるものではない。
【図面の簡単な説明】
【００２３】
【図１】本発明の実施形態に係る風切音解析装置（以下、「本実施装置」又は「実施装置」と称する。）が有する計算装置を示すブロック図である。
【図２Ａ】非定常ＣＦＤシミュレーションの対象となる車両モデル及びフロードメインの模式図である。
【図２Ｂ】時刻ｔ１における圧力変動の時刻歴データを車両モデルにプロットした図である。
【図３Ａ】車両モデルの解析表面における表面圧力変動の分布図である。
【図３Ｂ】表面圧力変動（風切音）の解析対象となる観測点の位置を示す図である。
【図４Ａ】着目評価点の選択方法を説明するために用いる図であり、ＰＤＳの参照値が１０３である等値面を示す図である。
【図４Ｂ】ＰＤＳの参照値が１０４である等値面を示す図である。
【図４Ｃ】ＰＤＳの参照値が１０５である等値面を示す図である。
【図５Ａ】原因パラメータ（乱れ流速）の分布図が作成される平面の位置を示す図である。
【図５Ｂ】乱れ流速の流速レベルの分布図である。
【図６】ＰＤＳ計算式の導出方法を説明するために用いる図であり、剛体モデル及びフロードメインの模式図である。
【図７】剛体モデル及びフロードメインの模式図である。
【図８】風切音解析方法の手順を示すフローチャートである（その１）。
【図９】風切音解析方法の手順を示すフローチャートである（その２）。
【図１０Ａ】車両モデルの形状を変更したときの表面圧力変動（風切音）の改善効果を説明するために用いる図であり、形状変更前の右側サイドミラーの概形を示す。
【図１０Ｂ】形状変更後の右側サイドミラーの概形を示す。
【図１１Ａ】形状変更前の右側サイドミラーの表面圧力変動の分布図である。
【図１１Ｂ】形状変更前の右側サイドミラーの表面圧力変動の分布図である。
【図１２Ａ】指標としてのＰＤＳの信頼性及び汎用性を説明するために用いる図であり、車両モデルの模式図である。
【図１２Ｂ】フォアステップモデルの模式図である。
【図１３Ａ】車両モデルの右側フロントサイドガラスにおける、ＰＤＳを用いて予測される表面圧力変動（予測表面圧力変動）の分布図である。
【図１３Ｂ】車両モデルの右側フロントサイドガラスにおける、ソフトウエアの非定常ＣＦＤシミュレーションにより算出される表面圧力変動の分布図である。
【図１４Ａ】フォアステップモデルのステップの上面における、ＰＤＳを用いて予測される表面圧力変動（予測表面圧力変動）の分布図である。
【図１４Ｂ】フォアステップモデルのステップの上面における、ソフトウエアの非定常ＣＦＤシミュレーションにより算出される表面圧力変動の分布図である。
【図１５】表面圧力変動に対する予測表面圧力変動の誤差を車両モデル及びフォアステップモデルについて示したグラフである。
【発明を実施するための形態】
【００２４】
＜本実施装置の概要＞
まず、実施形態に係る風切音解析装置（以下、「本実施装置」とも称する。）の概要を説明する。車両走行時に発生する風切音は、車両表面における圧力変動の振幅（表面圧力変動）と強く相関している。表面圧力変動は、車両周囲の流れ場（流速及び渦度）の影響を受けて変化する。本実施装置には、車両周囲の流れ場と表面圧力変動との関係を定量的に示す指標（別言すれば、「車両周囲の任意の位置における流れ場（の状態）」が「車両表面の任意の位置における表面圧力変動」にどれくらい寄与しているかを示す指標）の計算式が予め格納されている。指標の値が大きいほど寄与度は高くなる。本実施装置は、非定常ＣＦＤ（数値流体力学）計算により、走行中の車両モデルの表面圧力変動を、車両モデル表面の位置毎、及び、周波数帯毎、に算出する。作業者は、当該算出結果に基づいて、表面圧力変動が比較的に大きい位置（別言すれば、大きな風切音の発生原因となっている位置）を選択する。本実施装置は、上述した指標の計算式に基づいて、「当該選択された位置における表面圧力変動に対する流れ場の寄与度を示す指標」を流れ場の位置毎に算出する。なお、車両モデルは「構造物モデル」の一例に相当する。
【００２５】
作業者は、算出結果に基づいて、指標の値が比較的に大きい流れ場の位置（別言すれば、上記選択された位置における表面圧力変動に対する寄与度が比較的に高い位置）を選択し、当該選択された流れ場の位置における指標を構成する各物理量（典型的には、平均流速、乱れ流速及び平均渦度）のうち、指標の値が大きい原因となっている物理量を特定する。作業者は、当該特定された物理量の値に基づいて、上記選択された流れ場の位置における指標値が低減するように車両モデルの形状を検討、変更する。この構成によれば、指標値が低減するように車両モデルの形状を変更することにより、「当該指標値を有する位置の流れ場の、上記選択された位置における表面圧力変動への寄与度」を低減することができる。このため、当該選択された位置における表面圧力変動を低減でき、結果として、風切音の低減を実現することができる。
【００２６】
＜本実施装置の具体的な構成＞
以下、図面を参照して本実施装置について具体的に説明する。図１に示すように、本実施装置は計算装置１０を備える。計算装置１０は、入力部１２、演算部１４及び出力部１６を有する。出力部１６は、作業者が視認可能な位置に表示画面１８を有する。入力部１２は、風切音の解析対象となる車両の三次元モデルのデータ、解析範囲を示すフロードメインのデータ、及び、非定常ＣＦＤシミュレーションに必要なデータ等を入力する。加えて、入力部１２は、処理の途中で作業者により選択される座標のデータ、指標の参照値及び原因パラメータ等を入力する（後述）。
【００２７】
図２Ａは、車両の三次元モデル２０（以下、単に「車両モデル２０」と称する。）及びフロードメイン２１を模式的に示す図である。図２Ａに示すように、車両モデル２０の空間座標系はｅ１軸、ｅ２軸及びｅ３軸からなり、当該座標系は、車両モデル２０の前後方向における前方が＋ｅ１方向に一致し、且つ、左右方向における左方向が＋ｅ２方向に一致するように設定される。加えて、原点は、車両モデル２０の基準点（例えば、重心）から所定の方向に所定の距離（ゼロ値を含む）だけ離れた位置に設定される。
【００２８】
フロードメイン２１は、車両モデル２０の表面及びその周囲の空間のうち、非定常ＣＦＤ計算が実施される範囲を示す領域である。本実施形態では、車両モデル２０の右側サイドミラー２２及び右側フロントピラー２４の近傍で発生する風切音の評価が行われる。このため、フロードメイン２１は、「これらの部材２２、２４の表面を含み、且つ、車両モデル２０の周囲の流れ場のうち、これらの部材２２、２４の近傍における表面圧力変動に影響を与える可能性のある部分の流れ場を含んだ領域」として設定される。具体的には、フロードメイン２１は、右側サイドミラー２２、右側フロントピラー２４、フロントサイドガラス及びフロントガラスの右半分等を含む略直方体形状の空間から車両モデル２０と重複している部分がくりぬかれた形状を有する。フロードメイン２１は、風切音の解析が行われる車両モデル２０の位置に応じて自由に設定可能である。なお、フロードメイン２１は、「所定領域」の一例に相当する。
【００２９】
入力部１２にて入力されたこれらのデータは、演算部１４のＲＡＭ（後述）に格納される。なお、以下では、非定常ＣＦＤ計算を単に「ＣＦＤ計算」とも称する。
【００３０】
図１に戻って説明を続ける。演算部１４は、ＣＰＵ、ＲＯＭ、ＲＡＭ等からなるマイクロコンピュータを主要構成部品として有する。演算部１４は、各種インターフェースを有しており、このインターフェースを介して入力部１２及び出力部１６と信号を入出力可能に接続されている。
【００３１】
演算部１４は、ＲＡＭに格納されている上述したデータを用いて車両モデル２０を走行させる非定常ＣＦＤシミュレーションを実行する。ＣＦＤ計算処理は、時刻歴データの算出処理、平均化処理及び高速フーリエ変換処理を含む。この計算処理は、フロードメイン２１によって区画された範囲で行われる。以下、具体的に説明する。
【００３２】
（ＣＦＤ計算処理）
Ａ．時刻歴データの算出処理
演算部１４は、車両モデル２０の表面のうちフロードメイン２１と接している部分の表面（以下、「解析表面」とも称する。）における各表面節点ｘ（以下、単に「節点ｘ」と称する。）について、圧力変動の時刻歴データｐ（ｘ，ｔ）を所定時間に亘って算出する（ｔ：時刻）。解析表面は、右側サイドミラー２２、右側フロントピラー２４、フロントサイドガラス及びフロントガラスの右半分等を含む面である。なお、ルーフの一部もフロードメイン２１と接しているが、解析表面からは除外している。演算部１４は、圧力変動の時刻歴データｐ（ｘ，ｔ）を節点ｘの座標及び時刻ｔに関連付けてそのＲＡＭに格納する。
【００３３】
図２Ｂは、時刻ｔ１における圧力変動の時刻歴データｐ（ｘ，ｔ１）（即ち、時刻ｔ１のときの各節点ｘにおける圧力変動の振幅）を車両モデル２０にプロットした図である。圧力変動の振幅が大きいほど濃色で表示されている。図２Ｂによると、時刻ｔ１のときは、右側サイドミラー２２の車体側の一部、及び、フロントガラスの下部中央において圧力変動が比較的に高くなっている。演算部１４は、節点ｘ毎に圧力変動の振幅と色データとを対応づけたデータを所定時間における各時刻ｔについて作成し、作成したデータを出力部１６の表示画面１８にアニメーション表示させる表示指令を出力部１６に送信する。出力部１６は、当該表示指令を受信すると、その表示画面１８に当該データをアニメーション表示する。これにより、作業者は、所定時間における圧力変動の時間変化を把握することができる。
【００３４】
加えて、演算部１４は、フロードメイン２１内の各空間節点ｚについて、流れ場の流速の時刻歴データｕ（ｚ，ｔ）及び渦度の時刻歴データω（ｚ，ｔ）を所定時間に亘ってそれぞれ算出する。以下では、空間節点ｚを単に「節点ｚ」と称する。なお、フロードメイン２１内の節点ｚの位置は、節点ｚが解析表面上に存在することがないように（即ち、節点ｚが解析表面から離間するように）予め設定されている。演算部１４は、流速及び渦度の時刻歴データｕ（ｚ，ｔ）及びω（ｚ，ｔ）を、節点ｚの座標及び時刻ｔに関連付けてそのＲＡＭに格納する。
【００３５】
２．平均化処理
演算部１４は、ＲＡＭに格納されている流速及び渦度の時刻歴データｕ（ｚ，ｔ）及びω（ｚ，ｔ）をそれぞれ平均化処理して、フロードメイン２１内の各節点ｚについて、平均流速Ｕ−（ｚ）及び平均渦度Ω−（ｚ）を算出する。なお、Ｕ−（ｚ）及びΩ−（ｚ）の「−」は何れも平均を表す表記であり、後述する式中の文字Ｕ及びΩの上に付されたバーに替わるものである。演算部１４は、平均流速Ｕ−（ｚ）及び平均渦度Ω−（ｚ）を、節点ｚの座標に関連付けてそのＲＡＭに格納する。平均流速Ｕ−（ｚ）及び平均渦度Ω−（ｚ）は、指標を算出する際に用いられる（後述）。
【００３６】
３．高速フーリエ変換処理
演算部１４は、ＲＡＭに格納されている圧力変動の時刻歴データｐ（ｘ，ｔ）を高速フーリエ変換（FFT: Fast Fourier Transform）処理して、解析表面の各節点ｘについて、周波数（角周波数）毎に圧力変動の振幅ｐθ（ｘ）（即ち、表面圧力変動ｐθ（ｘ））を算出する（θ：角周波数。θ＝２πｆ）。表面圧力変動ｐθ（ｘ）は、風切音の評価対象となる周波数帯において算出される（以下、当該周波数帯毎を「着目周波数帯」とも称する。）。着目周波数帯は作業者によって選択可能であり、例えば、中心周波数ｆｍが５００Ｈｚ、１ｋＨｚ、２ｋＨｚ又は４ｋＨｚの周波数帯が選択され得る。以下では、中心周波数ｆｍがｎＨｚの周波数帯を、単に、「ｎＨｚの周波数帯」と称する。
【００３７】
例えば、着目周波数帯として１ｋＨｚ（下限周波数ｆｌ＝７１０Ｈｚ、上限周波数ｆｈ＝１４２０Ｈｚ）が選択された場合、表面圧力変動ｐθ（ｘ）は、７１０Ｈｚから１４２０Ｈｚの範囲内で、解析表面の各節点ｘについて算出される。演算部１４は、表面圧力変動ｐθ（ｘ）を、節点ｘの座標及び角周波数θに関連付けてそのＲＡＭに格納する。表面圧力変動ｐθ（ｘ）は、着目周波数帯の中心周波数ｆｍにおける表面圧力変動ｐθ（ｘ）の分布図を作成する際に用いられる（後述）。
【００３８】
加えて、演算部１４は、ＲＡＭに格納されている流速の時刻歴データｕ（ｚ，ｔ）を高速フーリエ変換処理して、フロードメイン２１内の各節点ｚについて、着目周波数帯における乱れ流速（厳密には、乱れ流速の振幅）ｕ〜θ（ｚ）の自己相関関数を算出する。なお、ｕ〜θ（ｚ）の「〜」は、流速ｕを平均化処理したときの乱れ成分を表す表記であり、後述する式中の文字ｕの上に付されたチルダに替わるものである。加えて、任意の関数Ｆθの自己相関関数は、下記式（１）によって定義される。
【数２】
2019191110.tif 000004
【００３９】
ここで、上記式（１）中のθｌ及びθｈは、着目角周波数帯の下限角周波数及び上限角周波数をそれぞれ表し、Ｆθ＊は、関数Ｆθの複素共役を表す。即ち、本明細書では、自己相関関数は、「任意の関数Ｆθとその複素共役関数Ｆθ＊との積（後述）を着目角周波数帯において積分した値」として定義される。
【００４０】
演算部１４は、乱れ流速ｕ〜θ（ｚ）の自己相関関数を、節点ｚの座標に関連付けてそのＲＡＭに格納する。乱れ流速ｕ〜θ（ｚ）の自己相関関数は、指標を算出する際に用いられる（後述）。演算部１４は、以上の処理をＣＦＤ計算として実行する。
【００４１】
（表面圧力変動表示指令）
演算部１４は、ＣＦＤ計算を終了すると、ＲＡＭに格納されている「各節点ｘにおける着目周波数帯の表面圧力変動ｐθ（ｘ）」の中から「各節点ｘにおける中心周波数ｆｍの表面圧力変動ｐθ（ｘ）」を抽出する。そして、「抽出した表面圧力変動ｐθ（ｘ）にその大きさに応じた色データを対応づけたデータ」を各節点ｘについて作成し、作成したデータを表示画面１８に画像表示させる表示指令を出力部１６に送信する。表面圧力変動ｐθ（ｘ）が大きいほど濃色の色データが対応づけられる。出力部１６は、当該表示指令を受信すると、その表示画面１８に当該データを画像表示する。
【００４２】
図３Ａは、当該表示指令に基づいて表示画面１８に画像表示された表面圧力変動ｐθ（ｘ）の分布図の一例である。着目周波数帯は１ｋＨｚ（即ち、ｆｍ＝１ｋＨｚ）に設定されている。この例では、解析表面のうち右側サイドミラー２２の近傍の部分のみが表示されている。作業者は、この分布図に基づいて表面圧力変動ｐθ（ｘ）が比較的に大きい節点ｘを選択する。図３Ａによると、着目周波数帯が１ｋＨｚのときは右側サイドミラー２２近傍の領域２６において表面圧力変動ｐθ（ｘ）が顕著であるため、作業者は、領域２６内の節点ｘの中から任意の節点を選択する。作業者は、以後、当該節点における表面圧力変動の低減を目的として当該節点についての指標を算出し（後述）、車両モデル２０の形状を変更して当該節点における表面圧力変動の挙動を観測する。このため、以下では、上記選択された節点を「観測点ｘ」とも称する（図３Ｂ参照）。観測点ｘは、「解析表面のうち、着目周波数帯における大きな風切音の発生原因となっている位置」を表す。作業者により観測点ｘが選択されると、入力部１２は、当該観測点ｘの座標データを入力する。入力された座標データは演算部１４のＲＡＭに格納される。なお、観測点ｘは、「解析対象点」の一例に相当する。
【００４３】
（ＰＤＳの算出）
演算部１４のＲＡＭには、フロードメイン２１内の流れ場と観測点ｘにおける表面圧力変動との関係を定量的に示す指標（別言すれば、「フロードメイン２１内の任意の節点ｚにおける流れ場（の状態）」が「観測点ｘにおける表面圧力変動ｐθ（ｘ）」にどれくらい寄与しているかを示す指標）の計算式が予め格納されている。この指標は、フロードメイン２１内の各節点ｚについて算出可能である。この指標をフロードメイン２１内の全節点ｚで空間積分した値は、「観測点ｘにおける表面圧力変動ｐθ（ｘ）の自己相関関数」の近似値である（後述）。このため、各節点ｚにおける指標の値は、フロードメイン２１の単位体積当たりの「観測点ｘにおける表面圧力変動ｐθ（ｘ）の自己相関関数」の値を表していると解釈することができる。従って、以下では、各節点ｚにおける指標を「ＰＤＳ（Pressure Density Source（圧力ソース密度））」とも称する。ＰＤＳ（ｘ，ｚ）は、「フロードメイン２１内の任意の節点ｚにおける流れ場」の「観測点ｘにおける表面圧力変動ｐθ（ｘ）」への寄与度を示す指標であり、その値が大きいほど寄与度が高いことを示す。なお、ＰＤＳの計算式の表記及び導出方法については後述する。
【００４４】
演算部１４は、ＲＡＭに格納されているＰＤＳの計算式に基づいて、フロードメイン２１内の各節点ｚについてＰＤＳ（ｘ，ｚ）を算出する。ＰＤＳ（ｘ，ｚ）の値は、節点ｚ毎に評価されるため、以下では、節点ｚを「評価点ｚ」とも称する。演算部１４は、ＰＤＳ（ｘ，ｚ）の値を観測点ｘの座標及び評価点ｚの座標に関連付けてそのＲＡＭに格納する。ＰＤＳは、作業者により設定される値を有するＰＤＳの等値面を作成する際に用いられる（後述）。
【００４５】
（ＰＤＳ等値面表示指令）
演算部１４は、ＰＤＳの算出を終了すると、ＲＡＭに格納されているＰＤＳ（ｘ，ｚ）の中から、作業者により設定された任意の値（以下、「参照値」とも称する。後述。）を有するＰＤＳ（ｘ，ｚ）に関連付けられた評価点ｚの座標を抽出する。そして、抽出された複数の評価点ｚ同士を接続した等値面を表す等値面データを作成し、作成した等値面データを表示画面１８に画像表示させる表示指令を出力部１６に送信する。出力部１６は、当該表示指令を受信すると、その表示画面１８に「参照値を有するＰＤＳの等値面データ」を画像表示する。即ち、演算部１４は、ＰＤＳの算出を終了すると、その旨を示すメッセージを表示画面１８に表示させ、作業者にＰＤＳの算出が終了したことを認識させる。これにより、作業者により任意の参照値が入力される。演算部１４は、上記表示指令の送信処理を、作業者により参照値が入力される毎に繰り返し行う。
【００４６】
図４Ａ〜図４Ｃは、参照値が互いに異なるＰＤＳの等値面を示す。図４Ａは、参照値が１０３である等値面２８を示し、図４Ｂは、参照値が１０４である等値面３０を示し、図４Ｃは、参照値が１０５である等値面３２を示す。着目周波数帯は何れも１ｋＨｚに設定されている。図４Ａ〜図４Ｃに示すように、等値面２８、３０、３２は何れも閉曲面である。加えて、各等値面２８、３０、３２は、複数の分離した等値面を有している。ＰＤＳ（ｘ，ｚ）は、観測点ｘ（図３Ｂ参照）から評価点ｚまでの距離ｒが減少するほど増加する性質を有する。即ち、評価点ｚが観測点ｘに接近するほどＰＤＳの値が増加する。このため、図４Ａ〜図４Ｃに示すように、参照値が増加するほど等値面２８、３０、３２の表面積が減少している。加えて、等値面３０は等値面２８の内部に含まれており、等値面３２は等値面３０の内部に含まれている。即ち、ある参照値を有する等値面の内部には、当該参照値より大きい参照値を有する等値面が包含されている。このため、参照値が徐々に増加すると等値面の表面積が徐々に減少するとともに等値面の個数も減少していき、ある時点で等値面の個数が２個となる。
【００４７】
図４Ｃにその一例を示す。図４Ｃに示すように、等値面３２は２個の等値面３２ａ、３２ｂを有している。図３Ｂを参照すると、等値面３２ａのほうが等値面３２ｂよりも観測点ｘに近い（即ち、距離ｒが短い）。このため、等値面３２ａのＰＤＳは、等値面３２ｂのＰＤＳよりも距離ｒの影響を強く受けている（即ち、距離ｒのＰＤＳ値への寄与度が高い）と言える。ここで、ＰＤＳ（ｘ，ｚ）は、パラメータとして以下の物理量、即ち、平均流速Ｕ−（ｚ）、乱れ流速ｕ〜θ（ｚ）、及び、平均渦度Ω−（ｚ）を有する。ＰＤＳの値は、これらの物理量が大きいほど大きくなる。このため、「等値面３２ａよりも距離ｒが長い等値面３２ｂ」が等値面３２ａと同一のＰＤＳ値を有しているのは、等値面３２ｂのＰＤＳが等値面３２ａのＰＤＳよりも上記物理量（Ｕ−（ｚ）、ｕ〜θ（ｚ）、Ω−（ｚ））の少なくとも１つの影響を強く受けている（即ち、上記物理量のＰＤＳ値への寄与度が高い）ためであると解釈することができる。
【００４８】
従って、作業者は、「距離ｒよりも上記物理量の影響を相対的に強く受けているＰＤＳ」を有する等値面３２ｂ上の任意の評価点ｚを、ＰＤＳ値が比較的に大きい（ほぼ最大である）評価点ｚとして選択する。以下、当該評価点ｚを「着目評価点ｚ」とも称する（図５Ａ参照）。着目評価点ｚは、「フロードメイン２１内の流れ場のうち、観測点ｘにおける表面圧力変動に対する寄与度が比較的に高い（ほぼ最大である）位置」を表す。作業者により着目評価点ｚが選択されると、入力部１２は、当該着目評価点ｚの座標データを入力する。入力された座標データは演算部１４のＲＡＭに格納される。
【００４９】
（ＰＤＳパラメータ表示指令）
演算部１４は、ＲＡＭに格納されている着目評価点ｚにおけるＰＤＳ（ｘ，ｚ）を構成している複数のパラメータ（物理量）の中から平均流速Ｕ−（ｚ）、乱れ流速ｕ〜θ（ｚ）、及び、平均渦度Ω−（ｚ）を選択し、これらの数値データを表示画面１８に表示させる表示指令を出力部１６に送信する。出力部１６は、当該表示指令を受信すると、その表示画面１８に選択されたパラメータの数値データを表示する。作業者は、表示された数値データに基づいて、これらのパラメータの中でＰＤＳ値への寄与が相対的に大きいパラメータ（即ち、ＰＤＳ値が大きい原因となっているパラメータ）を特定する。このとき、必要に応じて他の評価点ｚにおけるＰＤＳを構成するパラメータの数値データが表示されてもよい。作業者は、これらの数値データを比較しながらパラメータを特定してもよい。作業者によりパラメータが特定されると、入力部１２は、当該パラメータを入力する。入力されたパラメータは演算部１４のＲＡＭに格納される。以下では、このパラメータを「原因パラメータ」とも称する。
【００５０】
（物理量分布図表示指令）
演算部１４は、ＲＡＭに格納されている評価点ｚ群の座標データから、「着目評価点ｚを通りｅ１ｅ２平面（図２Ａ参照）に平行な平面（図５Ａの破線Ｌ参照）」上に位置する評価点ｚ群を抽出する。そして、抽出された評価点ｚ群に関連付けられた原因パラメータの分布図を表す分布図データを作成し、作成した分布図データを表示画面１８に画像表示させる表示指令を出力部１６に送信する。出力部１６は、当該表示指令を受信すると、その表示画面１８に当該データを画像表示する。
【００５１】
図５Ｂは、原因パラメータとして乱れ流速ｕ〜θ（ｚ）が特定された場合における乱れ流速ｕ〜θ（ｚ）の流速レベル（ｄＢ）の分布図を示す。流速レベルが大きいほど濃色で表示されている。図５Ａに示すように、この例では着目評価点ｚが右側サイドミラー２２の近傍に位置しているため、分布図の範囲は、右側サイドミラー２２及びその周囲のみが設定されている。なお、右側サイドミラー２２の位置には、車両モデル２０を平面視したときの右側サイドミラー２２の画像がはめ込み表示されている。着目評価点ｚは、図５Ｂの領域３４に含まれている。
【００５２】
作業者は、分布図に基づいて着目評価点ｚにおいて原因パラメータ（本例では乱れ流速ｕ〜θ（ｚ））が大きくなっている原因を分析する。図５Ｂの例では、右側サイドミラー２２の近傍においては、領域３４の乱れ流速ｕ〜θ（ｚ）が顕著である。これは、車両走行時において空気の流れが領域３４で剥離しているためと考えられる。従って、作業者は、空気の流れが領域３４で剥離する可能性を低減できるような右側サイドミラー２２の形状を検討してその形状を変更する。形状変更後の表面圧力変動ｐθ（ｘ）の具体的な改善効果については後述する。
【００５３】
出力部１６は、演算部１４からの表示指令を受信して表示画面１８に各種データを表示する。詳細は上述したとおりである。以上が、本実施装置の具体的な構成についての
説明である。
【００５４】
＜ＰＤＳ計算式の導出＞
続いて、上述したＰＤＳの計算式の導出方法を剛体モデルを用いて説明する。ＰＤＳは、公知のパウエル（Powell）式に基づいて算出される「剛体モデルにおける表面圧力変動とその周囲の流れ場との関係式」を時間ドメインから周波数ドメインに変換した式に基づいて定義される。但し、上記パウエル式は、ＰＤＳ計算式の導出に適した形式で表現される。このため、以下では、まず、図６を参照して剛体モデル及びフロードメインについて説明し、続いて、ＰＤＳ計算式の導出に適した形式のパウエル式の導出について説明する。次いで、時間ドメインにおける表面圧力変動と流れ場との関係式の導出について説明し、その後、ＰＤＳ計算式の導出について説明する。
【００５５】
（剛体モデル及びフロードメイン）
図６は、剛体モデル３６及びフロードメイン３８の模式図を示す。剛体モデル３６（以下、単に「剛体３６」とも称する。）は、その表面Ｓ上に観測点ｘを有する。フロードメイン３８は、剛体３６の周囲の流れ場のうち、観測点ｘにおける表面圧力変動に影響を与える可能性のある部分の流れ場である。フロードメイン３８は、その内部に複数の評価点ｚを有する。ＰＤＳ（ｘ，ｚ）は、観測点ｘにおける表面圧力変動に対する評価点ｚにおける流れ場の寄与度を表す指標である。なお、図６中のその他の文字については後述する。
【００５６】
（パウエル式の導出）
下記式（２）は、体積力がゼロの場合のナビエ−ストークス方程式（Navier-Stokes equation）を示す。
【数３】
2019191110.tif 000005
ここで、ρは流れ場の密度であり、ｕは流れ場の流速（ベクトル量）であり、Ｄ／Ｄｔは物質微分であり、ｐは剛体３６の表面における圧力（圧力変動、表面圧力変動）であり、ηは流れ場の粘性係数である。
【００５７】
上記式（２）の物質微分を偏微分を用いて書き直し、その後、以下の関係式；
【数４】
2019191110.tif 000006
を用いて式変形すると、下記式（３）が得られる。なお、上記の関係式中のＡは任意のベクトルを表す。
【数５】
2019191110.tif 000007
ここで、νは流れ場の動粘性係数であり、ν＝η／ρである。
【００５８】
続いて、上記式（３）を以下の関係式；
【数６】
2019191110.tif 000008
を用いて式変形すると、下記式（４）が得られる。なお、ωは流れ場の渦度（ベクトル量）であり、Ｂはエンタルピーである。
【数７】
2019191110.tif 000009
【００５９】
本実施形態の流れ場は非圧縮性流体又は微圧縮性流体であるが、この場合、以下の関係式；
【数８】
2019191110.tif 000010
が成立するため、上記式（４）は下記式（５）によって近似できる。
【数９】
2019191110.tif 000011
【００６０】
上記式（５）の両辺に流れ場の密度ρを乗算して発散をとると、下記式（６）が得られる。
【数１０】
2019191110.tif 000012
なお、上記式（６）中の「×」は乗算記号を表す。本明細書では基本的に乗算記号は省略して記載するが、明示したほうが理解し易い場合は例外的に記載する。
【００６１】
ここで、流れ場の流速ｕと流れ場の密度ρとの間には、下記式（７）に示す関係式が成立する。
【数１１】
2019191110.tif 000013
「流れ場の密度ρ」と「流れ場の流速ｕの時間についての偏微分」との積の発散（即ち、下記式（８）の左辺）を上記式（７）を用いて式変形すると、下記式（８）の右辺のようになる。
【数１２】
2019191110.tif 000014
ここで、ｃは流れ場の音速である。表面圧力変動ｐ、流れ場の密度ρ、及び、流れ場の音速ｃとの間には以下の関係が成立する。
【数１３】
2019191110.tif 000015
【００６２】
ところで、本実施形態の流れ場は高レイノルズ数流れ（即ち、流速が比較的に大きい流れであり、例えば、１００km/h〜１４０km/h。）である。高レイノルズ数流れでは、以下の関係式；
【数１４】
2019191110.tif 000016
が成立する。なお、レイノルズ数Ｒｅは、Ｒｅ＝ｕＬ／νとして定義され、Ｌは流れ場の代表長さを表す。上記式（８）の右辺の括弧内の項を以下に示すように式変形して上記関係式を適用すると、式（８）の右辺の括弧内の項は下記式（９）によって近似される。
【数１５】
2019191110.tif 000017
【００６３】
上記式（６）を上記式（８）及び式（９）を用いて式変形（近似）すると、下記式（１０）が得られる。
【数１６】
2019191110.tif 000018
【００６４】
ここで、本実施形態の流れ場は低マッハ数流れ（即ち、音速に対する流体の流速の比が比較的に低い流れ場）である。低マッハ数流れでは以下の関係式；
【数１７】
2019191110.tif 000019
が成立するため、上記式（１０）をこれらの関係式を用いて式変形すると、下記式（１１）に示すパウエル式が得られる。なお、上記の関係式のｃ０は静止状態のときの流れ場の音速であり、ρ０は静止状態のときの流れ場の密度である。
【数１８】
2019191110.tif 000020
【００６５】
（時間ドメインにおける表面圧力変動と流れ場との関係式の導出）
グリーン（Green）関数Ｇは、下記式（１２）の関係式を満たす。
【数１９】
2019191110.tif 000021
ここで、観測点ｘ及び評価点ｚは上述したとおりであり（図６参照）、τは過去の時刻である。なお、δ（ｘ）はディラックデルタ関数である。加えて、上記式（１２）中のグリーン関数Ｇは、自由空間においては下記式（１３）のように表される。
【数２０】
2019191110.tif 000022
【００６６】
物体境界条件を考慮すると、上記式（１１）に示すパウエル式は、下記式（１４）のように表される。
【数２１】
2019191110.tif 000023
ここで、Ｈ（ｚ）はヘビサイド関数（Heaviside Unit Function）であり、下記式（１５）によって定義される。
【数２２】
2019191110.tif 000024
ここで、図６に示すように、Ｖはフロードメイン３８を表すかたまりであり、Ｓ＋は剛体３６の表面Ｓを取り囲む仮想表面である。上記式（１５）中の「z in V」とは、評価点ｚがフロードメイン３８の内部に位置している場合を示し、「z inside S+」とは、評価点ｚが仮想表面Ｓ＋の内部に位置している場合を示す。
【００６７】
このヘビサイド関数Ｈ（ｚ）は、下記式（１６）の積分方程式を満たす。
【数２３】
2019191110.tif 000025
ここで、式中の（・）は任意の値との内積を表す。加えて、図６に示すように、剛体３６は、その表面Ｓのうち観測点ｘとは異なる位置に評価点ｙを有する。式中のｎは、当該評価点ｙの法線ベクトルｎ（ｙ）を表し、ｄＳは、評価点ｙ近傍の表面面積ｄＳ（ｙ）を表す。
【００６８】
上記式（１４）の左辺の１項目、２項目及び右辺は、それぞれ下記式(14-1)、(14-2)及び(14-3)のように変形できる。
【数２４】
2019191110.tif 000026
【００６９】
上記式(14-1)〜(14-3)を上記式（１４）に代入すると、下記式（１７）が得られる。
【数２５】
2019191110.tif 000027
【００７０】
上記式（１７）の右辺の２項目（−符号は省略）は、下記式(17-1)のように式変形できる。
【数２６】
2019191110.tif 000028
【００７１】
上記式(17-1)を上記式（１７）に代入すると、下記式（１８）が得られる。
【数２７】
2019191110.tif 000029
【００７２】
上記式（１３）を用いると、上記式（１８）の積分表示は下記式（１９）のようになる。
【数２８】
2019191110.tif 000030
【００７３】
上記式（１９）の右辺の１項目を上記式（１８）を用いて展開すると、下記式(19-1)が得られる。なお、以下では、表記の都合上、上記式（１９）の右辺の各項の式変形を記載する際は、時間積分を省略して記載する。
【数２９】
2019191110.tif 000031
【００７４】
上記式（１９）の右辺の２項目及び３項目を上記式（１８）を用いて展開すると、下記式(19-2)が得られる。
【数３０】
2019191110.tif 000032
【００７５】
上記式（１９）の右辺の４項目は、下記式(19-3)のように式変形できる。
【数３１】
2019191110.tif 000033
【００７６】
上記式(19-1)〜(19-3)を上記式（１９）に代入すると、下記式（２０）が得られる。
【数３２】
2019191110.tif 000034
【００７７】
仮想表面Ｓ＋を剛体３６の表面Ｓに収束させると、上記式（２０）は下記式（２１）のように書き直すことができる。
【数３３】
2019191110.tif 000035
【００７８】
ここで、高レイノルズ数流れでは、上記式（２１）の右辺の３項目は４項目に対して無視できるほど小さい。加えて、剛体３６が静止している場合は流速がないため、当該右辺の２項目はゼロに等しい。従って、上記式（２１）は、下記式（２２）のように簡略化できる。
【数３４】
2019191110.tif 000036
【００７９】
静止物体表面上では流速がないため、エンタルピーＢは下記式のように表される。
【数３５】
2019191110.tif 000037
【００８０】
上記式の右辺ｐ／ρをｐ／ρ＝ｐ＊と定義してエンタルピーＢを上記式（２２）に代入すると、下記式（２３）の関係式が得られる。この関係式が、「時間ドメインにおける表面圧力変動と流れ場との関係式」である。
【数３６】
2019191110.tif 000038
【００８１】
（ＰＤＳ計算式の導出）
上記式（２３）の関係式は、表面圧力変動ｐに基づく値（関数）ｐ＊（＝ｐ／ρ）をその両辺に含むため、観測点ｘにおける表面圧力変動ｐを流速ｕ及び渦度ωから直接的に計算することができない。このため、以下では、式変形により右辺からｐ＊を含む項を消去し、「観測点ｘにおける表面圧力変動ｐ」と「流れ場の流速ｕ及び渦度ω」との関係式を導出する。
【００８２】
図７に示すように、剛体３６の表面Ｓ上の観測点ｘの近傍領域を表面領域Ｓｘと定義し、表面領域Ｓｘの法線ベクトルをｎ（ｘ）とする。表面領域Ｓｘは十分に大きい平面である。加えて、法線ベクトルｎ（ｘ）上に、観測点ｘからの距離が等しい２個の点ｘ＋、ｘ−を定義する。点ｘ＋は剛体３６の外部に位置しており、点ｘ−は剛体３６の内部に位置している。なお、その他の文字については後述する。点ｘ＋、ｘ−をそれぞれ観測点として上記式（２２）に代入すると、下記式(22-1)、(22-2)が得られる。
【００８３】
【数３７】
2019191110.tif 000039
【００８４】
式(22-1)と式(22-2)の両辺をそれぞれ足し合わせてｘ−とｘ＋についてｘに近づくときの極限をとると、下記式（２４）が得られる。なお、以下の式変形では、Ｈ（ｘ＋）＝１及びＨ（ｘ−）＝０を利用している（式（１５）参照）。
【数３８】
2019191110.tif 000040
【００８５】
上記式（２４）の右辺の１項目及び２項目は、それぞれ下記式(24-1)及び(24-2)のように式変形できる。なお、式中のＳ＼Ｓｘは、剛体３６の表面Ｓのうち表面領域Ｓｘを含まない部分を表す（図７参照）。
【数３９】
2019191110.tif 000041
【数４０】
2019191110.tif 000042
【００８６】
上記式(24-1)、(24-2)を上記式（２４）に代入すると、下記式（２５）が得られる。
【数４１】
2019191110.tif 000043
【００８７】
フーリエ変換は下記式によって定義される。
【数４２】
2019191110.tif 000044
ここで、式中のＡθ（ｘ）は角周波数θにおける関数Ａ（ｘ，ｔ）の大きさ（振幅）を表す。Ａθ（ｘ）は複素数であり、その複素共役Ａθ＊（ｘ）は下記（２６）の関係式を満たす。
【数４３】
2019191110.tif 000045
【００８８】
上記式（２４）にフーリエ変換を作用し、畳み込み積分のフーリエ変換結果を用いると、下記式（２７）が得られる。
【数４４】
2019191110.tif 000046
【００８９】
ここで、グリーン関数Ｇをフーリエ変換したＧθ（ｘ，ｚ）は下記式（２８）のように表され、下記式（２９）の関係式（いわゆるヘルムホルツ(Helmholtz)方程式）を満たす。
【数４５】
2019191110.tif 000047
ここで、ｋ０は角周波数θ成分の波数（いわゆる音響波数）であり、ｋ０＝θ／ｃ０である。ｒは観測点ｘと評価点ｚとの距離である。
【００９０】
Ｇθ（ｘ，ｚ）の勾配は、上記式（２８）を用いて下記式（３０）のように表すことができる。なお、音響コンパクト条件（ｋ０ｒ＜＜１）においては、Ｇθ（ｘ，ｚ）の勾配は角周波数θに依存しておらず、距離ｒの関数のみ近似できる。加えて、式中のベクトルｒは、観測点ｘから評価点ｚへ向かうベクトルである。
【数４６】
2019191110.tif 000048
【００９１】
上記式（２７）の左辺である「観測点ｘにおける表面圧力変動ｐθ（ｘ）に基づく関数ｐ＊θ（ｘ）」の自己相関関数（上記式（１）参照）は、下記式（３１）のように表される。
【数４７】
2019191110.tif 000049
【００９２】
上記式（３１）は、Ｇθ（ｘ，ｚ）の勾配（式（３０）参照）の音響コンパクト条件を用いて下記式（３２）のように表される。
【数４８】
2019191110.tif 000050
【００９３】
ここで、着目角周波数帯[θｌ，θｈ]における任意の観測点ｘの物理量の相関距離ｌｘ,ｃｒ、相関面積Ｓｘ,ｃｒ（図７参照）及び相関体積Ｖｘ,ｃｒをそれぞれ下記式(33-1)、(33-2)及び(33-3)によって定義する。なお、式中のθｍは着目角周波数帯の中心角周波数である（θｍ＝２πｆｍ）。
【数４９】
2019191110.tif 000051
【００９４】
即ち、相関面積Ｓｘ,ｃｒは平均流速の二乗に比例しており、相関体積Ｖｘ,ｃｒは平均流速の三乗に比例している。
【００９５】
本実施形態の着目角周波数帯は比較的に大きいため、相関距離ｌｘ,ｃｒは比較的に小さい。従って、相関面積Ｓｘ,ｃｒ及び相関体積Ｖｘ,ｃｒも比較的に小さい。このため、図７に示すように、相関面積Ｓｘ,ｃｒは表面領域Ｓｘ内に含まれている（Ｓｘ,ｃｒ⊂Ｓｘ）。ここで、上記式（３２）の右辺の１項目の被積分項には剛体３６の表面Ｓ上の評価点ｙが含まれており、表面Ｓのうち表面領域Ｓｘを含まない領域（Ｓ＼Ｓｘ）で積分されている。このため、観測点ｘの物理量は評価点ｙの物理量と無相関となり、この結果、上記式（３２）の右辺の１項目はゼロに等しくなる。
【００９６】
加えて、以下では、任意の関数Ｆθとその複素共役関数Ｆθ＊との積を「関数Ｆθのノルム」と称し、下記式（３４）のように定義する。
【数５０】
2019191110.tif 000052
即ち、自己相関関数（式（１）参照）は、任意の関数Ｆθのノルムを着目角周波数帯において積分した値である。
【００９７】
これにより、上記式（３２）は、下記式（３５）のように表すことができる。
【数５１】
2019191110.tif 000053
【００９８】
上記式（３５）によると、着目角周波数帯における観測点ｘの表面圧力変動ｐθ（ｘ）は、フロードメイン３８内の流れ場の渦度ωと流速ｕの外積に依存している。本願発明者はこの点に着目し、上記式（３５）の右辺の空間積分の自己相関関数の値を低減すれば、観測点ｘにおける表面圧力変動ｐθ（ｘ）を低減できると考えた。この自己相関関数は、下記式（３６）に示すように式変形できる。
【数５２】
2019191110.tif 000054
【００９９】
ここで、式中のｚ１、ｚ２は、フロードメイン３８内の任意の評価点であり、ｚ１≠ｚ２である。着目角周波数帯における評価点ｚ１の物理量と評価点ｚ２の物理量との相関は、評価点ｚ１が評価点ｚ２の相関体積Ｖｚ２,ｃｒに含まれていない場合は無視することができる。このため、上記式（３６）は、更に下記式（３７）のように変形できる。
【数５３】
2019191110.tif 000055
【０１００】
上記式（３７）の右辺のうち、「着目角周波数帯について積分される部分」は、上記式（３０）を用いて下記式（３８）のように式変形できる。
【数５４】
2019191110.tif 000056
【０１０１】
上記式（３８）の流速ｕは、下記式(39-1)、(39-2)に示すように、平均成分としての平均流速Ｕ−と、乱れ成分としての乱れ流速ｕ〜に分解でき、渦度ωは、平均成分としての平均渦度Ω−と、乱れ成分としての乱れ渦度ω〜に分解できる。
【数５５】
2019191110.tif 000057
【０１０２】
加えて、  剛体モデル３６の空間座標系において、ベクトルｒを下記式（４０）によって規定する。
【数５６】
2019191110.tif 000058
【０１０３】
ここで、ｒ１、ｒ２及びｒ３は、観測点ｘを原点としたときのベクトルｒのｅ１成分、ｅ２成分及びｅ３成分である。
【０１０４】
これらの式(39-1)、(39-2)及び（４０）を上記式（３８）の右辺に代入すると、当該右辺は、下記式（４１）のように式変形できる。なお、乱れ成分同士の積は無視できるほど小さいため、下記式（４１）では、その記載を省略している。
【数５７】
2019191110.tif 000059
【０１０５】
ここで、２つの任意な関数Ｅθと関数Ｆθの相関関数を着目角周波数帯において下記式（４２）のように定義する。
【数５８】
2019191110.tif 000060
【０１０６】
上記式（１）及び式（４２）を用いると、渦一様等方性仮説より、下記式(43-1)〜(43-6)に示す関係式が得られる。
【数５９】
2019191110.tif 000061
ここで、ｋは流速ｕθと渦度ωθとの関係を規定した定数である。
【０１０７】
上記式(43-1)〜(43-6)を用いて上記式（４１）の右辺を整理すると、上記式（３８）の左辺は、最終的に下記（４４）のように表される。
【数６０】
2019191110.tif 000062
【０１０８】
上記式（４４）の右辺は、上記式（３７）の右辺のうちの「着目角周波数帯について積分される部分」に相当する。このため、式（４４）の右辺を式（３７）の該当部分に代入すると、下記式（４５）が得られる。
【数６１】
2019191110.tif 000063
【０１０９】
上記式（４５）の右辺の被積分項に「流れ場の密度の二乗ρ２」を乗算し、評価点ｚ２を一般の表記ｚに書き換えると、下記式（４６）に定義されるＰＤＳの計算式が得られる。本明細書では、このようにして導出された式（４６）を、「剛体３６の観測点ｘにおける表面圧力変動」に対する「フロードメイン３８内の評価点ｚにおける流れ場」の寄与度を表す指標として定義している。なお、ρ２を乗算するのは、ＰＤＳ計算式の導出過程でｐ＊＝ｐ／ρを用いたことによるものである。
【数６２】
2019191110.tif 000064
【０１１０】
なお、観測点ｘにおける表面圧力変動ｐθ（ｘ）とＰＤＳ（ｘ，ｚ）との間には、下記式（４７）に示す近似式が成立する（式（３５）、（４５）及び（４６）参照）。
【数６３】
2019191110.tif 000065
【０１１１】
即ち、「観測点ｘにおける表面圧力変動ｐθ（ｘ）のノルム」の着目角周波数帯における積分値（別言すれば、観測点ｘにおける表面圧力変動ｐθ（ｘ）の自己相関関数）は、「観測点ｘについてのＰＤＳをフロードメインで空間積分した値」として近似することができる。このため、ＰＤＳは、観測点ｘにおける表面圧力変動ｐθ（ｘ）に対する評価点ｚにおける流れ場の寄与度を表す指標として高い精度で機能する。以上が、ＰＤＳ計算式の導出についての説明である。
【０１１２】
＜風切音解析方法＞
続いて、走行中の車両モデル２０の風切音解析方法について具体的に説明する。計算装置１０の演算部１４のＣＰＵは、図８にフローチャートにより示したルーチンを実行するようになっている。ＣＰＵは、図８のステップ８００から処理を開始して以下に述べるステップ８０２〜ステップ８０８の処理を順に行う。
【０１１３】
ステップ８０２：ＣＰＵは、車両モデル２０を走行させる非定常ＣＦＤシミュレーションを開始し、解析範囲（解析表面及びフロードメイン２１内）において、圧力変動の時刻歴データｐ（ｘ，ｔ）、流速の時刻歴データｕ（ｚ，ｔ）及び渦度の時刻歴データω（ｚ，ｔ）を所定時間に亘ってそれぞれ算出する。時刻歴データｐ（ｘ，ｔ）は節点ｘ毎に算出され、時刻歴データｕ（ｚ，ｔ）及びω（ｚ，ｔ）は節点ｚ毎に算出される。
ステップ８０４：ＣＰＵは、ステップ８０２で算出された圧力変動の時刻歴データｐ（ｘ，ｔ）に基づいて、解析表面における圧力変動の振幅の時間推移を表すデータを作成し、当該データの表示指令を出力部１６に送信する。これにより、圧力変動の振幅の時間推移が表示画面１８にアニメーション表示される（図２Ｂ参照）。
【０１１４】
ステップ８０６：ＣＰＵは、ステップ８０２で算出された流速及び渦度の時刻歴データｕ（ｚ，ｔ）及びω（ｚ，ｔ）をそれぞれ平均化処理して、節点ｚ毎に、平均流速Ｕ−（ｚ）及び平均渦度Ω−（ｚ）を算出する。
ステップ８０８：ＣＰＵは、ステップ８０２で算出された圧力変動の時刻歴データｐ（ｘ，ｔ）を高速フーリエ変換処理して、着目角周波数帯における表面圧力変動ｐθ（ｘ）を算出する。表面圧力変動ｐθ（ｘ）は節点ｘ毎に算出される。加えて、ＣＰＵは、ステップ８０２で算出された流速の時刻歴データｕ（ｚ，ｔ）を高速フーリエ変換処理して、着目角周波数帯における乱れ流速ｕ〜θ（ｚ）の自己相関関数を算出する。この自己相関関数は節点ｚ毎に算出される。ステップ８０２〜ステップ８０８までの処理が非定常ＣＦＤ計算処理に相当する。
【０１１５】
ステップ８１０：ＣＰＵは、ステップ８０８で算出された表面圧力変動ｐθ（ｘ）に基づいて、「着目角周波数帯の中心角周波数θｍにおける表面圧力変動ｐθ（ｘ）の分布図を表すデータ」を作成し、当該データの表示指令を出力部１６に送信する。これにより、表面圧力変動ｐθ（ｘ）の分布図が表示画面１８に画像表示される（図３Ａ参照）。その後、ＣＰＵは、ステップ８１２に進む。
【０１１６】
ステップ８１２では、ＣＰＵは、作業者により観測点ｘが選択されたか否かを判定する。この観測点ｘは、表面圧力変動ｐθ（ｘ）の分布図（ステップ８１０参照）に基づいて作業者により選択される。観測点ｘが選択されたと判定した場合、ＣＰＵは、ステップ８１２にて「Ｙｅｓ」と判定し、後述するステップ８１４の処理を行う。これに対し、観測点ｘが選択されていないと判定した場合、ＣＰＵは、ステップ８１２にて「Ｎｏ」と判定し、再度ステップ８１２の判定を行う。ＣＰＵは、観測点ｘが選択されたと判定するまでこの処理を所定の演算時間が経過する毎に繰り返す。
【０１１７】
ステップ８１４：ＣＰＵは、演算部１４のＲＡＭに格納されているＰＤＳ計算式（式（４６）参照）を用いて、ステップ８１２で選択された観測点ｘについてのＰＤＳ（ｘ，ｚ）を節点ｚ毎に算出する。この計算式には、ステップ８０６で算出された平均流速Ｕ−（ｚ）及び平均渦度Ω−（ｚ）並びにステップ８０８で算出された乱れ流速ｕ〜θ（ｚ）の自己相関関数等が代入される。その後、ＣＰＵは、以下のステップ８１６の処理を行う。
【０１１８】
ステップ８１６：ＣＰＵは、ＰＤＳの算出を終了したことを示すメッセージの表示指令を出力部１６に送信する。これにより、当該メッセージが表示画面１８に表示されるため、作業者は算出処理が終了したことを認識できる。その後、ＣＰＵは、ステップ８１８に進む。
【０１１９】
ステップ８１８では、ＣＰＵは、作業者によりＰＤＳの参照値が入力されたか否かを判定する。参照値が入力されたと判定した場合、ＣＰＵは、ステップ８１８にて「Ｙｅｓ」と判定し、以下のステップ８２０の処理を行う。
【０１２０】
ステップ８２０：ＣＰＵは、ステップ８１８で入力された参照値を有するＰＤＳの等値面を表すデータを作成し、当該データの表示指令を出力部１６に送信する。これにより、作業者が入力した参照値を有するＰＤＳの等値面が表示画面１８に画像表示される（図４Ａ〜図４Ｃ参照）。その後、ＣＰＵは、ステップ８２２に進む。
【０１２１】
ステップ８２２では、ＣＰＵは、作業者により着目評価点ｚが選択されたか否かを判定する。この着目評価点ｚは、ＰＤＳの等値面（ステップ８２０参照）に基づいて作業者により選択される。着目評価点ｚが選択されたと判定した場合、ＣＰＵは、ステップ８２２にて「Ｙｅｓ」と判定し、後述するステップ８２４（図９）の処理を行う。
【０１２２】
これに対し、着目評価点ｚが選択されていないと判定した場合（具体的には、作業者が現時点におけるＰＤＳ等値面からは着目評価点ｚを選択できない場合、又は、作業者が着目評価点ｚを選択中の場合）、ＣＰＵは、ステップ８２２にて「Ｎｏ」と判定し、再度ステップ８１８の判定を行う。
【０１２３】
そして、ステップ８１８にて参照値が入力されたと判定した場合（即ち、作業者によりＰＤＳの参照値が再入力された場合）、ＣＰＵは、ステップ８１８にて「Ｙｅｓ」と判定してステップ８２０の処理を行い、ステップ８２２にて着目評価点ｚが選択されたか否かを再度判定する。一方、ステップ８１８にて参照値が入力されていないと判定した場合、ＣＰＵは、ステップ８１８にて「Ｎｏ」と判定してステップ８２２に進み、ステップ８２２にて着目評価点ｚが選択されたか否かを再度判定する。ＣＰＵは、着目評価点ｚが選択されたと判定するまでステップ８１８〜ステップ８２２の処理を繰り返す。
【０１２４】
これに対し、ステップ８１６から直接ステップ８１８に進んだ場合において参照値が入力されていないと判定した場合（即ち、初回の参照値が入力されていない場合）、ＣＰＵは、ステップ８１８にて「Ｎｏ」と判定してステップ８２２に進む。上述したように、着目評価点ｚはＰＤＳの等値面に基づいて選択される（即ち、ステップ８２０を経て選択される）ため、初回の参照値が入力されていない場合（即ち、ステップ８２０を経ていない場合）、ＣＰＵはステップ８２２にて「Ｎｏ」と判定する。
【０１２５】
ステップ８２４（図９）：ＣＰＵは、ステップ８１４で算出されたＰＤＳのうち、ステップ８２２で選択された着目評価点ｚにおけるＰＤＳ（ｘ，ｚ）を構成するパラメータの数値データの表示指令を出力部１６に送信する。これらのパラメータは、具体的には、着目評価点ｚにおける平均流速Ｕ−（ｚ）、乱れ流速ｕ〜θ（ｚ）、及び、平均渦度Ω−（ｚ）である（ステップ８０６、８０８参照）。これにより、ＰＤＳ（ｘ，ｚ）を構成するパラメータの数値データが表示画面１８に表示される。その後、ＣＰＵは、ステップ８２６に進む。
【０１２６】
ステップ８２６では、ＣＰＵは、原因パラメータが特定されたか否かを判定する。この原因パラメータは、ＰＤＳ（ｘ，ｚ）を構成するパラメータの数値データ（ステップ８２４参照）に基づいて、作業者により特定される。原因パラメータが特定されたと判定した場合、ＣＰＵは、ステップ８２６にて「Ｙｅｓ」と判定し、後述するステップ８２８の処理を行う。これに対し、原因パラメータが特定されていないと判定した場合、ＣＰＵは、ステップ８２６にて「Ｎｏ」と判定し、再度ステップ８２６の判定を行う。ＣＰＵは、原因パラメータが特定されたと判定するまでこの処理を所定の演算時間が経過する毎に繰り返す。
【０１２７】
ステップ８２８：ＣＰＵは、ステップ８２２で選択された着目評価点ｚを通る所定の平面（典型的には、ｅ１ｅ２平面に平行な平面）上における原因パラメータの分布図を表すデータを作成する。ＣＰＵは、当該データの表示指令を出力部１６に送信する。これにより、原因パラメータの分布図が表示画面１８に表示される（図５Ｂ参照）。その後、ＣＰＵは、ステップ８３０に進み、本ルーチンを終了する。
【０１２８】
＜本実施装置による表面圧力変動（風切音）の改善効果＞
次に、本実施装置を用いて車両モデル２０の解析範囲の形状を変更したときの表面圧力変動ｐθ（ｘ）の改善効果について図１０Ａ〜図１１Ｂを参照して具体的に説明する。図５Ａ及び図５Ｂの例では、空気の流れが領域３４で剥離する可能性を低減できるような右側サイドミラー２２の形状が検討される。図１０Ａは、車両モデル２０の右側サイドミラー２２の概形を示す。図１０Ａに示すように、右側サイドミラー２２は、ミラーが設けられているミラー部２２ａと、車体に固定されている基部２２ｂと、ミラー部２２ａと基部２２ｂとを接続する接続部２２ｃと、を有する。領域３４は、接続部２２ｃの近傍に位置している。接続部２２ｃは比較的に肉薄であり、接続部２２ｃとミラー部２２ａとの境界には段差が形成されている。本願発明者は、この段差により空気の流れが領域３４で剥離していると考え、段差を緩和できる形状を検討した。図１０Ｂは、上記検討に基づいて形状が変更された車両モデルの右側サイドミラー１２２を示す。図１０Ｂに示すように、右側サイドミラー１２２の接続部１２２ｃは、接続部２２ｃよりも肉厚となっている。これにより、接続部１２２ｃとミラー部１２２ａとの境界がなだらかになっており段差が緩和されている。
【０１２９】
本願発明者は、上記の形状変更による表面圧力変動ｐθ（ｘ）の改善効果を検証するため、各モデルの解析表面における表面圧力変動ｐθ（ｘ）の分布図を比較した。図１１Ａは、車両モデル２０の分布図を示し、図１１Ｂは、形状変更後の車両モデルの分布図を示す。図１１Ａと図１１Ｂとを比較すると、右側サイドミラー２２の形状を変更することにより右側サイドミラー１２２近傍の表面圧力変動ｐθ（ｘ）が大幅に低減している。以上より、本実施装置によれば、ＰＤＳを利用して車両モデルの表面圧力変動に大きく寄与している流れ場の位置を特定することにより、車両モデルの表面圧力変動（即ち、風切音）を効率的に且つ確実に低減できることが示された。
【０１３０】
＜指標としてのＰＤＳの信頼性及び汎用性＞
上述したように、「節点ｘについてのＰＤＳをフロードメインで空間積分した値」は、「当該節点ｘにおける表面圧力変動ｐθ（ｘ）の自己相関関数」の近似値である（式（４７）参照）。以下では、この近似精度を検証してＰＤＳの指標としての信頼性について検討する。具体的には、「ＰＤＳをフロードメインで空間積分した値」に基づいて算出（予測）される表面圧力変動（以下、「予測表面圧力変動」とも称する。）を、非定常ＣＦＤ計算用のソフトウエアによって算出される表面圧力変動ｐθと比較することにより、予測表面圧力変動の近似精度を検証する。なお、本実施形態では上記ソフトウエアとしてＰｏｗｅｒＦＬＯＷが使用されるが、ソフトウエアの種類はこれに限られない。
【０１３１】
加えて、以下では、車両モデルだけではなく、車両モデルとは異なる形状を有するフォアステップモデルについても上記近似精度を検証することにより、ＰＤＳの指標としての汎用性についても検討する。
【０１３２】
図１２Ａは、上記の説明で用いられた車両モデル２０である。図１２Ｂは、フォアステップモデル４０の模式図である。フォアステップモデル４０は、基準面４２と、基準面４２上に配置された直方体形状を有するステップ４４と、基準面４２上にステップ４４から離間して配置されたノズル４６と、を有する。ステップ４４の上面４４ａは基準面４２と平行である。ステップ４４の側面４４ｂ（ノズル４６側の面）と基準面４２との間には段差が形成されている。ノズル４６は、ステップ４４の側面４４ｂと平行に対向している送風口４６ａを有する。送風口４６ａは、ステップ４４に向かって風を送る。送風方向と直交する方向を「幅方向」と規定すると、幅方向におけるステップ４４の長さａは、幅方向におけるノズル４６の長さｂよりも長い。ステップ４４の高さはノズル４６の高さよりも低い。ノズル４６の幅方向における中心は、ステップ４４の幅方向における中心と、幅方向において一致している。本検証では、車両モデル２０の右側フロントサイドガラス２５とフォアステップモデル４０のステップ４４の上面４４ａが解析対象とされる。
【０１３３】
図１３Ａ及び図１３Ｂは、それぞれ、車両モデル２０の右側フロントサイドガラス２５における予測表面圧力変動及び表面圧力変動ｐθの分布図である。着目周波数帯は１ｋＨｚに設定されている。図１４Ａ及び図１４Ｂは、それぞれ、フォアステップモデル４０のステップ４４の上面４４ａにおける予測表面圧力変動及び表面圧力変動ｐθの分布図である。着目周波数帯は５００Ｈｚに設定されている。予測表面圧力変動及び表面圧力変動ｐθが大きいほど濃色で表示されている。図１３Ａと図１３Ｂを比較すると、予測表面圧力変動の大きさ及び位置は、表面圧力変動ｐθの大きさ及び位置と精度良く一致している。同様に、図１４Ａと図１４Ｂを比較すると、予測表面圧力変動の大きさ及び位置は、表面圧力変動ｐθの大きさ及び位置と精度良く一致している。以上より、ＰＤＳに基づいて算出される表面圧力変動（予測表面圧力変動）は、ソフトウエアのＣＦＤ計算により算出される表面圧力変動に対して高い近似精度を有することが示された。これは、ＰＤＳが指標として高い信頼性を有することを示している。
【０１３４】
図１５は、「表面圧力変動ｐθに対する予測表面圧力変動の誤差」を車両モデル２０及びフォアステップモデル４０についてそれぞれ示したグラフである。車両モデル２０についての誤差は、図１３Ａ及び図１３Ｂに示す２つの分布図の対応する各節点ｘにおける「表面圧力変動ｐθに対する予測表面圧力変動の差」の平均を圧力変動レベル（dB）に変換した値である。フォアステップモデル４０についての誤差は、図１４Ａ及び図１４Ｂに示す２つの分布図の対応する各節点ｘにおける「表面圧力変動ｐθに対する予測表面圧力変動の差」の平均を圧力変動レベル（dB）に変換した値である。車両モデル２０の誤差とフォアステップモデル４０の誤差を比較すると、何れも３ｄＢ強でありほぼ一致している。以上より、表面圧力変動ｐθに対する予測表面圧力変動の誤差は、解析対象モデルの形状に関係なくほぼ一定であることが示された。これは、ＰＤＳが指標として高い汎用性を有することを示している。
【０１３５】
本実施装置の作用効果について説明する。本実施装置では、車両モデルを走行させる非定常ＣＦＤシミュレーションにより算出された物理量（平均流速Ｕ−、平均渦度Ω−及び乱れ流速ｕ〜θ）に基づいてＰＤＳを算出する。このＰＤＳは、表面圧力変動ｐθ（ｘ）に対する流れ場の寄与度を示す指標である。このため、ＰＤＳを指標として用いることにより、表面圧力変動ｐθ（ｘ）に対する流れ場の寄与度を評価点ｚ毎に算出できる。従って、フロードメイン２１内のどの位置における流れ場が観測点ｘにおける車両モデル２０の表面圧力変動ｐθ（ｘ）に大きく寄与しているかを適切に特定することが可能となる。
【０１３６】
加えて、本実施装置では、観測点ｘについてのＰＤＳをフロードメイン２１で空間積分した値は、当該観測点ｘにおける表面圧力変動ｐθ（ｘ）の自己相関関数の近似値である。このため、ＰＤＳの挙動は、表面圧力変動ｐθ（ｘ）の挙動と良く一致する。即ち、ＰＤＳは、表面圧力変動ｐθ（ｘ）に対する流れ場の寄与度を表す指標として高い信頼性を有する（図１２Ａ〜図１５参照）。このため、車両モデルの表面圧力変動ｐθ（ｘ）に大きく寄与している流れ場の位置を適切に特定することができる。
【０１３７】
具体的には、従来は、上記関係を定量的に評価することができなかったため、技術者によって流れ場の解釈が異なることがあった。加えて、解析は、流れ場の断面における物理量の分布に基づいて行われるため、三次元空間における流れ場の情報が欠落し、表面圧力変動が大きい原因を見逃してしまう可能性があった。これに対し、本実施装置によれば、ＰＤＳを用いることにより表面圧力変動ｐθ（ｘ）と流れ場との関係を定量的に評価できるため、技術者間で流れ場の解釈が異なる可能性を大幅に低減でき、風切音の解析精度を向上することができる。加えて、本実施装置では、表面圧力変動への寄与度が高い流れ場の位置を特定してから当該位置に断面における物理量の分布を評価するため、流れ場の情報が欠落する可能性を大幅に低減でき、結果として、表面圧力変動が大きい原因を確実に発見することができる。
【０１３８】
更に、本実施装置では、表示画面１８にＰＤＳの等値面が表示されるため、作業者は、参照値を適宜設定することにより、所望のＰＤＳ値を有する等値面を視認することができる。この結果、所望のＰＤＳ値に対応している評価点ｚを効率的に選択することができ、車両モデル２０の表面圧力変動ｐθ（ｘ）に大きく寄与している流れ場の位置を効率的に特定することができる。
【０１３９】
更に、本実施装置では、観測点ｘは、非定常ＣＦＤ計算により算出された表面圧力変動の大きさに基づいて選択される。表面圧力変動は風切音と強い相関を有している。このため、この構成によれば、風切音をより適切に解析することができる。
【０１４０】
（変形例１）
次に、変形例１に係る風切音解析装置（以下、「第１変形装置」とも称する。）について説明する。第１変形装置は、原因パラメータの特定が作業者ではなく装置により行われる点で上記実施装置と異なっている。具体的には、第１変形装置は、図９のステップ８２４及び８２６の処理に代えて、以下の処理を実行する。即ち、第１変形装置の計算装置の演算部は、ステップ８２２で選択された着目評価点ｚにおけるＰＤＳを構成する３つのパラメータ（平均流速Ｕ−、平均渦度Ω−及び乱れ流速ｕ〜θ）の中からＰＤＳ値への寄与が相対的に大きいパラメータを原因パラメータとして特定する。具体的には、演算部は、着目評価点ｚにおける３つのパラメータの数値を、「着目評価点ｚにおけるＰＤＳ値よりも低いＰＤＳ値に対応している評価点ｚ」における３つのパラメータの数値と比較することにより原因パラメータを特定する。この構成によれば、原因パラメータを作業者が特定する必要がなくなるため、ＰＤＳが高い原因となっているパラメータを容易に把握でき、車両モデルの形状をより効率的に検討、変更することができる。
【０１４１】
（変形例２）
続いて、変形例２に係る風切音解析装置（以下、「第２変形装置」とも称する。）について説明する。第２変形装置は、着目評価点ｚの選択（抽出）が作業者ではなく装置により行われる点で上記実施装置と異なっている。具体的には、第２変形装置は、図８のステップ８１８〜ステップ８２２の処理に代えて、以下の処理を実行する。即ち、第２変形装置の計算装置の演算部は、ステップ８１６の処理が終了すると、任意の値を有するＰＤＳの等値面データを作成し、当該データにより示される等値面の個数が２個になるまでＰＤＳの値を変更しながら等値面データを作成する（厳密には、２個になってからも各等値面の表面積ができるだけ小さくなるまでＰＤＳの値を増加する。）。そして、観測点ｘまでの距離が長い方の等値面上の任意の点を着目評価点ｚとして抽出する。この構成によれば、着目評価点ｚを作業者が選択する必要がなくなるため、車両モデルの表面圧力変動に大きく寄与している流れ場の位置をより効率的に特定することができる。なお、第２変形装置では最大（ほぼ最大）のＰＤＳ値に対応している評価点ｚが着目評価点ｚとして抽出されるが、この構成に限られず、予め設定された条件を満たす着目評価点ｚが抽出される構成であっても良い。
【０１４２】
以上、実施形態及び変形例の風切音解析装置及び風切音解析方法について説明したが、本発明は上記実施形態及び変形例に限定されるものではなく、本発明の目的を逸脱しない限りにおいて種々の変更が可能である。
【０１４３】
例えば、上記実施形態及び変形例では解析対象となる構造物モデルとして車両モデルが用いられたが、この構成に限られず、例えば、航空機及び船舶等の構造物モデル等が用いられてもよい。
【０１４４】
加えて、上記実施形態及び変形例では、表面圧力変動が比較的に大きい観測点が解析対象とされたが、この構成に限られず、例えば、所定範囲における表面圧力変動の平均値を有する節点ｘを観測点として解析してもよい。
【０１４５】
更に、上記実施形態及び変形例では、等値面上の節点ｚが着目評価点ｚとして選択可能であったが、この構成に限られず、例えば、等値面内部の節点ｚが着目評価点ｚとして選択可能であってもよい。
【０１４６】
更に、上記実施形態及び変形例では、観測点ｘは作業者により選択されたが、この構成に限られず、例えば、表面圧力変動の値に基づいて風切音解析装置により選択されてもよい。
【０１４７】
更に、上記実施形態及び変形例では、観測点ｘは表面圧力変動の分布図に基づいて選択されたが、この構成に限られず、例えば、表面圧力変動以外の物理量に基づいて選択されてもよいし、任意の節点ｘが作業者により観測点ｘとして選択されてもよい。
【符号の説明】
【０１４８】
１０：計算装置、１２：入力部、１４：演算部、１６：出力部、１８：表示画面、２０：車両モデル、２１：フロードメイン、２２：右側サイドミラー、２４：右側フロントピラー、２５：右側フロントサイドガラス、２６：領域、２８，３０，３２（３２ａ，３２ｂ）：等値面、３４：領域、３６：剛体モデル、３８：フロードメイン、４０：フォアステップモデル、４２：基準面、４４：ステップ、４６：ノズル
【図１】
2019191110.tif 000067
【図２Ａ】
2019191110.tif 000068
【図２Ｂ】
2019191110.tif 000069
【図３Ａ】
2019191110.tif 000070
【図３Ｂ】
2019191110.tif 000071
【図４Ａ】
2019191110.tif 000072
【図４Ｂ】
2019191110.tif 000073
【図４Ｃ】
2019191110.tif 000074
【図５Ａ】
2019191110.tif 000075
【図５Ｂ】
2019191110.tif 000076
【図６】
2019191110.tif 000077
【図７】
2019191110.tif 000078
【図８】
2019191110.tif 000079
【図９】
2019191110.tif 000080
【図１０Ａ】
2019191110.tif 000081
【図１０Ｂ】
2019191110.tif 000082
【図１１Ａ】
2019191110.tif 000083
【図１１Ｂ】
2019191110.tif 000084
【図１２Ａ】
2019191110.tif 000085
【図１２Ｂ】
2019191110.tif 000086
【図１３Ａ】
2019191110.tif 000087
【図１３Ｂ】
2019191110.tif 000088
【図１４Ａ】
2019191110.tif 000089
【図１４Ｂ】
2019191110.tif 000090
【図１５】
2019191110.tif 000091