サイクロローターの気流を可視化する
最近注目されているサイクロローター式eVTOL機について、回転翼がどういう動きをするものか不思議に思ったので少し調べてみました。元々はオーストリアの技術者エルンスト・シュナイダー氏が実用化したもので、フォイト・シュナイダープロペラと呼ばれ船舶分野で先行していた推力偏向技術です。理論そのものは古くから知られていたもののあまり評判は良くなかったらしく、高度な機構制御ができる時代になってようやく航空分野に持ち込まれたようです。動作原理についてはWikipediaをご覧ください。Cyclorotor - Wikipediaen.wikipedia.org典型的な設計では4枚あるいは6枚の翼が円周上に配置され、クランク機構によって特定の角度で迎角が大きく振れるようになっています。これによって水泳のクロールのように空気を力強く掻き込むような動きを実現しており、推力に寄与しないストロークポジションでは抗力を最小化するよう翼列が水平となります。ネット検索すると論文がいくつか出てくるのですが、興味深い研究が見つかったので読んでみました。それがポーランドのシレジア工科大学動力工学ターボ機械学科による研究です。Influence of Geometrical Parameters on the Shape of the Cycloidal Function Curve of a Fan with a Cycloidal Rotorこの研究では翼弦長50mmの翼型CLARK Yを使用して風洞実験を行っています。諸元がすべて明らかになっているので、OpenFOAMでのシミュレーションに利用できそうです。諸元の確認クランク機構を理解するためには、各部の角度と寸法を規定しなければいけません。この論文では以下のように定義してあります。Oは回転中心、Pは翼の旋回中心で翼弦長50%位置、Rはローター半径、MはコンロッドLの取り付け位置です。eがオフセットリンクで、この角度が推力方向の制御に用いられます。eの角度と機構全体が生ずる推力/揚力方向との関係性を求めるには多くの計算が必要のように思います。aは途中計算に必要なP-E間の距離を示し、以下の数式で求まります。途中式をすっ飛ばして、必要な翼の角度θは以下のサイクロイド関数で求まります。回転数が与えられたとき、ある時刻における4枚の翼それぞれの迎角が求まります。論文によるとLabVIEWで各部寸法を検討した結果、R=70mmにおいてe=13mm、ε=292deg、d=27mm、L=74mmとなったとのことです。これらの数値を無からひねりだす過程が彼らの研究のハイライトと言えるでしょう。残るaは式(3)にローター角度ψを与えることで求まります。式(6)に出揃った変数を代入すれば、欲しかった角度θが求まります。シミュレーションの準備各要素の位置関係が分かったので準備を始めます。今回はESI版のOpenFOAM v2406を使います。インストールの方法はこちらをご覧ください。ゼロからケースファイルを作成するのは大変なので、テンプレートとなるケースを導入します。イタリアはジェノバ大学からのスピンオフ企業Wolf Dynamics社が公開しているオーバーセット・メッシュのチュートリアルをこちらから辿って入手します。ファイル名はdynamicMeshes_2025_clean.tar.gzです。この圧縮ファイルを解凍し、ケースファイルfoil_shmmeshと座標生成ユーティリティgen6DoFを$FOAM_RUNにコピーして動作することを確認しておきます。このケースはふらつきのある翼に生じる力係数を求めるためのものですが、オーバーセット・メッシュ対応の非定常圧縮性流体ソルバーpimpleDyMFoamを利用しており今回のシミュレーションの出発点とするのに最適であろうと判断しました。なお、ケースファイルの実行に先立ちrun_solver.shの最後にreconstructParコマンドが抜けているので追記する必要があります。このチュートリアルでは領域左側のinletから流速1m/sの空気を流していますが、今回のシミュレーションでは流速Uを0m/s、すなわち無風状態でローターがどのように空気をかき混ぜるかを見ます。もちろん、流速Uを書き換えればローターに風を当てることもできます。参照した論文ではローター回転速度1000rpmで風洞実験していますが、迎角を計算しやすくするため1200rpmでシミュレーションします。他の論文ではほぼ同じ大きさのローターを1400rpmで回転させて模型を飛ばしていることから、この辺りの回転数で検証するのが妥当であろうと思います。翼の周速度V[m/s]は以下の公式で求めます。直径D[mm]、回転数N[rpm]とすると、V = π * D * N / (1000 * 60) = π * 140 * 1200 / (1000 * 60) = 8.79[m/s]となります。翼の代表長さ0.05m、高度0kmのとき、レイノルズ数は30094となります。ちなみにヤマトホールディングス社とサイクロテック社が共同開発している中型eVTOL機CCY-01の資料によれば、ローター半径0.25m、回転数2600rpm、巡航速度120km/hとのことです。この場合翼の周速度は68.07m/sとなり、これが120km/hすなわち33.33m/sで飛行すると翼の速度は101.4m/sにもなるのでいよいよ空気を圧縮性流体とみなす必要があります。翼型STLモデル作成FreeCADを使ってこちらの手順でCLARK Y翼型を作成します。翼弦上の中央をモデル原点とし、翼弦長50mm、TOP方向すなわちZ軸方向に±125mm(翼幅250mm)押し出しておきます。TOP方向から見て0deg、180deg、90deg、270degの順番で回転させた4パターンを作成し、それぞれclarky1.ast、clarky2.ast、clarky3.ast、clarky4.astを書き出し、拡張子を.stlに変更します。できたSTLファイルをエディタで開き、行頭のsolidと行末のendsolidをMeshからairfoilに書き換えます。この辺りの事情はこちらとこちらをご覧ください。ケースファイルの編集フォルダfoil_shmmeshをcycloidにリネームし、以後このケースファイルcycloidを編集していきます。フォルダairfoilをairfoil1にリネームし、これと同じ内容のフォルダairfoil2、airfoil3、airfoil4を複製します。続いて、これらのフォルダ配下のconstant/triSurfaceに入っていたairfoil.stlを削除し、先ほど作成したclarky*.stlとそれぞれ番号を合わせて差し替えます。airfoil*の編集まずairfoil*の初期位置を示します。水色がairfoil1、灰色がairfoil2、橙色がairfoil3、茜色がairfoil4です。これらはセルセット番号、ゾーンIDとも数字を一致させます。airfoil*のblockMeshDictでは、翼弦長50mmの翼型を覆うメッシュを切ります。長辺80mm、短辺40mmとします。airfoil1とairfoil2は横長、airfoil3とairfoil4は縦長のそれぞれ同じメッシュ寸法になります。例)airfoil1, airfoil2vertices( (-0.04 -0.02 -0.025) ( 0.04 -0.02 -0.025) ( 0.04 0.02 -0.025) (-0.04 0.02 -0.025) (-0.04 -0.02 0.025) ( 0.04 -0.02 0.025) ( 0.04 0.02 0.025) (-0.04 0.02 0.025));blocks( hex (0 1 2 3 4 5 6 7) (40 20 1) simpleGrading (1 1 1));例)airfoil3, airfoil4vertices( (-0.02 -0.04 -0.025) ( 0.02 -0.04 -0.025) ( 0.02 0.04 -0.025) (-0.02 0.04 -0.025) (-0.02 -0.04 0.025) ( 0.02 -0.04 0.025) ( 0.02 0.04 0.025) (-0.02 0.04 0.025));blocks( hex (0 1 2 3 4 5 6 7) (20 40 1) simpleGrading (1 1 1));続いて、airfoil*のsnappyHexMeshが参照するSTLファイル名とnameを書き換えます。nameはwingからwing*となります。例)geometory{ "clarky1.stl" { type triSurfaceMesh; name wing1; }}また、locationInMeshの値はairfoil*全てが異なります。これらは軸対称に配置されるので、翼型を覆うメッシュの隅をどれも同じように狙っています。airfoil1: (-0.039 0 0)airfoil2: (0.039 0 0)airfoil3: (0 -0.039 0)airfoil4: (0 0.039 0)フォルダairfoil*配下のファイル編集はこれで終わりです。できたメッシュは後の工程でしか見られませんが、参考までにメッシュを切ったairfoil1を見てみます。メッシュの中心がY軸方向に70mmオフセットされていることが分かります。backGround/0, backGround/0_orgの編集先ほどサーフェスメッシュ名wingをwing*と書き換えました。これに伴い、フォルダ0と0_orgで定義されているwingについて、境界条件はそのままにwing*を書き加えます。ファイルを見れば何をすべきか分かると思いますが、修正箇所が多いので見落としのないようにしてください。なおファイルzoneIDは他のファイルと定義の仕方が異なっているので、以下のように修正します。"(inlet|outlet|topAndBottom|wing1|wing2|wing3|wing4)"{ type zeroGradient;}backGround/systemの編集領域backGroundはX軸、Y軸とも0.5mの正方形、Z軸は0.05mのメッシュとします。blockMeshDictを以下のように修正します。vertices( (-0.25 -0.25 -0.025) ( 0.25 -0.25 -0.025) ( 0.25 0.25 -0.025) (-0.25 0.25 -0.025) (-0.25 -0.25 0.025) ( 0.25 -0.25 0.025) ( 0.25 0.25 0.025) (-0.25 0.25 0.025));blocks( hex (0 1 2 3 4 5 6 7) (100 100 1) simpleGrading (1 1 1));続いて、controlDictを以下のように編集します。修正すべき行はバラバラですが列挙します。endTime 1;writeInterval 0.001;maxCo 4; // 最大クーラン数。ここでは変更しませんが必要に応じて変更のことpatches ("wing1" "wing2" "wing3" "wing4"); // この行は2か所ありますlibs ("libutilityFunctionObjects.so"); // 242行目、旧書式のため新書式に更新最大クーラン数は、条件を試行錯誤するうちは6程度まで上げて計算を早く済ませ、うまく解析できそうなら段階的に下げます。時間がかかっても信頼できる結果を得たいのであれば1まで下げることになります。続いて、setFieldDictとtopoSetDict_movingZoneを編集します。チュートリアルではゾーンは背景の0と翼の1しかありませんでしたが、翼が4つに増えたのでzone IDは4まで必要です。もちろんセルセットもc4、movingZoneも4まで必要です。これらの修正方法はファイルを見れば分かるので例示しません。続いて、topoSetDictを編集します。このファイルでzone IDを割り当てる領域を特定します。insidePointsは翼の初期角度を多少振っても領域を逸脱しない場所を指定しています。actionsだけ抜き出します。actions( { // 全領域をc0とする name c0; type cellSet; action new; source regionToCell; sourceInfo { insidePoints ((-0.249 -0.249 0)); } } { // c0を複製してc1を作成する name c1; type cellSet; action new; source cellToCell; sourceInfo { set c0; } } { // c0以外の領域をc1とする name c1; type cellSet; action invert; } { // c1の中にc2を作成する name c2; type cellSet; action new; source regionsToCell; sourceInfo { // c2領域を指定 insidePoints ((0 -0.065 0)); set c1; } } { // c1からc2を除外する name c1; type cellSet; action delete; source cellToCell; sourceInfo { set c2; } } { // c1の中にc3を作成する name c3; type cellSet; action new; source regionsToCell; sourceInfo { // c3領域を指定 insidePoints ((-0.065 0 0)); set c1; } } { // c1からc3を除外する name c1; type cellSet; action delete; source cellToCell; sourceInfo { set c3; } } { // c1の中にc4を作成 name c4; type cellSet; action new; source regionsToCell; sourceInfo { // c4領域を指定 insidePoints ((0.065 0 0)); set c1; } } { // c1からc4を除外する name c1; type cellSet; action delete; source cellToCell; sourceInfo { set c4; } });以上でairfoil*の編集で示した画像の通りzone IDが割り当てられます。backGround/constantの編集まず、dynamicMeshDictを開くとmovingZoneが設定してあります。翼4つ分必要ですので、movingZone1〜movingZone4まで用意します。また、各翼のモーション制御を個別で行いますので、後に生成するモーション座標データ6DoF1.dat〜6DoF4.datを各movingZoneに指定します。設定項目内のCofGはCenter of Gravityすなわち重心座標ですが、ここではモーションの原点座標になります。これは各movingZoneごとに座標が異なります。movingZone1: CofG (0 0.07 0)movingZone2: CofG (0 -0.07 0)movingZone3: CofG (-0.07 0 0)movingZone4: CofG (0.07 0 0)それ以外の編集箇所はファイルを見れば分かりますので例示しません。チュートリアルで使用したモーションデータ6DoF.datは必要ありませんので削除してください。代わりにここに置く翼4つ分のモーションデータをこれから生成します。gen6DoFの編集フォルダcycloidと共に扱うフォルダgen6DoFについて説明します。gen6DoFは、6自由度のモーションデータを生成するプログラムです。C++言語でユーザーオリジナルのプログラムを組んでコンパイル&ビルドすることができます。このgen6DoFを実行するとデータファイルを書き出します。フォーマットは行頭でデータ点数を宣言した後、各行は時刻 - 並進移動量 - 回転量の順となります。テキストエディタでデータファイルを見るとこのようになっています。数学関数にOpenFOAM独自のものを使用するため、その流儀に従ってプログラミングする必要があります。このソースコードはローター回転数1200rpmを前提として1mSecごと1秒間分のモーションデータ1000点を出力します。このケースで使用するソースコードgen6DoF.cを以下に示します。#include "List.H"#include "vector.H"#include "Vector2D.H"#include "Tuple2.H"#include "OFstream.H"#include "mathematicalConstants.H"using namespace Foam;// Variables for this case:// Rotor radius [m]const scalar R = 0.070;// Tension member length [m]const scalar L = 0.074;// Tension member fixed length from blade pivoting point [m]const scalar d = 0.027;// Shift arm length [m]const scalar e = 0.013;// Eccentric phase angle ε [deg]const scalar epsilon = 292;// Blade inclination angle θ [deg]scalar theta;// Position of the rotor blade ψ [deg]//scalar psi;// Distance between shift arm and blade pivoting point [m]scalar a;scalar angleCalc1(const scalar& rotAmpZ, const scalar& t){ scalar psi_rad = ((rotAmpZ * t * 20) + 270) * M_PI/180; scalar epsilon_rad = epsilon * M_PI/180; scalar theta_rad; scalar a_sq; a_sq = e*e + R*R - 2*e*R * Foam::cos(psi_rad + epsilon_rad + M_PI/2); a = Foam::sqrt(a_sq); theta_rad = (M_PI/2) - Foam::asin((e/a) * Foam::cos(psi_rad + epsilon_rad)) - Foam::acos((a*a + d*d - L*L)/(2*a*d)); theta = theta_rad * 180/M_PI; return rotAmpZ * t * 20 + theta;}scalar angleCalc2(const scalar& rotAmpZ, const scalar& t){ scalar psi_rad = ((rotAmpZ * t * 20) + 90) * M_PI/180; scalar epsilon_rad = epsilon * M_PI/180; scalar theta_rad; scalar a_sq; a_sq = e*e + R*R - 2*e*R * Foam::cos(psi_rad + epsilon_rad + M_PI/2); a = Foam::sqrt(a_sq); theta_rad = (M_PI/2) - Foam::asin((e/a) * Foam::cos(psi_rad + epsilon_rad)) - Foam::acos((a*a + d*d - L*L)/(2*a*d)); theta = theta_rad * 180/M_PI; return rotAmpZ * t * 20 + theta;}scalar angleCalc3(const scalar& rotAmpZ, const scalar& t){ scalar psi_rad = ((rotAmpZ * t * 20) + 0) * M_PI/180; scalar epsilon_rad = epsilon * M_PI/180; scalar theta_rad; scalar a_sq; a_sq = e*e + R*R - 2*e*R * Foam::cos(psi_rad + epsilon_rad + M_PI/2); a = Foam::sqrt(a_sq); theta_rad = (M_PI/2) - Foam::asin((e/a) * Foam::cos(psi_rad + epsilon_rad)) - Foam::acos((a*a + d*d - L*L)/(2*a*d)); theta = theta_rad * 180/M_PI; return rotAmpZ * t * 20 + theta;}scalar angleCalc4(const scalar& rotAmpZ, const scalar& t){ scalar psi_rad = ((rotAmpZ * t * 20) + 180) * M_PI/180; scalar epsilon_rad = epsilon * M_PI/180; scalar theta_rad; scalar a_sq; a_sq = e*e + R*R - 2*e*R * Foam::cos(psi_rad + epsilon_rad + M_PI/2); a = Foam::sqrt(a_sq); theta_rad = (M_PI/2) - Foam::asin((e/a) * Foam::cos(psi_rad + epsilon_rad)) - Foam::acos((a*a + d*d - L*L)/(2*a*d)); theta = theta_rad * 180/M_PI; return rotAmpZ * t * 20 + theta;}int main(int argc, char *argv[]){ // End-time of the table [sec] const scalar endTime = 1; // Number of entries in the table const label nTimes = 1000; // Amplitude of the translation [m] const vector transAmp(0.07, 0.07, 0); // Frequency of the translation [rad/s] const vector transOmega(125.6637, 125.6637, 0); // 1200 rpm // Amplitude of the rotation [deg] const vector rotAmp(0, 0, 360); // Frequency of the rotation [rad/s] const vector rotOmega(0, 0, 0); List<Tuple2<scalar, Vector2D<vector>>> timeValues(nTimes); forAll(timeValues, i) { scalar t = (endTime*i)/(nTimes - 1); timeValues[i].first() = t; timeValues[i].second()[0] = vector ( transAmp.x()*Foam::sin(transOmega.x()*t) * -1, transAmp.y()*Foam::cos(transOmega.y()*t) - 0.07, transAmp.z()*Foam::sin(transOmega.z()*t) ); timeValues[i].second()[1] = vector ( rotAmp.x()*Foam::sin(rotOmega.x()*t), rotAmp.y()*Foam::sin(rotOmega.y()*t), angleCalc1(rotAmp.z(), t) ); } { OFstream dataFile("6DoF1.dat"); dataFile.precision(6); dataFile.setf(std::ios::fixed); dataFile << timeValues << endl; } forAll(timeValues, i) { scalar t = (endTime*i)/(nTimes - 1); timeValues[i].first() = t; timeValues[i].second()[0] = vector ( transAmp.x()*Foam::sin(transOmega.x()*t), transAmp.y()*Foam::cos(transOmega.y()*t) * -1 + 0.07, transAmp.z()*Foam::sin(transOmega.z()*t) ); timeValues[i].second()[1] = vector ( rotAmp.x()*Foam::sin(rotOmega.x()*t), rotAmp.y()*Foam::sin(rotOmega.y()*t), angleCalc2(rotAmp.z(), t) ); } { OFstream dataFile("6DoF2.dat"); dataFile.precision(6); dataFile.setf(std::ios::fixed); dataFile << timeValues << endl; } forAll(timeValues, i) { scalar t = (endTime*i)/(nTimes - 1); timeValues[i].first() = t; timeValues[i].second()[0] = vector ( transAmp.x()*Foam::cos(transOmega.x()*t) * -1 + 0.07, transAmp.y()*Foam::sin(transOmega.y()*t) * -1, transAmp.z()*Foam::sin(transOmega.z()*t) ); timeValues[i].second()[1] = vector ( rotAmp.x()*Foam::sin(rotOmega.x()*t), rotAmp.y()*Foam::sin(rotOmega.y()*t), angleCalc3(rotAmp.z(), t) ); } { OFstream dataFile("6DoF3.dat"); dataFile.precision(6); dataFile.setf(std::ios::fixed); dataFile << timeValues << endl; } forAll(timeValues, i) { scalar t = (endTime*i)/(nTimes - 1); timeValues[i].first() = t; timeValues[i].second()[0] = vector ( transAmp.x()*Foam::cos(transOmega.x()*t) - 0.07, transAmp.y()*Foam::sin(transOmega.y()*t), transAmp.z()*Foam::sin(transOmega.z()*t) ); timeValues[i].second()[1] = vector ( rotAmp.x()*Foam::sin(rotOmega.x()*t), rotAmp.y()*Foam::sin(rotOmega.y()*t), angleCalc4(rotAmp.z(), t) ); } { OFstream dataFile("6DoF4.dat"); dataFile.precision(6); dataFile.setf(std::ios::fixed); dataFile << timeValues << endl; } Info<< "End\n" << endl; return 0;}このプログラムでは、翼が70mmの半径で回転する動作にZ軸の回転を使わず、X軸とY軸の並進運動で表現しています。Z軸の回転は翼の角度を決めるのに使い、サブルーチンangleCalc*()で迎角θを計算して翼の回転角度ψに加算しています。angleCalc3()で計算したairfoil3の絶対角度をグラフに描くと以下のようになります。翼の水平姿勢との相対角度すなわち迎角θに注目した論文のグラフでは-cos波とそう変わらない波形に見えましたが、領域との絶対角度ψ+θで見るとクセのある動きをしていることが分かります。gen6DoF.cは以下のコマンドでコンパイル&ビルドします。$ cd $FOAM_RUN$ cd gen6DoF$ wclean$ wmakeこれで$FOAM_USER_APPBINに実行ファイルgen6DoFが生成されます。$FOAM_USER_APPBINとはディレクトリで言うと~/OpenFOAM/user-v2406/platforms/linux64GccDPInt32Opt/binとなっています。では、このプログラムを実行します。$ cd $FOAM_USER_APPBIN$ ./gen6DoF同じ場所に6DoF1.dat, 6DoF2.dat, 6DoF3.dat, 6DoF4.datが一気に生成されます。これをbackGround/constantにコピー&ペーストします。decomposeParDictの編集ケースファイルの至る所にCPUコアに領域を割り付けるdecomposeParDictがあります。これについてはユーザーの計算機によってnumberOfSubdomainsの並列数(=CPUコア数)は違いますし、XYZ軸それぞれの分割数nの設定も違ったものになります。解析にはとにかく時間が必要ですので、計算機の性能を活かせるように設定してください。スクリプトの編集コマンドをバッチ処理するBashシェルスクリプトについて説明します。初期状態でrun_all.sh, run_mesh.sh, run_solver.shの3つのスクリプトがあります。run_all.shは全てのプロセスを通しで実行するスクリプトで、実態はただrun_mesh.shとrun_solver.shを呼び出しているだけす。run_mesh.shはメッシュ作成のみ、run_solver.shは解析のみ行います。スクリプトを分けることで、うまく解析できない時に工程を切り分けて状況を確認することができます。run_solver.shでは変数nprocsで計算に使うコア数を指定しています。これもできるだけ多くのコア数を設定しておくと良いでしょう。run_mesh.shのコマンド列はこのケース独特のものになりますので以下に示します。#!/bin/bashcd airfoil1blockMeshsnappyHexMesh -overwrite | tee log.snappyHexMeshextrudeMeshcreatePatch -overwritetransformPoints -translate '(0 0.07 0)'cd ..cd airfoil2blockMeshsnappyHexMesh -overwrite | tee log.snappyHexMeshextrudeMeshcreatePatch -overwritetransformPoints -translate '(0 -0.07 0)'cd ..cd airfoil3blockMeshsnappyHexMesh -overwrite | tee log.snappyHexMeshextrudeMeshcreatePatch -overwritetransformPoints -translate '(-0.07 0 0)'cd ..cd airfoil4blockMeshsnappyHexMesh -overwrite | tee log.snappyHexMeshextrudeMeshcreatePatch -overwritetransformPoints -translate '(0.07 0 0)'cd ..cd backGroundblockMeshmergeMeshes . ../airfoil1 -overwritemergeMeshes . ../airfoil2 -overwritemergeMeshes . ../airfoil3 -overwritemergeMeshes . ../airfoil4 -overwritetopoSettopoSet -dict system/topoSetDict_movingZonerm -r 0cp -r 0_org 0checkMesh | tee log.checkMeshsetFields | tee log.setFieldscd ..スクリプトの編集は以上です。airfoil*はそれぞれメッシュを作成して領域を初期位置上下左右に70mm移動し、その結果をbackGroundのメッシュにマージします。このブログの最後で述べますが、このスクリプトで本来しておくべきtransformPoints -rotate-z <deg>は技術的な理由で省略しています。解析の実行今回の解析では手始めにinletから風を送り込まず無風状態0m/sでどう空気が流れるかを見ますので、backGround/0/UとbackGround/0_org/Uの2つのファイルを開き、以下のように編集します。internalField uniform (0 0 0);これでケースファイルの準備が全て整いました。以下のコマンドで解析を開始し、ParaViewで可視化します。$ cd $FOAM_RUN$ cd cycloid$ sh run_all.sh$ cd backGround$ paraFoamparaFoamで解析結果を可視化すると、無風の中1200rpmで回転するローターを1mSec刻みで観察することができます。流速Uを表示する際、ColoringのRescale to custom rangeを0〜20[m/s]に狭めると見やすくなります。しかしながら、翼をよく見ると本体が青っぽく表示され、周囲の色に溶け込んで外形が見づらくなっています。画像にいくつかフィルタを入れて見やすく加工します。Pipeline BrowserのbackGround.foamを選択して、SliceをZ Normal方向、ColoringをUとして追加します。OpenFOAMはこの解析を2Dとして計算しましたが、ParaViewはこれを3Dモデルと認識しており、何も処置しなければユーザーは六面体の側面を見ることになります。このまま翼形状をくり抜くと遠近法で翼の断面が見えてしまいますので、純粋に2D表示に戻すためにSliceしなければいけません。次いでThresholdを追加し、ScalarsをcellTypes、Lower Thresholdを0、Upper Thresholdを1.5、ColoringをUとします。ParaViewはcellTypesを表示することができますが、翼の部分のcellTypeが2となっていることが分かります。そこでThresholdを1.5以下とすることでcellType2の翼を表示から除外することができます。最終的にThreshold1だけを表示させると、翼がくり抜かれて流れがとても見やすくなります。このときの設定を以下に示します。では、ローターの回転を観察してみましょう。以下は1200rpmでローターを回して1秒が経過した時の様子です。流速分布流速グリフ表示 Scale Factor 0.004迎角はローターが領域左側で最大となるので空気は左側から流入して右側に流出するものと予想していましたが、シミュレーションによれば下側から多く流入して少し右上に流出する傾向にあることが分かりました。下流の流れは論文のFigure 16.に近い結果であり、CFD、レーザードップラー流速計による実測結果とも傾向が一致しています。ところで、この解析の2ステップ目、つまりシミュレーション開始後1mSecの状態を見ると、airfoil3とairfoil4の回りに大きな流速が生じています。これは1ステップ目のairfoil3とairfoil4は初期配置の時に迎角をつけておらず、たった1mSecで角度が30度変化したために計算上このような結果が出たものと思われます。これを防ぐにはスクリプトにtransformPoints -rotate-z <deg>を追記して初期配置の段階で翼に迎角をつけておけばスムーズに始動できます。しかしその後の角度データ全てにオフセットをかけることになり、シミュレーションで最も重要なパラメータである迎角を直感的に捉えることができなくなります。とはいえローターを1秒間も回せばその影響は消え去るようですので、初期配置の段階で迎角をつけないことにしました。連続的にステップを観察すると、下流側の翼が気流を乱しがちなのが分かります。論文によればコンロッドの取り付け位置であるdを変化させると上流側と下流側の最大迎角に差をつけることができます。下流側の翼がもう少し良い角度で気流に相対するよう、dあるいは翼弦長の50%としたPを調整すれば効率を改善できるかも知れません。改めて論文を読むと、最適なPの位置は翼弦長の35〜50%としか書いていませんので彼らも翼弦長50%がベストではないと考えているようです。なお、様々な研究によりサイクロローターが生み出す推力は翼の揚力によるもので、その方向は気流の方向と一致しないことが分かっています。以上で、サイクロローターが生み出す気流を可視化することができました。簡略化のため論文で言及されていた領域中央のφ15mmシャフトの存在は考慮しませんでしたが、このケースの手法は6枚ローター、独特なサイクロイド関数、有人機クラスの機体にも適用できると思います。翼と力の働く方向の関係を整理しやすいように翼の初期位置を決めたにも関わらず複雑な気流を生み出すことから、プロペラやジェットエンジンと違って力の働く方向は実験的に求める必要があるだろうなという印象を受けました。追記後日、解析領域を1m四方に拡げて計算しました。左側より1m/sの風を入力し、ローターを1200rpmで10秒間回転させました。気流は左側から流入するようになりましたが、右斜め上に流出する傾向は変わりません。ローターの上側は至って穏やかなようです。こんなケースファイルでも欲しい方がいらっしゃるようですので、GitHubに置いておきました。これはあくまで趣味であり、シミュレーションの正確性は追求していませんのでご了承ください。ケースファイルの設定はこの記事と異なり、inletから1m/sの風を入力しています。openfoam-cases/cyclorotor at main · kirinote/openfoam-casesContribute to kirinote/openfoam-cases development by creating an account on GitHub.github.com