特許 7475103 シミュレータの開発方法、情報処理システム及びプログラム

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B1)

(11)【特許番号】

(24)【登録日】2024-04-18

(45)【発行日】2024-04-26

(54)【発明の名称】シミュレータの開発方法、情報処理システム及びプログラム

(51)【国際特許分類】

   G06F 30/28 20200101AFI20240419BHJP

   G06F 113/08 20200101ALN20240419BHJP

【ＦＩ】

G06F30/28

G06F113:08

【請求項の数】 9

(21)【出願番号】P 2024018919

(22)【出願日】2024-02-09

【審査請求日】2024-02-09

【早期審査対象出願】

(73)【特許権者】

【識別番号】304021417

【氏名又は名称】国立大学法人東京工業大学

(74)【代理人】

【識別番号】100103894

【弁理士】

【氏名又は名称】家入 健

(72)【発明者】

【氏名】青木 尊之

(72)【発明者】

【氏名】渡辺 勢也

(72)【発明者】

【氏名】イン イクイ

【審査官】堀井 啓明

(56)【参考文献】

【文献】特開２０１６－１３８１（ＪＰ，Ａ）

【文献】特開２０２１－８８９７４（ＪＰ，Ａ）

【文献】中国特許出願公開第１０２３２３９６３（ＣＮ，Ａ）

(58)【調査した分野】(Int.Cl.，ＤＢ名)

Ｇ０６Ｆ３０／００－３０／２８

(57)【特許請求の範囲】

【請求項1】

飛行体の胴体周辺の空間を、計算格子を用いて離散化し、離散化された流体方程式を解くことで空力計算の数値計算手法を適用するステップと、
前記飛行体のブレードに対してアクチュエータラインモデルを適用するステップと、
前記数値計算手法の適用結果と、前記ブレードに対する前記アクチュエータラインモデルの適用結果と、を用いることで前記飛行体の航行をシミュレートするステップと、
をコンピュータが実行するシミュレータの開発方法。

【請求項2】

前記コンピュータは、前記ブレードを構成するマーカー粒子の間隔、又はガウシアンフィルタの影響半径のパラメータ、の少なくともいずれかの情報の入力を受け付け、前記アクチュエータラインモデルの適用において受け付けた前記情報を使用する、
請求項１に記載のシミュレータの開発方法。

【請求項3】

前記コンピュータは、前記ブレードの直径の１０００分の１以上のサイズの計算格子を、前記ブレード周辺の計算領域において、前記アクチュエータラインモデルの適用に用いる、
請求項１又は２に記載のシミュレータの開発方法。

【請求項4】

飛行体の胴体周辺の空間を、計算格子を用いて離散化し、離散化された流体方程式を解くことで空力計算の数値計算手法を適用する第１適用部と、
前記飛行体のブレードに対してアクチュエータラインモデルを適用する第２適用部と、
前記数値計算手法の適用結果と、前記ブレードに対する前記アクチュエータラインモデルの適用結果と、を用いることで前記飛行体の航行をシミュレートするシミュレータ部と、を備える
情報処理システム。

【請求項5】

前記ブレードを構成するマーカー粒子の間隔、又はガウシアンフィルタの影響半径のパラメータ、の少なくともいずれかの情報の入力を受け付ける入力部を更に備え、
前記シミュレータ部は、前記アクチュエータラインモデルの適用において受け付けた前記情報を使用する、
請求項４に記載の情報処理システム。

【請求項6】

前記シミュレータ部は、前記ブレードの直径の１０００分の１以上のサイズの計算格子を、前記ブレード周辺の計算領域において、前記アクチュエータラインモデルの適用に用いる、
請求項４又は５に記載の情報処理システム。

【請求項7】

飛行体の胴体周辺の空間を、計算格子を用いて離散化し、離散化された流体方程式を解くことで空力計算の数値計算手法を適用するステップと、
前記飛行体のブレードに対してアクチュエータラインモデルを適用するステップと、
前記数値計算手法の適用結果と、前記ブレードに対する前記アクチュエータラインモデルの適用結果と、を用いることで前記飛行体の航行をシミュレートするステップと、
をコンピュータに実行させるプログラム。

【請求項8】

前記コンピュータに、前記ブレードを構成するマーカー粒子の間隔、又はガウシアンフィルタの影響半径のパラメータ、の少なくともいずれかの情報の入力を受け付け、前記アクチュエータラインモデルの適用において受け付けた前記情報を使用させる、
請求項７に記載のプログラム。

【請求項9】

前記コンピュータに、前記ブレードの直径の１０００分の１以上のサイズの計算格子を、前記ブレード周辺の計算領域において、前記アクチュエータラインモデルの適用に使用させる、
請求項７又は８に記載のプログラム。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、シミュレータの開発方法、情報処理システム及びプログラムに関する。

【背景技術】

【0002】

飛行機は、その開発段階において、安全な飛行が可能かどうかを検証する必要がある。特に、近年実用化されつつある空飛ぶクルマ等の新しいタイプの有人飛行体については、このような飛行体が新しい技術であり、かつ搭乗する人間の安全性を確保する必要があるため、より慎重な安全性の検証が必要となる。

【0003】

実験では設定が非常に困難な航行条件があったり、検証すべき航行条件の数が膨大であったりするため、実機又は模型を用いて現実の検証実験を行うことには限界がある。また、実験では測定できない内容もあるため、コンピュータ上の数値シミュレーションを用いた検証がなされることが多い。シミュレーションには、例えば、数値流体力学の手法が用いられる。

【0004】

特に、空飛ぶクルマは、無人機と比較して、非常に高い航行安全性が求められる。例えば、複数のマルチコプターが接近（特に上下方向において接近）する場合、マルチコプターが建物等に接近する場合、いずれかのローターが異常な回転（例えば停止）をする場合、予期せぬ強風が発生する場合等の状況が想定される。シミュレーションを実行することで、このような状況において、マルチコプターがどのようにして姿勢を立て直し安全な航行を維持するかを発見することが可能になる。また、シミュレーションを実行することで、安全性を向上させる新しいローターの配置（又は追加）等の検証、及び検証結果の設計への反映も可能となる。

【0005】

関連技術として、非特許文献１には、都市交通システムに用いられるクワッドローターのエアタクシーに対して、数値流体力学シミュレーションを適用する技術が開示されている。

【先行技術文献】

【非特許文献】

【0006】

【文献】Patricia Ventura Diaz and Seokkwan Yoon, “High-Fidelity Simulations of a Quadrotor Vehicle for Urban Air Mobility”, AIAA SciTech Forum 2022, January 3-7, 2022, San Diego, CA & Virtual

【発明の概要】

【発明が解決しようとする課題】

【0007】

マルチコプター等の飛行体では、飛行時に、ローターを構成する、薄く細いローターブレード（回転翼；以下、単にブレードとも記載する）が高速に回転する。このような飛行体に対して、シミュレーションを行うために計算格子を用いて空間を離散化し、流体方程式を直接数値計算する方法を適用することが考えられる。この場合、ブレード表面近傍には非常に細かい計算格子を配置する必要があることから、計算に用いる総計算格子点数が膨大となる。そのため、計算の各時間ステップにおいてブレードにかかる力（すなわち推力）および周囲の空気の流れの変化を求めるために膨大な計算量が必要になる。

【0008】

また、計算の１時間ステップは、非常に短い時間間隔となる。したがって、コンピュータは、機体の並進・回転・姿勢の変化といった航行の振舞いを確認できる程度の期間（例えば１～２分）までシミュレーションの計算をすることが極めて困難になる。例えば、スーパーコンピュータを使って上記の数値シミュレーションを実行した場合でも、せいぜいブレードが３０回転する時間（例えば２～３秒）までしか、シミュレーションする時間を進めることができないという課題があった。非特許文献１に記載の数値流体力学シミュレーションはブレードの回転により得られる推力を計算する目的でなされており、機体の航行を計算するものではなかった。

【0009】

本発明は、このような問題点に鑑みてなされたものであり、必要な計算量を削減することに寄与するシミュレータの開発方法、情報処理システム及びプログラムを提供するものである。

【課題を解決するための手段】

【0010】

本発明の一態様に係るシミュレータの開発方法は、
飛行体の胴体周辺の空間を、計算格子を用いて離散化し、離散化された流体方程式を解くことで空力計算の数値計算手法を適用するステップと、
前記飛行体のブレードに対してアクチュエータラインモデルを適用するステップと、
前記数値計算手法の適用結果と、前記ブレードに対する前記アクチュエータラインモデルの適用結果と、を用いることで前記飛行体の航行をシミュレートするステップと、
をコンピュータが実行するものである。

【0011】

本発明の一態様に係る情報処理システムは、
飛行体の胴体周辺の空間を、計算格子を用いて離散化し、離散化された流体方程式を解くことで空力計算の数値計算手法を適用する第１適用部と、
前記飛行体のブレードに対してアクチュエータラインモデルを適用する第２適用部と、
前記数値計算手法の適用結果と、前記ブレードに対する前記アクチュエータラインモデルの適用結果と、を用いることで前記飛行体の航行をシミュレートするシミュレータ部と、を備える。

【0012】

本発明の一態様に係るプログラムは、
飛行体の胴体周辺の空間を、計算格子を用いて離散化し、離散化された流体方程式を解くことで空力計算の数値計算手法を適用するステップと、
前記飛行体のブレードに対してアクチュエータラインモデルを適用するステップと、
前記数値計算手法の適用結果と、前記ブレードに対する前記アクチュエータラインモデルの適用結果と、を用いることで前記飛行体の航行をシミュレートするステップと、
をコンピュータに実行させる。

【発明の効果】

【0013】

本発明によれば、必要な計算量を削減することに寄与するシミュレータの開発方法、情報処理システム及びプログラムを提供することができる。

【図面の簡単な説明】

【0014】

【図1】図１は、情報処理システムの一例を示すブロック図である。

【図2】図２は、マルチコプターの一例を示す。

【図3】図３は、アクチュエータラインモデルを適用した場合のブレードの一例を示す。

【図4】図４は、ブレードに配置されたマーカー粒子が周囲の空気から受ける力の一例を示す。

【図5】図５は、情報処理システムの代表的な処理の一例を示すフローチャートである。

【図6】図６は、流れの速度の大きさを示す数値計算のデータである。

【図7】図７は、推力について、実験結果と、計算結果とを比較したグラフである。

【図8】図８は、アクチュエータラインモデルを用いて計算したブレードの推力の反作用が周囲の空気にダウンウォッシュとして与えられる状態を示す計算結果である。

【図9】図９は、実施の形態に係る情報処理装置のハードウェア構成例を示すブロック図である。

【発明を実施するための形態】

【0015】

以下、図面を参照して、本発明の実施の形態について説明する。なお、発明を実施するための形態における以下の記載及び図面は、説明の明確化のため、適宜、省略及び簡略化がなされている。また、本実施の形態では、明記のない限り、複数の項目について「複数の項目の少なくともいずれか」が定義された場合、その定義は、任意の１つの項目を意味しても良いし、全ての項目を含む、任意の複数の項目を意味しても良い。

【0016】

例示的な実施形態を説明するためにいずれか１つの図に示された特徴またはステップのすべてが必ずしも必須ではなく、任意の一部の特徴またはステップが省略されてもよい。また、明細書及び図面に記載されたステップの順序は、適宜変更されてもよい。

【0017】

［関連技術］
まず、複数のローターが取り付けられた空飛ぶクルマ等のマルチコプターに関する自由航行シミュレータの関連技術について説明する。マルチコプターは、重力下においてブレードが回転することにより、周囲の空気から推力を得て、自由に航行することが可能になる。自由航行シミュレータは、事前に与えられたマルチコプターの形状および重量などの情報に対して、時間的に変化するブレードの回転数だけを変数として与え、マルチコプターが空中でどのように航行するかをシミュレーションで予測するプログラムである。

【0018】

自由航行シミュレータでは、ブレードを含めた機体が周囲の空気から受ける力の計算と、物体（すなわち機体）が移動することで周囲の空気に対して与える影響の計算がなされる。この力の計算において、ブレード部分については、ブレードが回転することによって生じる力が考慮され、機体全体については、機体が飛行することによって受ける力が考慮される。

【0019】

コンピュータが、機体が周囲の空気から受ける力を精度高く計算するためには、粘性項付の流体方程式（Navier-Stokes方程式）を解く必要がある。高速で回転するブレードが受ける力をコンピュータが計算するための方法は、物体表面の境界層を解像できる高解像度の計算格子を用いるか、壁関数モデルを用いるか、ブレードを含んだ近傍の計算にＤＥＳ（Detached Eddy Simulation）を用いるか、等のいくつかの方法がある。しかしながら、いずれの方法であっても、コンピュータが実行する計算量は膨大なものとなる。上記の方法のうち、計算の負荷が最も軽いと考えられる壁関数を用いる方法においても、コンピュータは、計算において、常に高解像度の計算格子をブレード近傍に割り当てる必要がある。したがって、計算の負荷は非常に大きくなる。

【0020】

また、必要な計算量を軽減するため、コンピュータが粘性項のない流体方程式であるEuler方程式を用いて計算する方法も考えられる。この場合、空気流に関する境界層が生成されないため、計算格子の解像度に対する制約は物体形状のみになる。しかしながら、この方法を用いた場合、コンピュータは、ブレードの回転と垂直方向の揚力しか計算することができず、ブレードから生成されるダウンウォッシュ等の空気流が全く計算できない。そのため、この方法は、自由航行シミュレータには適用できない。

【0021】

本発明の目的の一例は、シミュレーション対象となる飛行体にブレードの回転条件だけを変数として与え、コンピュータの仮想空間の中で飛行体を自由航行させることのできるシミュレータを開発することである。そのためのアクチュエータラインモデルを導入（採用）することが、本発明の新規な特徴である。

【0022】

なお、本開示における「飛行体」は、ブレード部分と胴体部分を有する任意の飛行する物体であり、例えば、有人又は無人操縦のドローン、空飛ぶクルマ、飛行機、宇宙船などを含む。「飛行体」は、以下の実施の形態ではマルチコプターを例として説明するが、マルチコプター以外の回転翼機等を含んでもよい。

【0023】

実施の形態１
［構成の説明］
図１は、情報処理システムの一例を示すブロック図である。情報処理システム１０は、入力部１１、設定部１２、ブレード計算部１３、時間積分計算部１４、変化計算部１５、表示部１６及び記憶部１７を備える。以下、情報処理システム１０の各部について説明する。

【0024】

入力部１１は、ユーザからの入力情報を受け付ける構成ユニットである。入力部１１は、例えばタッチパネル、キーボード、マウス等の入力インタフェースで構成されていてもよい。ユーザは、入力部１１を用いることで、任意の内容の入力情報を情報処理システム１０に入力することができる。例えば、ユーザは、シミュレーション対象となるマルチコプターの各部の形状及び質量等についての情報や、マルチコプターのブレードの回転条件（例えば、回転速度）を、入力部１１を用いて入力してもよい。また、ユーザは、マルチコプターの周辺に存在する物体（以下、周辺物体とも記載する）のマルチコプターとの相対位置、大きさ及び形状についての情報や、計算に用いられる後述のパラメータについても、入力部１１を用いて入力してもよい。周辺物体は、例えば、地面や建物等である。

【0025】

設定部１２は、シミュレーション対象となるマルチコプターの全体を含む所定の空間の範囲を、シミュレーションにおける計算領域として設定する。空間の範囲は、例えば、入力部１１によって入力された、マルチコプターの大きさ及び形状に基づいて設定される。また、設定部１２は、安定して時間積分を行うことのできる時間ステップの間隔も設定する。

【0026】

また、シミュレーションにおいて必要である場合、設定部１２は、計算領域（所定範囲の空間）に周辺物体を配置する。周辺物体の配置は、例えば、入力部１１によって入力された相対位置、大きさ及び形状についての情報に基づいて設定される。

【0027】

また、設定部１２は、計算領域の境界条件（例えば流出境界条件）などを設定する。さらに、設定部１２は、マルチコプターを、ブレードと、マルチコプターの機体全体からブレードを除いた全ての部分である胴体とに分割する。例えば、機体におけるローターの支持部も、胴体に含まれる。設定部１２は、分割を、入力部１１から入力されたマルチコプターの情報に基づいて実行してもよい。

【0028】

ここで、設定部１２は、計算領域を、ブレード周辺の領域（以下、領域１とも記載）と、領域１以外の計算領域の領域（以下、領域２とも記載）に分割してもよい。領域１は、例えば、マルチコプターのブレード表面から所定長以内に存在する領域である。領域２は、マルチコプターの胴体周辺の領域と、それ以外の計算領域に分割されてもよい。ただし、設定部１２は、計算領域を領域１と領域２に分割しなくともよい。

【0029】

図２は、マルチコプターの一例を示す。マルチコプターＭ１は、ローターＲ１～Ｒ４と、胴体Ｄ１を備える。設定部１２は、マルチコプターＭ１に関しての計算領域を分割する場合には、計算領域を、ローターＲ１～Ｒ４周辺の領域である領域１と、胴体Ｄ１周辺を含む領域である領域２と、に分割してもよい。

【0030】

また、設定部１２は、計算領域を分割する計算格子の形状を設定する。計算格子の形状は、直方体等の六面体、角柱、四面体等、任意の形状を適用することができる。計算格子の形状は、ユーザが入力部１１を用いて入力したものであってもよいし、記憶部１７に記憶されたものであってもよいし、計算に関する情報に基づいて設定部１２が設定したものであってもよい。設定部１２は、十分な精度で効率よく計算できるような格子サイズを設定する。設定部１２は、胴体の近傍（例えば、胴体から所定の距離以内）では細かい（すなわち、小さい）格子を用い、胴体の近傍以外では粗い（すなわち、大きい）格子を用いてもよい。また、設定部１２は、発生する流れにおける渦の強さに応じて格子サイズを変えてもよい。

【0031】

図１に戻り、ブレード計算部１３について説明する。ブレード計算部１３は、以下の処理を実行することで、ブレード周辺の領域について、ブレードが回転することでブレードが周囲の空気から受ける力の計算と、ブレードが周囲の空気に対して与える運動量の計算を実行する。設定部１２が計算領域を領域１と領域２に分割する場合には、ブレード計算部１３は、設定された領域１について上記の計算を実行する。以下、ブレード計算部１３が実行する、アクチュエータラインモデル適用の処理の詳細について記載する。

【0032】

（Ａ）まず、ブレード計算部１３は、マルチコプターの３Ｄ（dimensions）形状モデルからブレードを削除する。

【0033】

（Ｂ）次に、ブレード計算部１３は、ブレードの形状に合わせた翼素データを、削除されたブレードの領域に対して適用する。翼素データは、ブレードの断面形状を示す２Ｄデータとその空力特性のデータである。ブレード計算部１３は、例えば記憶部１７に格納されている翼素データを用いてもよいし、情報処理システム１０の外部（例えばインターネットで提供されている）から、用いる翼素データを取得してもよい。

【0034】

（Ｃ）そして、ブレード計算部１３は、適用するアクチュエータラインモデルのパラメータを具体的に設定する。

【0035】

図３は、アクチュエータラインモデルを適用した場合のブレードの一例を示す。図３において、各ブレードＢ１～Ｂ３は、ブレードを構成する複数の粒子（以下、マーカー粒子とも記載）で代表される。マーカー粒子は、図３では点として示されており、マーカー粒子が並ぶことによって各ブレードが構成される。また、各ブレードＢ１～Ｂ３は、計算格子内に存在する。このとき、マーカー粒子が周囲の計算格子に対して及ぼす影響範囲を決定するガウシアンフィルタηは、以下の通り示される。

【数1】

（１）において、ｄ_ｉは、ｉ番目のマーカー粒子と、（x,y,z）の位置にある計算格子点との距離である。εは、ガウス分布の標準偏差であり、マーカー粒子の影響半径のパラメータである。図３には、マーカー粒子を中心とする半径εの円が示されている。なお、図３は、簡単のためブレード及び計算格子を２Ｄで表現しているが、実際にはこれらは３Ｄで表される。

【0036】

なお、後述の通り、時間積分計算部１４が計算を実行する際には、マルチコプターの物体境界条件が設定される。しかしながら、ブレード周辺の領域に対しては上記の通りアクチュエータラインモデルが設定されるため、ブレード計算部１３は、物体境界条件を設定する必要はない。

【0037】

図４は、ブレードに配置されたマーカー粒子が翼素運動量理論に基づいて周囲の空気から受ける力の一例を示す。ブレードＢ２のマーカー粒子ｉに対してかかる揚力及び抗力は、それぞれＦ_l及びＦ_ｄとして表される。θ、ｚはそれぞれ回転方向と推力の方向の座標軸を示す。γは、ブレードの回転方向であるθ軸に対するブレードの取り付け角度（局所ピッチ角）であり、αは迎角である。このとき、ブレードの相対速度ｕ_ｒｅｌは次のように表される。

【数2】

ここで、ｕ_ｚは周囲の空気の流れのｚ軸成分であり、ｕ_θは周囲の空気の流れのθ軸成分である。また、ｒはマーカー粒子ｉのブレードの回転中心からの距離であり、Ωはブレードの角速度である。ｕ_ｒｅｌを用いると、マーカー粒子の位置でブレードが周囲の空気から受ける力Ｆ_２Ｄは、以下のように表される。

【数3】

ここで、Ｃ_ａはブレードの翼弦長であり、ｅ_ｌ、ｅ_ｄはそれぞれ揚力、抗力が作用する方向の単位ベクトルである。Ｃ_ｌ、Ｃ_ｄはそれぞれ揚力係数、抗力係数であり、迎角αに対するデータセットとして事前に翼素データとして用意される。また、Δｒは、マーカー粒子同士の間隔を示すパラメータである。

【0038】

また、マーカー粒子ｉが、（ｘ,ｙ,ｚ）の位置における計算格子点に与える力（すなわち、ブレードが回転していることで周囲の空気に与える運動量）は、次のように算出される。

【数4】

ここで、Ｆ_２Ｄ,ｉはＦ_２Ｄのマーカー粒子ｉの成分であり、Ｎ_ｐはマーカー粒子の総数である。

【0039】

ブレード計算部１３は、パラメータε、Δｒを設定することで、上記の（１）～（４）の計算を実行できるようにする。パラメータε、Δｒの設定は、ユーザが入力部１１を用いることでなされてもよいし、記憶部１７に格納されている設定をブレード計算部１３が用いることでなされてもよい。（１）～（４）の計算において用いられるパラメータε、Δｒ以外の数値は、マルチコプターの情報等とシミュレーションによって得られるデータに基づいて決定される。

【0040】

（Ｄ）その後、ブレード計算部１３は、適用するアクチュエータラインモデルに対し、マルチコプターの機体の姿勢に応じて、機体とともに並進又は回転軸の回転の少なくともいずれかに関する座標変換を実行する。例えば、ブレード計算部１３は、クオータニオンを用いることで回転軸の回転の座標変換を実現してもよい。

【0041】

（Ｅ）次に、ブレード計算部１３は、マルチコプターの３Ｄモデルを計算領域に配置した場合にアクチュエータラインモデルを適用するブレードを決定する。ブレード計算部１３は、ブレード周辺の領域にブレードの直径の１００分の１の長さ程度の格子サイズとなる計算格子を配置する。格子サイズは任意であるが、格子サイズの大小は計算精度に影響する。計算格子の形状は、設定部１２が設定したものである。

【0042】

（Ｆ）そして、ブレード計算部１３は、アクチュエータラインモデルを適用することで数式（３）により計算されるブレードの受ける力（推力とトルク）を計算し、計算された受ける力を、機体のローターブレードを支持している部分に与える。一方、ブレード計算部１３は、その力の反作用の運動量を、数式（４）に示す通りブレード周囲の空気に与える。

【0043】

（Ｇ）ブレード計算部１３は、（Ｅ）で決定した領域１内の１枚のブレードについて、（Ｆ）でなされた処理の結果を用いて、上記の（１）～（４）の計算を実行する。計算で用いられる計算格子の設定は、（Ｅ）でなされた通りである。その結果、ブレード計算部１３は、機体が移動したり回転したりする影響下でブレードが回転することでブレードが周囲の空気から受ける力と周囲の空気に与える運動量を算出することができる。

【0044】

（Ｈ）ブレード計算部１３は、領域１内の各ブレードについて、（Ｇ）の計算を実行することで、ブレードが周囲の空気から受ける力と周囲の空気に与える運動量を、全ブレードについて算出する。

【0045】

図１に戻り、時間積分計算部１４について説明する。時間積分計算部１４は、胴体全長（すなわち、機体の全長）の３００分の１の長さ程度の格子サイズの計算格子を、機体表面のごく近傍に配置する。設定部１２が計算領域を領域１と領域２に分割する場合には、時間積分計算部１４は、設定された領域２（すなわち、胴体周辺の領域）に領域１と異なる格子サイズの計算格子を配置してもよい。格子サイズは任意であるが、ある程度の計算精度を確保するためには、所定の大きさ以下の細かい格子を配置する必要がある。ただし、格子サイズとして、機体から遠方になるに従い、大きいサイズの格子が用いられてもよい。計算格子の形状は、設定部１２が設定したものである。そして、時間積分計算部１４は、計算格子を配置した領域に、移動する物体（すなわちマルチコプター）の物体境界条件を設定する。

【0046】

時間積分計算部１４は、設定した物体境界条件を考慮しながら、計算格子を用いて離散化された流体方程式に対して、数値計算手法を適用して解くことで空力計算を行う。これにより、時間積分計算部１４は、機体全体からブレードを除いた全ての部分（すなわち胴体）が移動することで胴体が周囲の空気から受ける力を算出する。また、移動する物体を反映する境界条件が空気側に及ぼす影響も計算される。「空力計算の数値計算手法」は、アクチュエータラインモデル以外の公知の計算手法であり、例えば、有限体積法、有限要素法、格子ボルツマン法、粒子法などを含む。この計算手法の詳細についてはよく知られているため、説明を省略する。

【0047】

また、時間積分計算部１４は、必要であれば、計算量を削減するために、スライディング格子、重合格子等の移動計算格子法を用いてもよい。

【0048】

ここまでの過程で、ブレード計算部１３及び時間積分計算部１４は、ある時刻でのブレード及び胴体（すなわち機体全体）が周囲の空気から受ける力を計算することができる。ブレード計算部１３及び時間積分計算部１４は、この計算を、ブレードが回転する一定の時間について繰り返し実行することで、ブレードの回転に応じた空気流、および胴体が移動したことにより発生する周囲の空気流を算出することができる。これにより、マルチコプターの航行の振舞いを確認することが可能となる。なお、ブレード計算部１３は、ブレードの回転に応じた空気流を算出する際に、入力部１１から入力された回転数が時間的に変化する条件を用いてもよい。

【0049】

変化計算部１５は、ブレード計算部１３が計算した、ブレードが周囲の空気から受ける力と、時間積分計算部１４が計算した胴体が周囲の空気から受ける力と、に基づいて、マルチコプターが空中でどのように航行するかをシミュレートする。詳細には、変化計算部１５は、胴体の表面を細かく分割し、時間積分計算部１４で求めた胴体表面に働く力を積分することで、機体全体が受ける平行移動の力と回転のトルクを求める。変化計算部１５は、マルチコプターの各部の大きさ、形状及び質量等の情報や重心周りの慣性モーメントの情報を用いて運動方程式を積分することにより、機体の並進・回転・姿勢の時間変化を計算する。換言すれば、変化計算部１５は、胴体に対する空力計算の数値計算手法の適用結果と、ブレードに対するアクチュエータラインモデルの適用結果と、を用いることでマルチコプターの航行をシミュレートする。

【0050】

表示部１６は、変化計算部１５のシミュレート結果を表示する。表示部１６は、必要であれば、ブレード計算部１３の計算結果又は時間積分計算部１４の計算結果の少なくともいずれかについても、さらに表示してもよい。これらの計算結果は、ブレード周辺や機体の表面近くの空気の流れ、機体表面のさまざまな場所が空気から受ける圧力などである。表示部１６は、例えばディスプレイやタッチパネル等である。ユーザは表示部１６を視認することで、一定の時間経過におけるマルチコプターの自由航行シミュレーションの結果を把握することができる。

【0051】

記憶部１７は、時間積分計算部１４が実行する計算で得られた計算領域内の全ての空気の流れの時間変化、機体表面の圧力分布の時間変化、ブレード計算部１３が実行する計算で得られたブレードが受ける力（推力とトルク）の時間変化を記憶する。記憶する頻度は、設定部１２で指定される。また、記憶部１７は、設定部１２～表示部１６に処理を実行させるためのプログラムや、計算格子の形状、翼素データ、パラメータε及びΔｒの設定など、処理に必要な情報も格納する。

【0052】

［フローの説明］
図５は、情報処理システム１０の代表的な処理の一例を示すフローチャートである。図５のフローチャートによって、情報処理システム１０の処理が説明される。なお、各処理の詳細については上述の通りであるため、説明を省略する。

【0053】

入力部１１は、ユーザの操作に基づいて、シミュレーション対象となるマルチコプターの情報、ブレードの回転条件、周辺物体の情報等、計算に必要な情報を入力する（ステップＳ１１）。設定部１２は、計算領域、周辺物体の配置、境界条件、計算格子の形状といった計算条件を設定する（ステップＳ１２）。なお、設定部１２は、領域１及び２についてさらに設定してもよい。

【0054】

そして、ブレード計算部１３は、アクチュエータラインモデルを適用することで、周囲の空気から受ける力をブレードについて算出する（ステップＳ１３）。この詳細は（Ａ）～（Ｈ）において示した通りである。また、時間積分計算部１４は、空力計算の数値計算手法を適用することで、周囲の空気から機体が受ける力と周囲の空気への影響を算出する（ステップＳ１４）。この詳細も上記の通りである。ブレード計算部１３及び時間積分計算部１４は、ステップＳ１３及びＳ１４の計算を、ブレードが回転する一定の時間経過についてループして繰り返し実行する。ここで、ステップＳ１３とステップＳ１４はいずれが先に実行されてもよいし、両者が並行して実行されてもよい。

【0055】

変化計算部１５は、ブレード計算部１３の計算結果及び時間積分計算部１４の計算結果に基づいて、マルチコプターの航行をシミュレートする（ステップＳ１５）。表示部１６は、変化計算部１５の計算結果を表示する（ステップＳ１６）。

【0056】

［効果の説明］
以上に示した通り、情報処理システム１０は、飛行体のブレードの計算にはアクチュエータラインモデルを適用する一方、それ以外の機体の部分の計算には空力計算の数値計算手法を適用する。そのため、自由航行シミュレータに必要な計算量を劇的に削減することが可能となる。

【0057】

これまでの関連技術では、ブレードの計算に対しても粘性項を持つNavier-Stokes方程式を離散化して解く空力計算の数値計算手法が適用されていた。このとき、十分な精度の計算結果を得るために、ブレード近傍には、ブレードの直径の１０００分の１以下の長さの格子サイズの計算格子が配置される。ここで、ブレード近傍に配置される計算格子がブレードの直径の１０００分の１のサイズである場合（状況Ａ）、配置される計算格子がブレードの直径の１００分の１のサイズである場合（状況Ｂ）と比較して、計算対象として考慮すべき計算格子の数は１０^３（１０００）倍となる。さらに、状況Ａは、状況Ｂと比較して、計算の時間ステップの間隔も１０分の１とする必要がある。したがって、例えば領域１についての計算量は、状況Ａは、状況Ｂと比較して、１０^４（１００００）倍となる。また、ブレードの回転数が早いほど、計算の正確性を担保するために、計算格子を小さくする必要が生じる。そのため、ブレードについての計算量が膨大なものとなっていた。

【0058】

一方、本発明において、ブレード計算部１３はアクチュエータラインモデルを適応することで、ブレード近傍に配置される計算格子を、ブレードの直径の１０００分の１のサイズよりも大きくすることができる。そのため、関連技術と比較して、ブレードの１ステップ当たりの計算負荷を低下させることができる。例えば、関連技術における状況Ａを、本発明では状況Ｂに変化させることができる。その結果、領域１に関する計算量を１００００分の１以下に減少させることができ、領域１に関する計算量を例えば自由航行シミュレータ全体の計算量の１０分の１以下とすることも可能となる。

【0059】

したがって、本発明に係る情報処理システム１０は、ブレードが１０００回転する程度までの時間を十分計算できることから、機体全体の自由航行をシミュレーション可能とすることができる。例えば、飛行体の離着陸の過程についても、情報処理システム１０は、自由航行シミュレーションが可能である。

【0060】

さらに、以上に示した情報処理システム１０の手法は、計算の精度も担保することができる。以下、実際のシミュレータ結果のデータを参照して、情報処理システム１０の計算の正確性について説明する。

【0061】

図６は、参照論文（Knut Erik Teigen Giljarhus, Alessandro Porcarelli and Jorgen Apeland, “Investigation of Rotor Efficiency with Varying Rotor Pitch Angle for a Coaxial Drone”, Drones 2022, April 4, 2022, MDPI）で使用されたブレードと同じ形状と回転数のブレードを用いた場合の、周囲の空気が影響を受けることで生じる、ある時刻の流れの速度の大きさを示す数値計算のデータである。数値計算のデータは、ブレードの回転速度が１６００ｒｐｍ、１９００ｒｐｍ、２２００ｒｐｍ、２５００ｒｐｍの４種類の計算結果がある。図６から見てとれるように、アクチュエータラインモデルを用いた本発明に係る手法は、非常に細かい格子を用いてNavier-Stokes方程式を離散化して解く空力計算による結果（図６の上段）と非常によく一致する計算結果（図６の下段）を得ることができる。

【0062】

図７は、参照論文と同じブレードの形状と回転数を用いて情報処理システム１０が計算を実行した結果（すなわち、本発明に係るアクチュエータラインモデルを用いた計算結果）を表したグラフである。図７には、計算結果として、ブレードの回転速度が１６００ｒｐｍ、１９００ｒｐｍ、２２００ｒｐｍ、２５００ｒｐｍ、２６００ｒｐｍとなる５種類の計算結果が示されている。図７に示される通り、参照論文の実験データと、非常に細かい格子を用いてNavier-Stokes方程式を離散化して解く空力計算による結果と、が推力に関して非常によく一致する計算結果を得ることができる。

【0063】

図８は、空飛ぶクルマが離陸する際に、アクチュエータラインモデルを用いて計算したブレードの推力の反作用が周囲の空気にダウンウォッシュとして与えられ、それが地面に当たり横に広がる状態を示す計算結果である。図８では発生する渦の様子が示されていて、空気の流れの速度勾配テンソルの第二不変量の等値面が示す通り、情報処理システム１０は、ブレードから非常に細かい渦が放出されることを計算することが可能となる。

【0064】

なおユーザは、ブレードに対してアクチュエータラインモデルを適用するときの、マーカー粒子の間隔、又はガウシアンフィルタの影響半径のパラメータ、の少なくともいずれかの情報を入力してもよい。ブレード計算部１３は、アクチュエータラインモデルの適用において、入力された情報を使用することができる。したがって、ユーザは、アクチュエータラインモデルの計算量の削減度合い及び正確性を、入力パラメータを変化させることで、自由に調整することができる。

【0065】

また、ブレード計算部１３は、ブレード周辺の領域に適用される計算格子のサイズを、ブレードの直径の１０００分の１以上としてもよい。例えば、ブレード計算部１３は、計算格子のサイズを、ブレードの直径の１００分の１としても良い。この値にすることで、ブレード計算部１３は、アクチュエータラインモデルを用いた計算の精度を保ちつつ、計算量の削減の効果を図ることができる。しかしながら、ブレード計算部１３は、必要に応じて、計算格子の長さをそれ以外の大きさとしてもよい。例えば、ブレード計算部１３は、計算格子のサイズを、ブレードの直径の５０分の１の大きさ、２００分の１の大きさなど任意のサイズとしてもよい。

【0066】

本発明の情報処理システム１０は、単一のコンピュータ装置として構成されてもよいし、複数のコンピュータ装置を有する分散システムとして構成されてもよい。分散システムにおいて、情報処理システム１０が実行する処理は、複数のコンピュータ装置によって分担されて実行することができる。つまり、入力部１１～記憶部１７の構成要素は、２台以上のコンピュータ装置に分散されて搭載されてもよい。

【0067】

以上に示した実施の形態では、本発明に係る情報処理システムをハードウェアの構成として説明したが、本発明に係る情報処理システムは、これに限定されるものではない。本発明は、上述の実施形態において説明された情報処理システム１０を構成する各装置の処理を、コンピュータ内のプロセッサにコンピュータプログラムを実行させることにより実現することも可能である。

【0068】

図９は、実施の形態に係る情報処理システム（換言すればコンピュータ）のハードウェア構成例を示すブロック図である。図９を参照すると、情報処理システム９０は、信号処理回路９１、プロセッサ９２、メモリ９３、ストレージ９４及びインタフェース９５を含む。

【0069】

信号処理回路９１は、プロセッサ９２の制御に応じて、信号を処理するための各種の回路である。

【0070】

プロセッサ９２は、メモリ９３と接続されており、メモリ９３からコンピュータプログラムを読み出し、メモリ９３と通信しながら実行することで、上述の実施形態において説明されたシステムの処理を行う。プロセッサ９２の一例として、ＣＰＵ（Central Processing Unit）、ＭＰＵ（Micro Processing Unit）、ＧＰＵ（Graphics Processing Unit）、ＦＰＧＡ（Field-Programmable Gate Array）、ＤＳＰ（Demand-Side Platform）、ＡＳＩＣ（Application Specific Integrated Circuit）のうち一つを用いてもよいし、そのうちの複数を並列で用いてもよい。

【0071】

メモリ９３は、例えば揮発性メモリで構成される主記憶装置である。メモリ９３は、１個に限られず、複数設けられてもよい。なお、揮発性メモリは、例えば、ＤＲＡＭ (Dynamic Random Access Memory)等のＲＡＭ (Random Access Memory)であってもよい。

【0072】

メモリ９３は、１以上の命令とデータを格納するために使用される。ここで、１以上の命令は、プログラムとしてメモリ９３に格納される。プロセッサ９２は、これらのプログラムとデータをメモリ９３から読み出して実行することで、上述の実施形態において説明された処理を行うことができる。

【0073】

なお、メモリ９３は、プロセッサ９２の外部に設けられるものに加えて、プロセッサ９２に内蔵されているものを含んでもよい。また、メモリ９３は、プロセッサ９２を構成するプロセッサから離れて配置されたストレージを含んでもよい。この場合、プロセッサ９２は、Ｉ/Ｏ（Input/Output）インタフェースを介してメモリ９３にアクセスすることができる。

【0074】

ストレージ９４は、例えば不揮発性メモリで構成される補助記憶装置である。ストレージ９４は、１個に限られず、複数設けられてもよい。なお、不揮発性メモリは、例えば、ＨＤＤ（Hard Disk Drive）、ＰＲＯＭ (Programmable Random Only Memory)、ＥＰＲＯＭ (Erasable Programmable Read Only Memory) 等のＲＯＭ (Read Only Memory)、フラッシュメモリや、ＳＳＤ（Solid State Drive）であってもよい。ストレージ９４には、メモリ９３に供給されるプログラムが格納されている。また、ストレージ９４は、実施の形態１における記憶部１７としても機能し、上記に示した時間積分計算部１４、変化計算部１５の計算結果や、処理に必要な情報を記憶してもよい。

【0075】

インタフェース９５は、ネットワークを介して信号又はデータを送受信するための通信回路を含む。インタフェース９５は、例えばＮＩＣ（Network Interface Card）であってもよい。プロセッサ９２は、メモリ９３およびストレージ９４に格納されたデータを、インタフェース９５を介して他の情報処理システムに送信してもよいし、他の情報処理システムから送信されたデータを、インタフェース９５を介してメモリ９３およびストレージ９４に格納してもよい。

【0076】

以上に説明したように、上述の実施形態における各システムが有する１又は複数のプロセッサは、図面を用いて説明されたアルゴリズムをコンピュータに行わせるための命令群を含む１又は複数のプログラムを実行する。プログラムが実行されることで、各実施の形態に記載された情報処理が実現できる。

【0077】

プログラムは、コンピュータに読み込まれた場合に、実施形態で説明された１又はそれ以上の機能をコンピュータに行わせるための命令群又はソフトウェアコードを含む。プログラムは、非一時的なコンピュータ可読媒体又は実体のある記憶媒体に格納されてもよい。限定ではなく例として、コンピュータ可読媒体又は実体のある記憶媒体は、random-access memory（RAM）、read-only memory（ROM）、hard disk drive(HDD)、フラッシュメモリ、solid-state drive（SSD）又はその他のメモリ技術、compact disc read-only memory（CD-ROM）、digital versatile disk（DVD）、Blu-ray（登録商標）ディスク又はその他の光ディスクストレージ、磁気カセット、磁気テープ、磁気ディスクストレージ又はその他の磁気ストレージデバイスを含む。プログラムは、一時的なコンピュータ可読媒体又は通信媒体上で送信されてもよい。限定ではなく例として、一時的なコンピュータ可読媒体又は通信媒体は、電気的、光学的、音響的、またはその他の形式の伝搬信号を含む。一時的なコンピュータ可読媒体又は通信媒体は、電線及び光ファイバ等の有線通信路、又は無線通信路を介して、プログラムをコンピュータに供給できる。

【0078】

以上、実施の形態を参照して本開示を説明したが、本開示は上述の実施の形態に限定されるものではない。本開示の構成や詳細には、本開示のスコープ内で当業者が理解し得る様々な変更をすることができる。そして、各実施の形態は、適宜他の実施の形態と組み合わせることができる。

【符号の説明】

【0079】

１０ 情報処理システム
１１ 入力部 １２ 設定部
１３ ブレード計算部 １４ 時間積分計算部
１５ 変化計算部 １６ 表示部
１７ 記憶部
９０ 情報処理システム
９１ 信号処理回路 ９２ プロセッサ
９３ メモリ ９４ ストレージ
９５ インタフェース

【要約】

【課題】必要な計算量を削減することに寄与するシミュレータの開発方法、情報処理システム及びプログラムを提供すること。
【解決手段】本開示の一実施の形態にかかるシミュレータの開発方法は、飛行体の胴体周辺の空間を、計算格子を用いて離散化し、離散化された流体方程式を解くことで空力計算の数値計算手法を適用するステップと、飛行体のブレードに対してアクチュエータラインモデルを適用するステップと、数値計算手法の適用結果と、ブレードに対するアクチュエータラインモデルの適用結果と、を用いることで飛行体の航行をシミュレートするステップと、をコンピュータが実行するものである。
【選択図】図５

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】