きりまるのブログ -2ページ目

きりまるのブログ -2ページ目

ホームピグアメブロ

芸能人ブログ人気ブログ

前ページ
次ページ

01Mar
- 続・ベイパーコーンについて考えてみた
  前回の記事ではOpenFOAMでベイパーコーンの可視化に失敗しました。シミュレーションの前提条件に問題があった可能性があるため、もう少し条件を変えて同じように検証してみることにしました。チュートリアルの動作確認今回は遷音速領域のシミュレーションに適しているとされる圧縮性ソルバーを採用します。手始めにOpenFOAM11に付属しているチュートリアルprismを実行してみます。このチュートリアルは圧縮性流体用非定常ソルバーrhoPimpleFoamを使って三角柱周りの超音速流れを解析します。乱流モデルは標準k-εモデルです。解析空間は全長1mです。チュートリアルでの流速は650m/sおよそマッハ1.9ですが、このまま実行しても他のブログと同じ結果が出るだけです。折角ですから、遷音速領域335.28m/sまで減速して結果を見たいと思います。チュートリアルの保存場所から作業場所$FOAM_RUNにケースファイルをコピーします。$ cp -r /opt/openfoam11/tutorials/fluid/prism $FOAM_RUN$ cd $FOAM_RUN/prism/prism/0/Uを開き、ファイル内の数値650を全て335.28に修正し保存します。解析空間がどのようなブロックで構成されているかparaFoamで見てみます。$ paraFoam -block$ blockMesh$ rhoPimpleFoam$ paraFoamポスト処理メッシュ流速分布静圧分布温度分布三角柱の後方に圧力と温度の低い領域が円錐状に広がっています。流速に至っては、後方で470m/sおよそマッハ1.37と音速を超える状況となっています。遷音速で飛ぶ航空機では、機体周りの気流が部分的に音速を超える可能性があることが分かりました。モデルのスケール調整この検証の過程でSSAM-Gen5という素晴らしい航空機モデルを手に入れましたが、OpenFOAMで解析するためにはスケールを正しく調整しておく必要があります。これまではモデルをXSimに投入するさいに都度スケールダウンしていました。今回はチュートリアルをベースにケースファイルを編集しますので、モデルを使用する前にスケールを調整しておきます。gen5_full.stlをフォルダ$FOAM_RUNに配置し、surfaceConvertコマンドでスケールを1/1000に縮小します。$ cd $FOAM_RUN$ surfaceConvert -clean -scale 0.001 gen5_full.stl gen5_real.stl$ surfaceCheck gen5_real.stl生成されたgen5_real.stlは、OpenFOAMから1/1スケールのモデルとして認識されます。surfaceConvertコマンドを利用して生成されたSTLファイルをメモ帳で開くと、solidがpatch0という名称となっています。一方、XSimの[形状のインポート]-[形状の編集]タブでgen5_full.stlから縮小したSTLファイルは、solidがgen5_fullという名称となります。snappyHexMeshDictのregionsは、surfaceConvertを使用した時とXSimを使用した時とで名称が異なります。XSimで生成したsnappyHexMeshDictを流用する時は、regionsの定義をgen5_fullからpatch0に修正する必要があります。今後ケースファイルに組み込む航空機モデルにはgen5_real.stlを使用します。ケースファイルvpconeを作成先ほど流速Uを335.28m/sに変更したチュートリアルprismを、vpconeに名称変更します。/vpcone/constant配下にフォルダtriSurfaceを新規作成し、その中にgen5_real.stlを配置します。/vpcone/system/blockMeshDictを以下のように全面的に書き換えます。convertToMetersは1で、頂点verticesの座標がそのままメートルになります。prismの複雑なブロックメッシュを取り払い、代わりに全長50m、幅30m、高さ30mの簡素なブロックを作成します。ブロックの構造についてはこちらを参照してください。ベースメッシュは25cmのメッシュで、288万セルとなります。使わないprismWallを取り除き、代わりに解析空間の側面に上下と境界条件と同じsideWallを定義します。これで全体が壁で閉じたからっぽの風洞ができます。航空機モデルgen5_realのサーフェスは後ほど組み込みます。FoamFile{ format ascii; class dictionary; object blockMeshDict;}convertToMeters 1;vertices( (-10 -15 -15) // vertex number 0 ( 40 -15 -15) // vertex number 1 ( 40 15 -15) // vertex number 2 (-10 15 -15) // vertex number 3 (-10 -15 15) // vertex number 4 ( 40 -15 15) // vertex number 5 ( 40 15 15) // vertex number 6 (-10 15 15) // vertex number 7);blocks( hex (0 1 2 3 4 5 6 7) (200 120 120) simpleGrading (1 1 1));defaultPatch{ type empty;}boundary( inlet { type patch; faces ( (0 4 7 3) ); } outlet { type patch; faces ( (1 2 6 5) ); } bottomWall { type patch; faces ( (0 3 2 1) ); } topWall { type patch; faces ( (4 5 6 7) ); } sideWall { type patch; faces ( (0 1 5 4) (2 3 7 6) ); });sideWall用の定義を追加/vpcone/0配下にはT，U，alphat，epsilon，k，nut，pと７つのファイルでそれぞれ定数が定義されています。解析空間の側面にsideWallを追加したので、ここにも定数を全て定義します。また、prismWallはgen5_realに書き換えます。これらのファイル全ては手順が同じですので、ファイルTを例にとって説明します。なお、元のファイルTは温度が300Kすなわち摂氏26.85℃と定義されていますが、国際標準大気表より海面高度0mで摂氏15℃としたいのでuniform 300のところを全て288.15に置き換えます。・・・中略・・・ topWall { type inletOutlet; inletValue uniform 288.15; value uniform 288.15; } sideWall // topWallからコピー＆ペースト { type inletOutlet; inletValue uniform 288.15; value uniform 288.15; } gen5_real // prismWallから名称変更 { type zeroGradient; }・・・中略・・・見てのとおり、追加したsideWallはtopWallと同じ定義を設定しています。また、航空機モデルgen5_realのサーフェスはprismWallのゼロ勾配を引き継いでいます。残り６つのファイルもこれと同じ要領で定義しなおしてください。航空機モデルの組み込み続いて、解析空間に航空機モデルを組み込みます。/vpcone/system/surfaceFeaturesDictを新規作成し、以下のように記述します。FoamFile{ version 2.0; format ascii; class dictionary; location "system"; object surfaceFeaturesDict;}surfaces( "gen5_real.stl");includedAngle 165.0subsetFeatures{ nonManifoldEdges no; openEdges yes;}writeObj yes;このファイルにより、gen5_real.stlを読み込んで特徴線の抽出を行うことができます。STLモデルへのメッシュ適合続いて/vpcone/system/snappyHexMeshDictを新規作成し、以下のように記述します。locationInMeshでは、どこにメッシュを生成するかを指定するため解析空間の内側かつ航空機モデルの外側の座標（-5, 0, 0）を指定しています。regionsがpatch0となっていることに留意してください。このpatch0は、航空機モデルのASCII形式STLデータをエディタで開くと行頭solidと行末endsolidの２か所に定義されています。ケースファイルを実行したとき、STLファイル名が合っているのにFATAL ERRORが出る場合は、この定義がsnappyHexMeshDictのregionsと一致していない可能性があります。FoamFile{ version 2.0; format ascii; class dictionary; location "system"; object snappyHexMeshDict;}castellatedMesh true;snap true;addLayers true;geometry{ gen5_real.stl { type triSurfaceMesh; name gen5_real; regions { patch0 { name gen5_real; } } }}castellatedMeshControls{ features ( { file "gen5_real.eMesh"; level 1; } ); refinementSurfaces { gen5_real { level (0 0); regions { patch0 { level (0 0); patchInfo { type wall; } } } } } refinementRegions {} locationInMesh (-5.0 0.0 0.0); maxLocalCells 100000000; maxGlobalCells 100000000; minRefinementCells 1; nCellsBetweenLevels 3; resolveFeatureAngle 30; allowFreeStandingZoneFaces false;}snapControls{ nSolveIter 30; nSmoothPatch 3; tolerance 4.0; nRelaxIter 5; nFeatureSnapIter 10;}addLayersControls{ layers { } relativeSizes true; expansionRatio 1.0; finalLayerThickness 0.3; minThickness 0.03; nGrow 1; featureAngle 60; nRelaxIter 5; nSmoothSurfaceNormals 1; nSmoothNormals 3; nSmoothThickness 10; maxFaceThicknessRatio 0.5; maxThicknessToMedialRatio 0.3; minMedianAxisAngle 130; nBufferCellsNoExtrude 0; nLayerIter 50; nRelaxedIter 20;}meshQualityControls{ maxNonOrtho 65; maxBoundarySkewness 20; maxInternalSkewness 4; maxConcave 80; minFlatness 0.5; minVol 1.00E-13; minTetQuality -1e30; minArea -1; minTwist 0.05; minDeterminant 0.001; minFaceWeight 0.05; minVolRatio 0.01; minTriangleTwist -1; nSmoothScale 4; errorReduction 0.75; relaxed { maxNonOrtho 180; }}debug 0;mergeTolerance 1e-5;以上のように大規模なファイルですが、色々なケースファイルを見ると、ファイルを分割するテクニックを利用しているものを見かけます。慣れて来たらそのようなテクニックを真似ると良いでしょう。並列処理への対応現代のパソコン向けＣＰＵは、当然のようにコアを複数個搭載しています。チュートリアルのprismは小規模であったため、１コアでも大して時間はかかりませんでした。ですがそれなりの結果が出れば調子に乗ってどんどん大規模な解析をしてしまうので、今のうちにマルチコア対応しておきます。以下の内容で/vpcone/system/decomposeParDictを新規作成してください。ここでは例として解析空間を４つのサブドメインに分割し、ひいては４コアを使用します。FoamFile{ version 2.0; format ascii; class dictionary; location "system"; object decomposeParDict;}numberOfSubdomains 4;method scotch;たったこれだけの記述でマルチコアＣＰＵによる並列処理が可能になります。今回は細かい設定の要らないscotchメソッドを採用しましたが、他にも色々な分割方法があります。解析時間の設定最後に解析時間と間隔などを決めます。/vpcone/system/controlDictを新規作成し、以下のように記述します。特に重要な設定箇所は解析時間endTime、解析間隔deltaTです。これらはメッシュ密度と相まって解析結果に大きく影響します。計算実行中にクーラン数が出力されますが、最大クーラン数が１を下回るようにパラメータを調整します。記録間隔writeIntervalは、細かく刻むほどアニメーションがなめらかになります。FoamFile{ format ascii; class dictionary; location "system"; object controlDict;}application foamRun;solver fluid;startFrom startTime;startTime 0;stopAt endTime;endTime 0.08;deltaT 5e-05;writeControl runTime;writeInterval 0.008;purgeWrite 0;writeFormat ascii;writePrecision 6;writeCompression off;timeFormat general;timePrecision 6;runTimeModifiable true;vpconeの流速はprismのおよそ半分です。同じように解析しようとするとendTimeはprismの倍の時間が必要で、全長も50倍あるのでprismの100倍の0.04秒必要です。ですがこの間に航空機が進む距離はたったの13.4mです。せめて航空機の全長20mは空気を通過させたいので、倍の時間をとって0.08秒とします。writeIntervalを0.008秒とすると、解析時間全体で10コマ記録されます。deltaTはとりえず100倍の5e-05として、最大クーラン数がどうなるか見てみます。解析の実行ケースファイルが全て整ったら、解析をはじめます。並列処理は効率を上げるためバックグラウンドプロセスで行い、経過はlogファイルに書き出します。$ cd $FOAM_RUN$ cd vpcone$ blockMesh$ surfaceFeatures$ decomposePar$ mpirun -np 4 snappyHexMesh -overwrite -parallel > log &$ reconstructPar$ rm -rf ./processor*$ decomposePar$ mpirun -np 4 foamRun -solver fluid -parallel > log &$ reconstructPar$ rm -rf ./processor*$ paraFoamプロセスを並列処理するたびに解析領域の分割と結合をします。snappyHexMeshを実行すると290万セルとなりました。ソルバーを実行すると最大クーラン数は0.7となりました。クーラン数は１未満が適正ですので、いい感じに収まったと思います。並列処理中のメモリ使用量はおよそ8GBでした。バックグラウンドで動いているソルバーを止めたいときは、書き込みの続くlogファイルを見て、controlDictのendTimeを現在計算中の時刻より少し先に書き換えると安全に終了させることができます。ソルバーの完走に膨大な時間を費やす前に、取り急ぎ結果をチラ見するのに便利です。計算をただちに終了したい場合は、mpirunを実行した時にコンソールに表示されるプロセスIDを指定してkillコマンドを発行します。[1] 123456$ kill 123456ポスト処理解析空間メッシュ流速分布静圧分布温度分布参考画像出典：Wikimedia Commons, U. S. Navy with the ID 110606-N-TU221-408静圧分布（３Ｄ）温度分布（３Ｄ）考察これらの図から、機体の凸部後方から低温低圧の領域が円錐状に広がっているのがわかります。キャノピー後方、エアインテーク下部で低温低圧となっている他、不思議なことに尾翼後方に大きな温度勾配があります。この低温低圧が凝結雲の発生する条件を満たすとすれば、前回記事で紹介した写真に見える雲の出現位置と形状ともによく合致します。機体後部には高温のジェット気流もありますから、急激な温度上昇によって凝結雲が消えます。このため一見ボリュームがあるように見える雲も、特に機体後方で生じているのは円錐の表面だけと考えられます。機体の速度は音速に満たない遷音速335.28m/sですが最大流速は384.51m/sおよそマッハ1.1となっており、機体周辺に流速が音速を超える超音速流が見られます。ベイパーコーン形成の原因の一つとして考えていた「音」は、結果として生じたとしても関与はしていないことが分かりました。ただこれは学術的な研究ではなく素人の考えに過ぎませんので、現象の説明にはあくまで専門家の見解が必要です。OpenFOAMでは無理かなと思ったこともありましたが、なんとか遷音速で飛ぶ航空機に現れるベイパーコーンを可視化することができました。まだ空力的に洗練されていなかった頃の戦闘機モデルで試すとより派手な円錐が見られるかも知れません。計算機の仕様筆者の計算機は1CPU Xeon E5-2695 v4 128GBです。CPUクーラーはNH-P1、GPUはRTX3050 KalmX 6GB、電源はPX-500を使用し、完全ファンレスで静かに稼働します。今の所GPUをソルバーの並列処理に使うことはできませんが、NVIDIA RTXシリーズはparaView NVIDIA IndeX™プラグインに対応しておりデータセットを高速に表示することができます。
14Feb
- べイパーコーンについて考えてみた
  先日、空力を研究されている方から航空機が生み出すベイパーコーンについての話題を提供して頂き、あれこれ考えるうちに解析してみたくなりました。とくに興味を持ったのは、ベイパーコーンの発生する位置とその形状についてです。ベイパーコーンとは遷音速で飛ぶ航空機に見られる現象で、以下の画像のような円錐状の雲が生じます。面白いのは、雲が機体全体を覆っているのではないという点です。出典：Wikimedia Commons, U. S. Navy with the ID 110606-N-TU221-408ベイパーコーンは空気中の水分が断熱膨張冷却によって凝結されて雲になったものです。機体周りの一部で気圧が低下し、湿度の高い環境で露点を下回ると空気中の水分が凝結して雲が発生します。この現象は湿度の高い空気の入ったフラスコから注射器で空気を一気に抜くとできるのと同じ原理によります。速度と関係なく見られる現象ですが、機体周りにプラントル・グロワートの特異点が現れた時に急激な気圧勾配が表れます。この特異点を表す公式にマッハ数Ｍが入っているために、音速が関係しているように見えるのです。詳しい解説は以下のWikipediaを参照してください。プラントル・グロワートの特異点 - Wikipediaw.wiki実験でベイパーコーンを発生させるにはどんな規模の風洞設備を持ってしても困難と思われますので、シミュレーションでどうにかできないか検討します。仮想空間に条件を作り出して雲の発生を観察できるような超絶便利なソフトを筆者は知らないのですが、OpenFOAMで円錐状の流速勾配／圧力勾配が再現できれば雲ができたことになると考えました。が、結論から言いますと可視化できませんでした！以下の手順は、それが判明する前から備忘録を兼ねて編集していたもので、消すのが勿体ないので記録に残します。航空機モデルの用意まずは音速で飛びそうな航空機モデルを探索します。戦闘機の精密な３Ｄモデルならネット上に多くありますが、これらはエアインテークの穴もモデリングされており、ここに流れ込む気流が全体の計算結果に影響する恐れがあります。従って、ジェット機の空力性能を正確に評価するには、結果を比較するために標準化されたモデルが必要です。この分野では長らくカナダ国立研究評議会が考案した”標準空力モデルSDM”が使われてきました（下の画像）。しかし現在では時代遅れのモデルとなっており、これでやってみようという気になれません。出典：VSP Airshow, Standard Dynamics Model, DSTO-TR-2600そこで今回はオーストラリアのニューカッスル大学とシドニー大学の研究チームが考案した”シドニー標準空力モデルSSAM”を採用します。これは最新鋭機の空力性能を評価するためのモデルで、その形状は米軍の第５世代戦闘機F-22ラプターを模したものとなっています。このモデルで空力性能を解析した研究論文はこちらで読むことができます。SSAM-Gen5と呼ばれるこのモデルはエアインテークが塞がれたソリッドモデルとなっており、エアインテークの存在が空力計算に悪影響を及ぼさないよう配慮されています。見た目も洗練されており、これならいい結果が得られそうな気がします。出典：Aerospace 2023, 10(9), 746,A Generic Model for Benchmark Aerodynamic Analysis of Fifth-Generation High-Performance Aircraftこの論文で使われた機体モデルはこちらから入手することができます。モデルデータはSTEP形式で収録されています。FreeCADでgen5_full.stepをインポートして、そのままSTL形式でgen5_full.stlにエクスポートすると良いでしょう。Blenderでインポートするには別途プラグインが必要です。モデルを実寸に修正する公開資料によれば、F-22実機の寸法は以下の通りです。　全長 62 ft 1 in= 18.92m　翼幅 44 ft 6 in= 13.56mgen5_full.stlの機体全長は1となっていますので、BlenderでモデルをインポートしてXYZスケールを18.92倍します。同名のSTLファイルにエクスポートすると、全長18.92となったモデルが上書きされます。この翼幅を計測すると13.499となり、実機の翼幅にかなり近い数値になることが分ります。モデル原点はノーズコーンの先端です。XSimの起動こちらのチュートリアルに従って解析領域を決定します。モデルの全長をおよそ20mとして、全長160m、全幅50m、全高50mの風洞にモデルを配置します。XSimのプロジェクト名はvpconeとします。形状のインポートモデルファイルgen5_full.stlをドロップします。例によって全長18.92kmというとんでもない大きさで認識されますので、[形状の編集]タブにて拡大縮小率を0.001します。[形状の追加]タブ-[タイプ]：直方体[形状の追加]タブ-[パラメーター]-[最小座標]：(-40, -25, -25)[形状の追加]タブ-[パラメーター]-[最大座標]：(120, 25, 25)直方体を追加すると以下のようになります。メッシュモデルを覆うように再分割領域を追加します。目標ベースメッシュ数は空間分解能と計算時間のトレードオフになりますので、どちらも許容できる範囲に調整してください。[体積メッシュ設定]-[目標ベースメッシュ数]：80000[体積メッシュ設定]-[再分割設定]-[レベル間のメッシュ数]：3[体積メッシュ設定]-[再分割設定]-[範囲タイプ]：直方体[体積メッシュ設定]-[再分割設定]-[最小座標]：(-10, -10, -10)[体積メッシュ設定]-[再分割設定]-[最大座標]：(50, 10, 10)[体積メッシュ設定]-[再分割設定]-[再分割レベル]：2[体積メッシュ設定]-[再分割設定]-[範囲タイプ]：表面（領域）[体積メッシュ設定]-[再分割設定]-[領域：gen5_full[体積メッシュ設定]-[再分割設定]-[再分割レベル]：3[レイヤーメッシュ設定]-[厚み比率]：0.3[レイヤーメッシュ設定]-[領域]：gen5_full[レイヤーメッシュ設定]-[レイヤー層数]：3基本設定非定常解析を行います。[終了時刻]：1（秒）[時間間隔]-[クーラン数による自動設定]：オン[時間間隔]-[初期時間間隔]：0.0001（秒）[時間間隔]-[最大クーラン数]：0.9[解析タイプ]：乱流[乱流モデル]：Realizable k-ε モデル物性[物性名]：Air初期条件こちらにて海面高度75ft、750mphで飛行するF/A-18Cでベイパーコーンが観測されたとのことですので、流速は335.28m/sとします。海面高度はそれほど高くなく、空気の物性は標準値とほとんど変わらないので0mとします。海面高度0mにおける音速は340.29m/sです。[物理量]：速度[値]：(335.28, 0, 0)流れ境界条件流入口、流出口、モデル表面のみ条件を設定します。条件が設定されていない領域は断熱のすべり壁になります。[領域]：XMin[タイプ]：流速指定[流れ境界条件パラメーター]-[流速]：335.28[流れ境界条件パラメーター]-[流入温度]：288.15[領域]：XMax[タイプ]：自然流入出[領域]：gen5_full[タイプ]：静止壁計算設定[並列数]：お使いのコンピュータのCPUコア数出力代表面積は揚力方向を指定すると自動計算されます。[タイプ]：指定時刻ごと[間隔]：0.1（秒）[全体]タブ：残差，y+値[領域]タブ-[領域]：gen5_full[領域]タブ-[タイプ]：力係数[領域]タブ-[代表速度]：335.28[領域]タブ-[代表長さ]：18.92[領域]タブ-[代表面積]：109.59858633546901（揚力方向）エクスポート[フォーマット]：OpenFOAM 11エクスポートをクリックするとvpcone.zipがダウンロードされますので、これを展開して得られたフォルダvpconeを$FOAM_RUNに配置します。解析実行以下の手順で解析を実行します。$ cd $FOAM_RUN$ cd vpcone$ chmod 755 ./Allclean$ chmod 755 ./Allrun$ ./AllrunparaFoamの起動解析が終了したら、paraFoamを起動します。$ ./paraFoam結果の確認まずはモデル断面のメッシュを確認します。あまり綺麗ではありませんが、そのままです。なぜならば予想と違う結果が出たことでメッシュにこだわる理由を無くしたからです。流速のカラーマップは300m/s〜370m/sの範囲としています。メッシュ流速分布流速分布（3D Glyphs）静圧分布流線表示遷音速の気流中にモデルを置いたことでエッジの効いた形状部分に変化が見られますが、期待した円錐状の”何か”が見当たりません。流れもほぼX軸方向で整っており、とりたてて大きな乱流がモデル周りに生じているように見えません。あえて気になるところと言えば、機体形状の凹凸部分で静圧が低下しているというごく普通の結果だけです。これより、ベイパーコーンは気流とあまり関係なさそうであることが分かりました。考察のようなものOpenFOAMを使えばすべてを明らかにできると思いこんでいましたが、この考えは改めねばなりません。思いついたのは、この記事の冒頭で排除した「音（圧力変動）」が関係しているかも知れないことです。今回のシミュレーションでは音波の影響を考慮していませんので、当然、結果にも反映されません。ここで音の性質を確認しておきます。以下の図は音源が移動した場合に音波がどのように拡散するかを示したものです。音源が静止している場合は全周一様に音が伝わりますが、音源が移動するとその方向に音波が偏り、音速では音源から音が逃れることができずに大きな衝撃波が生じます。超音速（Supersonic）に入るともはや音が音源に追従できず、衝撃が円錐状に広がるようになります。円錐の鋭さは速度に依存し、音速の２倍で移動すると角度θは30°となります。ただし、記事冒頭のベイパーコーンの画像は音速未満の速度で観測されたものですから、雲の形が似ていようとも超音速の計算を当てはめることはできません。しかしながら、遷音速（Transonic）においても音の偏りによって衝撃波は生じていて、急激な気圧変化はすでに始まっているのかも知れません。円錐状の凝結雲がキャノピーなど機体凸部の後方に見えるのは、衝撃波が引き起こした気圧低下に機体横断面積の減少に伴う断熱膨張が重畳され、部分的に冷却効果が強まったためではないかと推測します。以上、ベイパーコーンに興味を持ったために余計なことをしてしまいました。この結論はあくまで推測ですので、正確な情報をお持ちの方がいらっしゃいましたらぜひコメント、ご連絡いただけたらと思います。追記もう少し調べてみたところ、遷音速領域の解析には圧縮性ソルバーsonicFoamやrhoCentralFoam/shockFluid等が適しているとの情報がありました。この手法を今後取り入れてみようと思います。
18Dec
- Ubuntu24.04のスワップ領域を設定する
  Ubuntu24.04をインストールすると、起動用ドライブに数GByteのスワップ領域swap.imgが作成されます。OpenFOAMで重たい解析を行うとメインメモリ領域を使い果たし、このスワップ領域を使うようになりました。そうそう買い換えることのない大容量SSDを頻繁に書き換えると寿命を著しく縮めてしまうので、スワップ領域を他のドライブに割り当てることにしました。今回、このパソコンにM.2接続のNVMeドライブを増設し、スワップ領域をそちらに移して起動用ドライブを保護しようと思います。増設したNVMeドライブを認識しているか確認するcd /sudo fdisk -lパーティションを作成するsudo fdisk /dev/nvme0n1fdiskの質問には以下の順番で回答します。署名を削除しますか？と訊かれたらyと答えます。　d → n →1→Enter→Enter →（y）→ p →w現行のスワップ領域を解放してスワップファイルを削除するsudo swapoff /swap.imgsudo rm /swap.imgNVMeドライブ全体をスワップ領域として確保するsudo mkswap /dev/nvme0n1p1sudo swapon/dev/nvme0n1p1sudo swapon -sfree -hスワップ領域を永続化させるテキストエディタVimを設定ファイルの編集に用います。使いかたは省略します。sudo apt install vimsudo vim /etc/fstabファイルfstabの最終行を以下のように変更して保存します。/dev/nvme0n1p1 swap swap defaults 0 0仕上げに再起動するreboot以上でスワップ領域を別ディスクに退避しつつ容量を拡大することができました。
28Nov
- プロペラを解析する
  前回APCプロペラをモデリングしましたので、7000rpmで回転した時の静的推力を解析してみます。OpenVSPで前回保存したプロペラモデルを読み込みます。[Analysis]-[Aero]-[VSPAERO...]をクリックすると、設定画面が現れます。OverviewタブVLM法で解析します。[Case Setup]-[Voltex Lattice (VLM)]：ONプロペラが地面に固定された状態と同じですので、Flow Conditionの設定による影響はないようです。[Flow Condition]-[Alpha Start]：0.000[Flow Condition]-[BetaStart]：0.000[Flow Condition]-[MachStart]：0.000[Flow Condition]-[ReCref Start]：500000Advancedタブお使いのパソコンのCPUコア数以下に設定してください。コア数が増えるほど解析時間が短縮されます。[Num CPU]：CPUコア数ビューワーで表示する際の航跡の密度を指定します。解析にはおそらく影響しません。[Wake Nodes]：64Actuator Diskは定常解析、Rotating Bladesは非定常解析です。通常はRotating Bladesを選択します。[Propeller Representation]-[Rotating Blades]：ONVinfは自由流速です。プロペラは回転するだけでその場から動きませんので0m/sとなります。[Advanced Flow Conditions]-[Vinf]：0.00Vrefはプロペラ先端の移動速度です。毎秒回転数[rps]にプロペラの円周[m]を掛けます。Vref = 7000rpm／60sec・0.5588m・π = 204.81[m/s][Advanced Flow Conditions]-[VRef]：204.81MachRefはプロペラ先端のマッハ数で、流速[m/s]と音速340.29[m/s]から計算で求めます。別途計算サイトを使うと便利です。[Advanced Flow Conditions]-[MachRef]：0.602Rhoは空気の密度で、海面高度0kmにおいて1.225[kg/m3]です。[Advanced Flow Conditions]-[Rho]：1.225Propellerタブ非定常解析ですので、いつまで解析するか決めます。プロペラが10回転するまで解析します。時間にするとおよそ85.7mSecです。[Time]-[Auto Time Step]：ON[Time]-[Num Revs]：10回転数がエンジン始動直後から7000rpmで安定しているとします。[Propellers]-[RPM]：7000.000ここまで設定できたらダイアログ左下の[Launch Solver]をクリックして解析をはじめます。解析が終わると結果マネージャが自動的に開きます。色々グラフが描かれますが、Unsteadyタブを開き、[Y-Data]からThrustを選択すると全時間分の推力グラフが表示されます。横軸は85.7mSecになっています。時折大きなノイズが入るものの、推力はおよそ100[N]辺りになっています。正確な数値は別のファイルに書き出されています。OpenVSPフォルダに、22x10-PERF_DegenGeom.rotor.1というファイルがありますのでこれをメモ帳などで開きます。行を下っていくと、次の項目が見つかります。左から三つ目の数値がThrust（推力）です。Average over last full revolution: 0.559 7000.001 111.645 -1.082 112.727 4627.911APCプロペラ22"×10"の7000rpmでの静止推力は111.645[N]=11.384[kgf]と解析されました。ここで、静止推力Fの一般的な計算式は以下の通りです（こちらで紹介したGabriel Staples氏考案の２枚羽根校正済み静止推力の方程式とは異なります）。F：静止推力[gf]D：プロペラの直径[inch]P：プロペラのピッチ[inch]N：プロペラの毎分回転数[rpm]K：比例定数。ここでは22とするこの式によれば、静止推力Fは11478[gf]=11.478[kgf]となります。もう少し色々なサイズのプロペラで検証が必要かも知れませんが、VSPAEROの解析で計算式に近い数字が得られたものと思います。
26Nov
- プロペラをモデリングする
  模型飛行機や無人機で使用しているプロペラを評価するために、正確な3Dモデルが欲しくなることがあります。ネットには出来合いの3Dモデルもありますが、手元にあるプロペラと形状が異なっていては評価する意味がありません。プロペラは複雑な形状をしているため、各部の寸法や角度を測って3DCADやOpenVSPに入力していくのは困難です。プロペラ形状の技術的な解説はネット上にありますが、ノウハウの塊なのか商品となっているプロペラについての情報は見つかりません。そんな中でも業界大手のAPCプロペラは形状データを公開しており、うまくやればOpenVPSで読めるフォーマットに変換することができます。仕様が明らかであるという点を買って、APCブランドを選ぶ人がいても不思議はないと思います。JuliaのインストールAPCプロペラが配布している形状データをOpenVPSで読める形式にコンバートするユーティリティAPC2BEMは、新興のプログラミング言語Juliaで作られています。こちらからJuliaのインストーラを入手してください。Windows11で何も変更しなければ以下の場所にインストールされます。C:\Users\ユーザー名\AppData\Local\Programs\Julia-1.11.1Juliaコンソールを開き、キーボード右側の右カッコ "]"キーを押してパッケージ管理モードに入ります。２つの追加パッケージをインストールします。(@v.1.11) pkg>add ArgParse(@v.1.11) pkg>add InterpolationsインストールしたらBackSpaceキーを押して元のモードに戻り、Juliaコンソールを閉じます。これでコンバートの準備ができました。APC2BEMの入手こちらよりコンバートプログラムAPC2BEMをダウンロードし、解凍して取り出したファイルapc2bem.jlを先ほどできたフォルダJulia-1.11.1\binに配置します。形状データの入手APCプロペラのホームページからプロペラ形状のデータを入手します。[TECHNICAL INFO]-[MORE]-[DOWNLOADS]-[Propeller Geometry Data]からPE0-FILES_WEB-202410.zipxをダウンロードします。拡張子zipxはWindowsで解凍できない圧縮形式なので、7-Zipを使って圧縮ファイルを解凍します。ここでは、22"×10"のプロペラ形状をコンバートします。ファイル名は22x10-PERF.PE0です。この.PE0ファイルをフォルダJulia-1.11.1\binに配置します。結果、\bin配下にこれらのファイルが揃います。形状データのコンバートでは早速BEMファイルにコンバートします。このとき、形状寸法をインチからメートルに変換します。> cdC:\Users\ユーザー名\AppData\Local\Programs\Julia-1.11.1> julia apc2bem.jl -s 0.254 22x10-PERF.PE0 22x10-PERF.BEM以上で同じフォルダに.BEMファイルが生成されました。このファイルをOpenVSPフォルダに配置します。形状データを開くOpenVPSを立ち上げ、[File]-[Import...]-[Blade Element (.bem)]より22x10-PERF.BEMを開きます。Geom BrowserからPropGeomを選択すると、Propeller: PropGeomが開きます。タブ：Design[Design]-[Construct X/C]：0.000[Design]-[Feather Axis]：0.125これらのパラメータは、APCプロペラの形状ファイルであれば全て同じです。ここまで編集できたら念のためファイルを保存しておきます。ファイル名は22x10-PERF.vsp3とします。この時点ではプロペラのスケールが間違ったままですが、後ほど修正します。翼型断面の確認さて、ここで一旦元ファイル.PE0の内容を確認します。APCプロペラが提供しているのは、プロペラの外形だけで翼型の断面形状は含まれていません。そこで、プロペラの断面にどんな翼型を採用しているのかの確認します。メモ帳で22x10-PERF.PE0ファイルを開き、最後の一文を確認すると以下のように書いてあります。----- AIRFOIL SECTIONS -----AIRFOIL1: 4.00, CLARK-Y (Transition Start, Airfoil 1)AIRFOIL2: 6.00, CLARK-Y (Transition End, Airfoil 2)AIRFOILS ARE SCALED BASED ON THICKNESS RATIO IN TABLE ABOVE.これよりプロペラの断面は全体がクラークＹであることが判りましたので、おなじみairfoiltools.comよりCLARK Y AIRFOILの翼型データを入手します。画面中央のDat fileにあるデータを全て選択し、これをメモ帳に貼り付けてclarky.af（UTF-8 BOM無し）として保存します。続いて、clarky.afを編集します。ファイル冒頭を以下のように書き換えます。GEOM AIRFOIL FILECLARK Y0 Sym Flag (0 - No, 1 - Yes)61 Num Pnts Upper61 Num Pnts Lower0.0000000 0.00000000.0005000 0.00233900.0010000 0.00372710.0020000 0.0058025・・・以下略・・・続いて、２カ所あるX座標1.0000000の行についてほぼ0である（が、ゼロではない）Y座標を0.0000000に書き換えます。[66行目] 1.0000000 0.0000000[128行目]1.0000000 -.0000000今回の場合は以上の編集で正しく読めるようになります。元データの形式は様々ですので、翼型データによっては追加の編集が必要になる場合があります。ファイルフォーマットの詳細はこちらをご確認ください。翼型断面の読み込みOpenVPSを起動し、プロペラファイル22x10-PERF.vsp3を開きます。[Propeller: PropGeom]の[XSec]タブで、[Type]-[Choose Type:]にAF_FILEを選択するとファイルを読み込めるようになります。ここに先ほど作成したclarky.afを登録します。翼根側の翼型クラークＹを記憶するため、[Cross Section]-[Copy]をクリックします。[Cross Section]を０→１に切り替えて[Paste]をクリックし、翼端側にもクラークＹを割り当てます。すると、プロペラ断面全体がCLARK Yの形に整います。さて、プロペラの先端を見るとポリゴンが予定した形に追従できていません。テッセレーションを強化して先端形状を予定の形に整えます。タブ：Gen[Tessellation]-[Num_U]：36タブ：More[Tessellation Control]-[Tip Clustering]：0.1これでプロペラの形が整いました。スケールの調整ここで、[Design]-[Diameter]を確認します。ファイル形式をコンバートした際にスケールを調整したので直径22inch=0.5588mになっているはずですが、なぜか直径5.588mとなっています。とりあえず、スケールを１／１０します。タブ：XForm[Scale Factor]-[Scale]：0.1 → Accept[Design]タブに戻り、[Diameter]が0.5588mとなっていることを確認してください。一見してハブ側の隙間、いわゆるシャンク部に違和感を覚えるかも知れませんが、22x10-PERF.PE0で仕様を確認すると以下のように記述されています。RADIUS: 11.00 PROPELLER RADIUS (IN)HUBTRA: 3.13 HUB TRANSITION (IN)BLADES: 2 NUMBER OF BLADESハブからの距離が3.13inch=0.0795mとあり、この部分は空力的に意味がなく実機では大部分がスピナーコーンに覆われます。片羽分の距離なので、直径は0.0795×2=0.159mであり実測値と合っています。以上でプロペラのモデルが完成しました。筆者は元ファイル名を踏襲し、22x10-PERF.vsp3としてモデルを保存しました。
19Nov
- ステルス性能を解析する
  無人機を製作したら、それがどれぐらいレーダーに映るか知りたくなります。ここでは、ブラジル陸軍工兵学校の学生さんが米海軍大学院教授David Jenn氏の開発したレーダー反射断面積解析プログラムPOfacetsを参考に開発したOpen RCSを使用してみます。このソフトは高周波近似解法の一種である物理光学近似（PO）と幾何光学近似（GO）を使用してRCS計算を行います。Python3で動作します。Python3のインストールこちらからPython3の最新版をダウンロードしてお使いのＰＣにインストールします。現時点での最新バージョンはPython 3.13.0です。Add python.exe to PATHにチェックを入れてどこでもPythonを実行できるようにしておきます。Install Nowでインストールが完了します。Open RCSのインストールこちらからopen-rcs-main.zipをダウンロードして解凍します。フォルダ名をopen-rcsとして作業用フォルダに移動します。Windows11の場合は、C:\Users\ユーザー名\に配置すると良いでしょう。open-rcsフォルダを開き、そこでresultsという名称のフォルダを新規作成します。現バージョンではこのフォルダがないと動作しません。Open RCSの起動コマンドプロンプトを立ち上げ、Open RCSを起動します。初回起動時のみ、動作に必要なライブラリをインストールします。> cd c:\Users\ユーザー名\open-rcs> pip install -r requirements.txt> python openrcs.pyOpen RCSの画面が表示されます。ブラジルで開発されたソフトなので、ユーザーインターフェースはポルトガル語となっています。モノスタティック測定モードとバイスタティック測定モードを切り替えることができます。モノスタティック測定とは、送受信レーダー１基で探知する測定方法です。バイスタティック測定とは、別の場所に配置した送信用、受信用のレーダー２基が連携して探知する測定方法です。Monoestático：モノスタティック測定モードBiestático：バイスタティック測定モードFrequência (GHz)：レーダー周波数Distância (m)：相関距離Desvio Padrão (m)：標準偏差Phi Inicial (º) ～Phi Final (º)：Φ角度範囲Passo Phi (º)：Φ角度ステップTeta Inicial (º) ～Teta Final (º)：θ角度範囲Passo Teta (º)：θ角度ステップPolarização：入射波の偏波。平行偏波TM-zまたは垂直偏波TE-zResistividade：材質の電気抵抗率。特定の材料または完全導体材料Upload Modelo (.stl)：STLモデルのアップロード画面下側のGerar Resultados do Input Fileボタンを押すと、input_filesフォルダ内の.datファイルが読み込まれて何かしら解析結果が表示されます。毎度これらのパラメータを入力するのが面倒なときは、.datファイルを編集することで時間を節約できます。解析結果の妥当性を確認するまずはOpen RCSを使って意味のある結果が得られるかを確認します。RCSの評価においては、理論値との比較を容易にするため、完全導体材料PEC（Perfect Electric Conductor）を定義します。完全導体球のモノスタティックRCSの理論値は、球の半径a=0.05mと置くと以下の式で与えられます。RCS = πa2= π(0.05)2= 0.007854[m2]これをdBsm単位に変換すると以下の数値になります。RCS = 10log10(0.007854) ≈ -21.049[dBsm]理想的な条件では周波数2.998GHz（波長0.1m）のとき、すべての角度において-21.049dBsmが得られます[1][2]。解析に利用する半径0.05mの球体を用意します。FreeCADで生成できるSTLモデルの精度には限界があるため、商用CADを使用して高精度なSTLモデルを作成し、解析を行いました。画面右下のDownload Arquivo de Dados（データファイルのダウンロード）ボタンをクリックすると、resultsフォルダに数値データファイルを出力します。この.datファイルのRCS Theta (dBsm)を確認すると、概ね-21.043dBsm近辺と理論値に近い値が得られました。続いて、NASAがレーダーターゲットのベンチマークとして考案した円錐球モデルの解析をします[3]。このモデル断面の各部寸法（単位：mm）を以下に示します。こちらも商用CADにて高精度なSTLモデルを作成しました。垂直偏波9GHzにおける円錐球のモノスタティックRCSを解析すると、以下のような結果が得られました。既知の測定結果[3]とよく一致しているように思えます。この結果を見てもまだレーダー波をどの方向からどの方向に向かって走査したのかは分りません。そこで、θおよびΦのパラメータと照射方向の関係性を把握するため３つのパターンでデータを取りました。　　以上より、θおよびΦは次の角度を示すことがわかりました。RCSを平方メートルで語ることが多々ありますので、以下の変換式を覚えておくと役立ちます。0dBsm=1m2です。1m2のPEC板に垂直偏波1GHzのレーダー波を照射した場合です。平板のRCSFPは縦a[m]、横b[m]、波長λ[m]として以下の計算式で求めます。RCSFP = 4πa2b2 ／ λ2 = 4π ／ (0.3)2 = 139.626[m2] = 21.449[dBsm]正面での解析結果は21.449dBsmであり、理論値と一致しました。最後に、NASAアーモンドと呼ばれるベンチマークターゲット[3]に垂直偏波7GHzのレーダー波を照射した場合です。NASAアーモンドは低RCS形状の代表的なモデルのひとつです。数値データが見つからないため、MATLABのヘルプページ[4]で紹介されているグラフと見比べます。物理光学/幾何光学近似法（PO/GO）の計算結果は正しそうですが、重い計算負荷と引き換えに高精度な結果が得られるモーメント法（MoM）の精度には及ばないようです。高精度な結果を求めるには、モーメント法など厳密解法に対応したソフトウェアが必要だということが分かりました。田中芳樹氏のＳＦ小説『銀河英雄伝説』では、ステルス技術の進歩で艦艇をレーダーで捕捉できなくなった遠い未来が描かれています。物語に登場する銀河帝国軍の最新鋭戦艦ブリュンヒルト（加藤直之氏デザイン）は、NASAアーモンドのような流線型をしています。これはレーダーRCSの観点からとても合理的なデザインだと言えるでしょう。モデルを解析するなかなか使えそうだという手応えがありましたので、Open RCSに付属しているモデルを解析してみます。旅客機や貨物機などの民間機は、用途に応じた形を存分に追求した機体形状なので大抵はレーダーにがっつり映りますし、そうでないと困ります。レーダー反射断面積を気にする航空機といえばやはり軍用機やミサイルです。Open RCSのモデルにもそのような傾向があります。同梱のF-35モデルを解析してみます。機体を上から見て後方から反時計周りに全周を水平走査します。グラフ中央が機体の正面になります。この機体のRCS解析結果はネットで検索するとよく出てきます。製造元がRCSを公表するとは思えないので、これらはどこかの研究機関が勝手に調べたものでしょう。世界最強のステルス機であるがゆえに、諸外国から徹底的にマークされていることが窺えます。得られる結果は機体形状と仮定の表面材質だけで計算したもので、その他のテクノロジーは反映しておらずあくまで参考程度です。一説にはここから電波吸収材RAM（Radar Absorbent Material）によって-10dBsmほど押し下げられ、実機のステルス性能は対Xバンドレーダーで-30dBsm程度になるとされています。ネットで見つけたB-2の無料モデルでも試してみます。全幅52m、全長21mの大型機でありながら、電波吸収材も特殊塗料もない状態ですでに平均-3dBsm=0.5m2のRCSしかありません。ズムウォルト級駆逐艦、水上艦艇もRCS解析の対象です。対するSバンドレーダーは海面反射と吸収の影響が少ない水平偏波です。喫水線下のモデルは削除しています。このブログで作成した無人機のSTLモデル、myUAV.stlもそのまま解析に使えます。以上、色々なモデルでRCSを調べてみました。RCSを考慮した航空機には独特の機能美が感じられます。隠しかたを書いたのなら、見つけかたについても触れておく必要があります。レーダーの探知能力を評価するために利用されるレーダー方程式がこちらです。Rmax：最大探知距離[m]P：送信パルス電力[W]G：アンテナの絶対利得[dB]λ：波長[m]σ：物体のレーダー反射断面積RCS[m2]Smin：最小受信感度MDS(Minimum Detectable Signal)[W]とりあえずこの式を使ってみます。以下のような仕様のＬバンド対空レーダー装置を仮定します。送信尖頭電力：5kW送受アンテナ利得：24dBiレーダー周波数：1.5GHz（波長0.2m）最小受信感度：-110dBm関数電卓がないと計算が難しいのですが、こちらの計算サイトにデータを入力すると簡単に求まります。たとえばRCSが-3dBsmの航空機は、この対空レーダーから23.7kmの距離で探知されます。レーダーに軍事技術のようなイメージを持つ方もいらっしゃると思いますが、航空管制や気象観測、近年は乗用車にもミリ波レーダーが搭載されるなどレーダー技術は身近な存在になりつつあります。様々な解像度の球体モデルを作成して比較した結果、良好な解析結果を得るにはできるかぎり高解像度のSTLモデルを作成するのがコツだとわかりました。機体の正確な形状を知っているのは、機体の設計者だけです。模型飛行機やドローンを作った人なら、その過程で使用したCADモデルが手元にあると思います。いまお持ちの機体のステルス性能を検討できる希な立場にあるので、ぜひ一度解析してみると良いでしょう。引用元および参考文献[1] 電子情報通信学会エレクトロニクスシミュレーション研究会http://www.ieice.org/es/est/activities/canonical_problems/[2] 「高周波回路設計におけるシミュレータ活用の勘所 Ⅰ, Ⅱ」Microwave Workshops & Exhibition (MWE) マイクロウェーブワークショップダイジェスト, WS08-10, pp.83-88, 2012年11月29日.[3] H.T. G. Wang, M. L. Sanders, A. C. Woo, and M. J. Schuh."Radar Cross Section Measurement Data, Electromagnetic Code Consortium Benchmark Targets".NWC TM 6985, May 1991.[4] https://jp.mathworks.com/help/antenna/ug/radar-cross-section-benchmarking.html
09Nov
- 無人機の主翼を解析する
  前々回モデリングした無人機をOpenFOAMで解析してみます。ウィングレットは無しとします。胴体や尾翼はほとんど揚力を生まないので、主翼だけを取り出して力係数を計算します。機体全体を仮想風洞に置いて解析したいところですが、その条件で起きる現象には今のところ用がないので、地味な解析になってしまいます。主翼モデルの準備OpenVSPでモデリングしたデータを読み込み、Geom Browserで主翼以外すべてをNoShow（非表示）にします。モデルの主翼は胴体に深く食い込まないよう機軸線から翼幅方向に0.1mオフセットした状態でした。このまま解析すると翼根側の隙間で意味のない気流の乱れが生じますので、一時的にオフセットを取り除いて翼を密着させます。[WingGeom]-[XForm]-[Transforms]-[YLoc]：0主翼の表面をなめらかに仕上げるため、テッセレーションの値を最大まで上げます。[WingGeom]-[Gen]-[Tessellation]：133[File]-[Export...]でファイルフォーマットにStereolith (.stl)を選択します。STL Optionsは全てOFFのままOKをクリックし、ファイル名をmyWing.stlとしてエクスポートします。FreeCADを起動してmyWing.stlを開き、視点をFrontに切り替えます。このときモデルの全長は2.3mmとなっていますが、OpenFOAMでは2.3mと認識されるためスケールはそのままで作業を進めます。テッセレーションの値を最大にしたことで、NACA4412の翼断面をきれいに再現できました。なお、この表面仕上げがシミュレーションの精度を左右するわけではありません。解析の前処理にてsnappyHexMeshが適切にメッシュを生成できるよう精密なモデルを渡すことが目的で、多段階のメッシングプロセスを経てSTLモデルに近似したポリゴンに変換されます。作業結果を確認するため、[View]-[Toggle axis cross]でXYZ軸原点を表示しておきます。Combo ViewのModelリストからMyUAVを選択し、[Property]-[Base]-[Placement]右側の三点リーダー[・・・]をクリックします。モデルを滑走路の上空2mに置いてシミュレーションしますので、モデルをZ軸方向に2mm移動させます。機体のピッチを変えたい時はこの時に傾けます。今回は機体ピッチ0度、つまり水平飛行とします。[Placement]-[Translation]-[Z]：2視点キューブ、XYZ原点との位置関係は以下のようになります。これで主翼モデルの準備が整いました。ファイル名をmyWing_2m.stlとして保存します。caseファイルの編集caseファイルには以前作成したmotorBikeSteadyの設定を流用します。motorBikeSteady/system/forceCoeffsファイルを開き、パラメータを以下のように設定します。対気速度27.78m/s、空力平均翼弦長0.36m、翼面積1.96m2です。magUInf 27.78;lRef 0.36; // Mean Aerodynamic ChordAref 1.96; // Estimated0/include/initialConditionsファイルを開き、以下のように編集します。flowVelocity (27.780 0);constant/geometryフォルダにmyWing_2m.stlをコピー＆ペーストします。端末を立ち上げ、STLファイルをOBJファイルに変換します。$ cd $FOAM_RUN/motorBikeSteady/constant/geometry$ surfaceConvert myWing_2m.stl motorBike.obj -clean$ rm myWing_2m.stl解析する全ての準備が整ったら、解析をはじめます。$ cd $FOAM_RUN$ ./Allclean$ ./Allrun結果をparaFoamで確認します。仮想風洞内に置いた翼モデルは以下の形状です。空力平均翼弦MAC位置、翼根から翼端方向に1.266m位置でのモデル断面です。３層のレイヤーメッシュで覆われています。高さ2mで流線表示すると以下のような表示になりました。素人目に翼端渦はそう激しくないように見えます。併せてVSPAEROでウィングレットが力係数にどのように影響するか調べましたが、小手先の工夫では特性の改善に結びつきそうにないためウイングレットは不要と判断しました。実際に翼を作って風洞実験をするよりシミュレーションを繰り返すほうがはるかに楽なので、翼端渦の影響が気になるようならいくつかモデルを作って様子を見ても良いかと思います。motorBikeSteady/postProcessing/forceCoeffs1/0/forceCoeffs.datを見てみます。最終行Time:500を見ると以下の数値が出ています。Cl：3.974416e-01Cd：2.499756e-02Cm：-2.157143e-01機体ピッチ0度で揚力係数0.397、揚抗比15.9と出ました。
05Nov
- 主翼にウィングレットを追加する
  前回作成した無人機の主翼にウィングレットを追加して、もう少しかっこよくしてみます。単に形状を整えることが目的で、空力的な検討はしません。OpenVSPで作成したモデルを表示し、[View]-[Left Iso]で作業しやすい視点にします。主翼以外のジオメトリを全て選択し、[Selection]-[NoShow]で主翼だけを表示させます。以降、主翼を編集します。現時点で、主翼にはひとつのセクションしかありません。まずはこれを主翼、トランジション部、ウィングレットの３つのセクションに分割します。タブ：Sect[Wing Section]-[Split]　// 分割[Wing Section]：2[Wing Section]-[Split]　// 分割[Wing Section]：1[Section Planform]-[Span]：2.8[Section Planform]-[Tip C]：0.25これで、主翼の部分は元の寸法に戻りました。次は翼端の形状を整えます。[Wing Section]：2[Section Planform]-[Span]：0.025[Section Planform]-[Tip C]：0.22[Sweep]-[Sweep Loc]：0.9[Dihedral]-[Dihedral]：30[Wing Section]：3[Section Planform]-[Span]：0.1[Section Planform]-[Tip C]：0.12[Sweep]-[Sweep Loc]：0.9[Dihedral]-[Rotate Section Tip Foil To Match Dihedral]：ON[Dihedral]-[Dihedral]：60途中経過は以下のようになります。折り目がまだ角張っています。ウィングレットに生じる揚力は形状抗力となり、支持部材に負荷がかかります。空力的に中立が望ましいため、翼根から翼端に向かってNACA4412→2412→0008の断面に小さくしていきます。タブ：Airfoil[Wing Section]：2[Type]-[Camber]：0.02[Wing Section]：3[Type]-[T/C]：0.008[Type]-[Camber]：0トランジション部にブレンドカーブを設定し、主翼とウィングレットを滑らかに接続します。後の設定項目に出てくる"LE"はLeading Edge（前縁）、"TE"はTrailing Edge（後縁）、"TRAP"はTrapezoid（台形）の意味です。トランジション部の翼根側翼断面Airfoil 1を編集します。翼端方向前縁（Outboard）の角度を、翼根方向前縁の台形角度（Inboard_LeadingEdge_Trapezoid）に無理矢理マッチさせます。また、後縁も同様に無理矢理マッチさせます。タブ：Blending[Wing Section]：2[Blend Airfoil]：1[Leading Edge]-[Outboard]-[Match]：IN_LE_TRAP[Leading Edge]-[Outboard]-[Strength]：0.4[Trailing Edge]-[Outboard]-[Match]：IN_TE_TRAP[Trailing Edge]-[Outboard]-[Strength]：0.4トランジション部の翼端側翼断面Airfoil 2を編集します。翼根方向前縁（Inboard）の角度を、翼端方向前縁の台形角度（Outboard_LeadingEdge_Trapezoid）に無理矢理マッチさせます。また、後縁も同様に無理矢理マッチさせます。[Blend Airfoil]：2[Leading Edge]-[Inboard]-[Match]：OUT_LE_TRAP[Leading Edge]-[Inboard]-[Strength]：0.4[Trailing Edge]-[Inboard]-[Match]：OUT_TE_TRAP[Trailing Edge]-[Inboard]-[Strength]：0.4曲面に変換されるのはトランジション部だけで、主翼とウィングレットの形状は元のままです。正面からウィングレットを見るとこうなっています。こんな感じになりました。が、空力的にどうかと聞かれると・・・当たり前のように作業を進めてきましたが、このブレンド機能はユーザーから望まれつつ実装が困難で、Version3.0でようやく追加されたとのことです。このように理想のウィングレットをデザインできるのは良いことですが、実際に翼の先端に作り込むのはなかなか大変そうです。
01Nov
- OpenVSPで無人機をモデリングする
  翼だけ解析してもつまらないので、機体設計ソフトOpenVSPで無人航空機をモデリングしてみようと思います。モデリング手順の紹介なので、模型飛行機ファンや鳥人間コンテストに出場する人が参考にするような記事ではありません。無人機の仕様は次の通りとします。経済的に、長時間滞空できる空撮動画を伝送できるこの記事を書き始めた当初は大きめのラジコンをイメージしていましたが、検討を進めるうちに不整合が生じ、寸法を見直すたびに機体が大きくなってしまいました。現行の法規では動画の伝送は困難なため、静止画を転送する形態になります。この点は現在ドローンメーカー各社から規制緩和の要望が出ているようなので、いずれは利用可能になるかと思います。ライセンスは一等無人航空機操縦士、機体は第一種機体認証が必要なレベル４飛行に該当する中高度滞空型無人機となります。出来上がりは以下のようになります。全幅5.8m、全長2.3mです。何も制約がないと自由すぎて検討を始められないので、まずは推進方式を決めます。速度は求められないので、プロペラ推進のレシプロ機とします。エンジンには構造的に振動が少ない61cc水平対向エンジンを選びます。エンジンユニットの重量は2.36kg、燃料はレギュラーガソリンです。燃費30cc/minだと仮定すると搭載燃料1000ccの９割を消費して30分ほど滞空できますが、機体の役割を考えると大容量の燃料タンクが欲しくなります。観測機材を機首に配置すると、エンジンは機体後部に配置するしかないのでプッシャー式航空機になります。滞空性能を重視するので、主翼のアスペクト比を大きくして揚力を稼ごうと思います。プロペラには22"×10"CCW２枚羽根を使います。直径Dは22inch、幾何ピッチは10inchです。エンジンの回転数7,000rpmのとき、計算上の巡航速度は106.68km/hになります。推力F[N]を以下の式で求めます。これらはUAVエンジニアGabriel Staples氏が２枚羽根での実測データをもとに導き出したものです。静止推力方程式動的推力方程式ここで、pitchは10inch、1min/60secは1/60、対気速度V0[m/s]です。これより静止推力Fは143.9N=14.6kgfと計算されます。参考までにV0が27.78m/s、すなわち100km/h巡航時の動的推力は0.91kgfと計算されます。プロペラは空気を加速する時に推力を生ずるため、巡航速度に近づいて加速度が小さくなると推力は低下します。ピッチをPとして、ピッチ角β3/4を以下の要領で計算します。β3/4 = arctan(P[inch]／(π・3／4・D[inch])) = arctan(10／(π・16.5))= arctan(0.1929) = 10.919[deg]構造強度の維持と揚力分布を改善するためテーパー翼を採用し、主翼の翼幅は片翼2.8m、翼根側翼弦長0.45m、翼端側翼弦長0.25mの中翼とします。主翼の面積Sは1.96m2よりアスペクト比は(2.8×2)2/1.96=16となります。テーパー翼における空力平均翼弦長MACは以下の式によって求めます。翼根の翼弦長t0[m]、翼端の翼弦長t1[m]（但しt0>t1）とすると、MAC= t1 + ((t0 - t1)・(2・t0 + t1) ／ (3・(t0 + t1)))= 0.25 + ((0.2・1.15) ／ (3・0.7)) ≓0.36[m]話が前後しますが、このテーパーの角度は前縁と後縁で異なります。後にOpenVSPで形状を検討した際にSweep Logの値を0.2に設定して翼端を機首側に少し寄せ、翼端前縁を主翼前縁から0.04m後方としました。従って、このMACが翼根側のどの位置にあたるかは、FreeCADで形を描いてみないと分りません。それが以下の図です。赤丸の部分を拡大してみたのが以下の図です。MACの前縁は主翼の前縁から18.1mmの位置にあります。この大雑把な設計では無視できるほど小さな数値になりますが、大抵の主翼形状では無視できない数値になるため、いちおう意識しています。この値は後ほど参考にします。胴体セクションは直径240mmの円形とします。胴体前部は電装品、中央下側は燃料タンク、上側はパラシュートを格納します。その他のスペースはLiFePO4充電池になります。国際標準大気の表によると高度13,000ftの気温は-10.8℃であり、軽量なリチウムイオン充電池は低温に弱いため使えません。全長は最後まで決まらず、最終的に2.3mに落ち着きました。カメラユニットには30倍光学ズーム、5kmレーザー距離計搭載タイプを採用します。このカメラは球体で、直径は138mmです。また、風向風速センサー、GPS/QZSS、ModeSトランスポンダ/ADS-Bも必要です。空気の粘性係数1.7894e-5Pa·s、密度ρ1.225kg/m3、対気速度27.78m/s、空力平均翼弦長0.36mの場合、レイノルズ数Reは6.846e+5となります。翼型は以下のような特性を持つNACA4412とします。広い迎角範囲で高い揚力を得られる、たいへん魅力的な翼型です（有力候補：NACA6412）。主翼取り付け角と上反角はとりあえず0度とします。データ出典：airfoiltools.com尾翼の翼幅は片翼600mmと決め打ちします。翼根側の翼弦長300mm、翼端側を150mmのテーパ翼とします。重心の安定をイメージして尾翼の前縁部は主翼の前縁部から0.9m後方に配置し、主翼の前縁部は機首から1mに配置します。手元に構造設計の資料がないので重心の設計はできませんが、安定性を確保するため主翼の重心位置を23%MACとします。0.36×0.23=0.083[m]、MAC前縁は主翼前縁の18.1mm後方にあることを加味して主翼の重心は主翼前縁の後方およそ0.1mとなります。タンク内燃料の影響が不明ですが、とにかく機首から1.1mの位置より少し後ろ（正確には主翼空力重心の後ろ、かつ全機空力中心の前ですが計算は省略）に重心が来るよう機材の重量を配分する必要があります。エンジンの吸気ダクトを避け、UAVで主流になりつつあるV字型、水平から30度と鉛直方向に３枚配置します。翼型はNACA0012とします。下側の尾翼は着陸装置の邪魔になるため、翼幅300mmで設計します。吸気ダクトは解析の障害になるのでモデリングしません。おおよそ翼周りの形状ができたので、OpenVSPを使って形を作っていきます。各部の寸法は全てパラメータの調整によって決まるため、空力に関わる機体形状を精確にモデリングすることができます。作業に取りかかる前に注意事項があります。ファイルを自動でセーブする機能はないので、編集中はこまめにデータを保存してください。また、起動直後Unnamed.vsp3というファイル名で保存する用意をしています。このためファイルを開くときに誤って"保存"を選ぶと、空のデータであっさり上書きされてしまいますので慎重に操作してください。モデルの作成を開始したら、早い段階でファイルを保存します。今回はファイル名をmyUAV.vsp3としました。このソフトは単位系を自由に決めて良いことになっていますが、解析段階では単位系が決まっていますので最初からm, N, kg, sといったMKS単位系で設計します。OpenVSPを起動したら、主翼を作るためGeom BrowserでWINGを追加します。デフォルトで後退翼が配置されますので、各部のパラメータを調整して欲しい形に整えます。タブ：Sect[Section Planform]-[Span]：2.8[Section Planform]-[RootC]：0.45[Section Planform]-[Tip C]：0.25[Sweep]-[Sweep]：0[Sweep]-[Sweep Loc]：0.2タブ：XForm[Transforms]-[YLoc]：0.1主翼の翼型NACA4412を定義します。翼根側、翼端側のセクション両方を設定します。一方のセクションを設定すれば、もう一方のAirfoil Sectionに設定をコピー＆ペーストすることができます。タブ：Airfoil[Airfoil Section]：Section 0 および Section 1[Type]-[T/C]：0.12[Type]-[Camber]：0.04[Type]-[CamberLoc]：0.4タブ：Plan[Tip Treatment]-[Tip Cap Type]：Round続いて左右尾翼を追加します。Geom BrowserでWINGを追加します。タブ：XForm[Transforms]-[XLoc]：0.9[Transforms]-[YLoc]：0.06タブ：Sect[Section Planform]-[Span]：0.6[Section Planform]-[RootC]：0.3[Section Planform]-[Tip C]：0.15[Sweep]-[Sweep]：0[Sweep]-[Sweep Loc]：0.8[Dihedral]-[Dihedral]：30尾翼の翼型NACA0012を定義します。翼根側、翼端側のセクション両方を設定します。タブ：Airfoil[Airfoil Section]：Section 0 および Section 1[Type]-[T/C]：0.12[Type]-[Camber]：0[Type]-[CamberLoc]：0.1タブ：Plan[Tip Treatment]-[Tip Cap Type]：Round下部尾翼を追加します。Vehicle>WingGeom配下にWINGを追加します。タブ：XForm[Symmetry]-[Planer]：OFF[Transforms]-[XLoc]：0.9[Transforms]-[ZLoc]：-0.06[Transforms]-[XRot]：-90タブ：Sect[Section Planform]-[Span]：0.3[Section Planform]-[RootC]：0.3[Section Planform]-[Tip C]：0.15[Sweep]-[Sweep]：0[Sweep]-[Sweep Loc]：0.8同じく尾翼の翼型NACA0012を定義します。翼根側、翼端側のセクション両方を設定します。タブ：Airfoil[Airfoil Section]：Section 0 および Section 1[Type]-[T/C]：0.12[Type]-[Camber]：0[Type]-[CamberLoc]：0.1タブ：Plan[Tip Treatment]-[Tip Cap Type]：Round胴体を追加します。Geom BrowserでTransportFuseを追加します。タブ：Design[Design]-[Length]：2.3[Design]-[Diameter]：0.24[Design]-[Aft Mult]：2.5[Design]-[Aft Center]：0[Design]-[Aft Width]：0.5[Design]-[Aft Height]：0.5タブ：XForm[Transforms]-[XLoc]：-1プロペラを追加します。Geom BrowserでPROPを追加します。プッシャ式航空機でもトラクタ式と推力方向は変わりませんが、エンジンが後ろ向きのため回転方向が逆になります。このためプロペラは羽根が逆向きのCCWタイプを使います。今回は最低限のパラメータしか設定しませんが、多くの項目があり、資料さえあれば厳密に定義することができます。タブ：Degin[Design]-[Diameter]：0.5588[Design]-[Num Blades]：3[Design]-[Cylindrical Foils]-[Rev]：ON[Design]-[Cylindrical Foils]-[Beta 3/4]：10.919[Modeling]-[Prop Mode]：Bothタブ：XForm[Transforms]-[XLoc]：1.3最後にカメラポッドを追加します。Geom BrowserでELLIPSOIDを追加します。タブ：XForm[Transforms]-[XLoc]：-0.8[Transforms]-[ZLoc]：-0.15タブ：Design[Design]-[A Radius]：0.069[Design]-[B Radius]：0.069[Design]-[C Radius]：0.069以上でモデリングが完了しました。記事の冒頭で紹介したモデルと全く同じものができたと思います。ここまでできるとどんなモデルでも作れるように思えますが、実際はそうでもないようです。たとえば高揚力装置などのモデリングは案外難しく、当のNASAでさえも既存の機能を組み合わせて実現した方法がこちらの論文になっています。最後に、付属の解析ソフトVSPAEROでごく簡単な解析をしてみます。解析時にプロペラが邪魔になるため、プロペラを一時的に消します。Geom BrowserのリストからPropGeomを選択し、[Selection]-[NoShow]をクリックするとプロペラモデルが非表示になります。メニューから[Analysis]-[Aero]-[VSPAERO...]を選択し、VSPAEROを起動します。パラメータは最低限のものしか調整しません。タブ：Overview[Reference Area, Lengths]-[From Model]：ON　// 基準面積をモデルから取得[Reference Area, Lengths]-[Ref. Wing]：0_WingGeom　// 主翼を指定[Case Setup]-[Geometry Set:]：Shown　// 表示中の翼モデルのみ解析[Frow Condition]-[Alpha Start]：0　// 機体ピッチ0度から[Frow Condition]-[End]：10　// +10度まで[Frow Condition]-[Npts]：3　// 3点を解析[Frow Condition]-[ReCref Start]：684600　// 解析対象のレイノルズ数Launch Solverをクリックすると解析が始まります。解析が終わるとResults Managerが開きます。Load Dist.タブの[Y-Data]：cl*c/crefを選択すると、以下のグラフが表示されます。よく似たデータclと異なるのは、翼全体の揚力に対してそのスパン位置で生じた揚力の貢献具合を加味したものです。clに対してそのスパンにおける翼弦長cを平均翼弦長crefで割った係数を掛けたもので、局所的にcl値が大きくてもその翼弦長が短いとあまり役に立たないと分ります。いま、Flow Conditionでは姿勢の異なる３つのケース全てを選択、グラフ表示しています。Legendにある通り、機体のピッチについて赤色が0度、黄色が5度、緑色が10度での翼幅方向の揚力係数分布を示しています。この機体の設計では主翼取り付け角を0度としたので、機体のピッチは主翼の迎角でもあります。このグラフより、主翼の揚力係数CLは0.35程度とわかりました。また別のグラフ[Y-Data]：cdによると、抗力係数CDは0.004程度でした。これより揚抗比は0.35/0.004=87.5となります。尾翼の特性に少し気になるところがありますが、ここでは追求しません。以上の解析により、高度0m、機体のピッチ0度、対気速度が100km/hのとき、揚力Lは以下のように計算できます。L = 1／2・空気密度ρ[kg/m3]・(対気速度U[m/s])2・翼面積S[m2]・揚力係数CL= 1／2・1.225・(27.78)2・1.96・0.35= 324.26[N] = 33.06[kgf]同様に高度3000mにおける空気密度0.909kg/m3では24.53kgf、高度4000mでは22.11kgfとなります。高度が上がるほど空気密度が低下して揚力も減少し、またエンジンの燃焼効率も低下します。中高度を滞空する航空機として高度3000mは譲れませんし、高度4000mでは搭載したセンサーの仕様を超えますので、実用上昇限度はその辺りとなります。これまでまっすぐ飛ぶ前提で考えてきましたが、水平定常旋回には余分に揚力が必要です。旋回するたびに高度が下がると観測に支障が出かねません。バンク角20度のとき翼の揚力はcos20°=0.94倍に減じるため、揚力をそのぶん割り引いて見る必要があります。これらの事情より、目標とする機体重量Wは24.53×0.94=23kgfとします。このとき推力重量比は14.6/23=0.63、翼面荷重は23/1.96=11.73kgf/m2です。推力重量比0.63は一般的なモーターグライダーの数値0.1～0.3と比べてかなり大きめです。では、もう一段小さい41ccのエンジンではどうでしょうか。先に示した静止推力方程式で計算すると、18"×10"のプロペラ7,000rpmにおける静止推力は7.27kgfです。それでも推力重量比7.27/23=0.31を確保しているので、41ccでも操縦性に問題はなさそうです。41ccのエンジンユニットは1.6kgで、61ccのエンジンより0.76kg軽量です。バンク角Φ20度における旋回半径rと旋回率ωは、真対気速度VTAS=27.78m/s=54ktとして以下のように求めます。r = VTAS2／(68579・tanΦ) = (54)2／(68579・tan20°)= 0.1168[NM] = 187.97[m]ω = 1092・tanΦ／VTAS = 1092・tan20°／54 = 7.36[deg/sec]最低飛行速度Vminは、VSPAEROの解析結果より機体のピッチ+10度での揚力CL1.3として以下のように求めます。Vmin2 = (2／(ρ・CL))・(W／S) = (2／(1.225・1.3))・(23／1.96) = 14.737Vmin = 3.83[m/s] = 13.78[km/h]ただしこの時の回転数を計算すると904rpmであり、エンジンの実用回転数1,500rpmを下回るため意図的にこの速度で飛ぶことはできません。近年UAV向けに発電機を追加したエンジンも出てきているため、うまく活用すれば充電池の容積を減らせるかもしれません。まだまだOpenVSPの全容を把握しきれていませんが、試行錯誤しながらモデルを調整していくのだろうなという感覚はつかめました。
24Oct
- OpenFOAMで翼型を三次元解析する
  今回は実寸の翼モデルを作成して三次元解析をしてみます。三次元解析では、二次元解析にはなかった誘導抗力が生じます。計算領域を狭めるため翼幅を短く定義すると、アスペクト比が小さくなるため誘導抗力の影響が大きくなります。このため、計算コストはかかりますが実寸の翼幅で解析することにします。翼型NACA2412の翼弦長を1.49m、翼幅9.8m、対気速度60kt、迎角10度における揚力と抗力を求めます。翼弦長1の翼型データNACA2412.datは、前回作成したものを使います。翼モデルの作成FreeCADでNACA2412.datを読み込み、[View]-[Toggle axis cross]をオンにします。ワークベンチ：Draft[Scale]コマンドを実行します。Global X, Y, Z：0, 0, 0X factor：1490Y factor：1490Modify original：オン[Base]-[Placement]-[...]をクリックします。Translation X：-372.5Center X：372.5Axis Z：-10ここでOKし、姿勢変更をいったん確定します。次に、サーフェスの向きをTOP面からFRONT面に倒します。[Base]-[Placement]-[...]をクリックします。Apply incremental changes：オンAxis X：90モデリング画面右上の視点キューブをFRONT視点に切り替え、作業を続行します。ワークベンチ：PartModelからWireを選択します。[Part]-[Extrude...]をクリックします。[Direction]-[Custom direction]を選択します。X, Y, Z：0, -1, 0Length-Along：9800Symmetric：オンCreate solid：オン風洞の寸法は全長40m、全幅16m、高さ16mとし、翼モデルを高さ2mに配置します。ModelからExtrudeを選択します。[Base]-[Placement]-[...]をクリックします。Translation Z：2000ワークベンチ：DraftModelからExtrudeを選択します。[Draft]-[Downgrade...]を２回実行すると、Solidを経てFaceに細かく分解されます。ここで、全てのFaceが選択されていることを確認してください。[File]-[Export...]にて一旦AST形式としてNACA2412.astを書き出します。これは、STL形式でエクスポートするとバイナリ形式となってしまうためです。OpenFOAMで扱う時はASCII形式である必要があります。書き出したNACA2412.astをNACA2412.stlに拡張子を書き換えます。以上でFreeCADによる作業は終わりです。チュートリアルファイルのコピーバイク周りの流れを解析するOpenFOAM11のチュートリアルcaseファイルを流用します。このcaseは72km/hで走行するバイクをシミュレートしたものです。lowerWallすなわち地面も72km/hで流れる移動壁となっており、この解析条件は滑走路の上空を飛行する航空機をシミュレートするのに適しています。乱流モデルはk-ω SSTです。ファイル一式を作業用ディレクトリ$FOAM_RUNに取り込んだ後、motorBikeSteady/0フォルダ内のU.origファイルをUにリネームします。$ cd /opt/openfoam11/tutorials/incompressibleFluid$ cp -r motorBikeSteady $FOAM_RUN$ cd $FOAM_RUN/motorBikeSteady/0$ mvU.orig U$ cd $FOAM_RUNチュートリアル内に収録されているバイクモデルは読み込む必要がないので、motorBikeSteady/Allrunの以下の行を#でコメントアウトします。# cp $FOAM_TUTORIALS/resources/geometry/motorBike.obj.gz constant/geometry/motorBikeSteady/Allcleanを開き、以下の行を#でコメントアウトします。# rm -f constant/geometry/motorBike.obj.gz作成した翼モデルNACA2412.stlをconstant/geometryフォルダに配置し、翼モデルの単位系をmに変更するため以下のコマンドで翼モデルの寸法を1/1000縮小します。同時に、ファイル形式をSTLからOBJにコンバートします。$ cd $FOAM_RUN/motorBikeSteady/constant/geometry$ surfaceConvert NACA2412.stl NACA2412.obj -clean -scale 0.001$ rm NACA2412.stlファイル名をmotorBike.objに変更します。なぜなら、motorBikeという文字列はcaseファイルのいたる所で使用されており、これらを見つけて全て書き換える手間を省くためです。編集点を最小限にとどめることで、ファイル比較ソフトで元ファイルと編集済みファイルの異なる箇所を見比べやすくなります。$ mv NACA2412.obj motorBike.objfindコマンドを使って文字列を一気に書き換えるテクニックはありますが、ここでは説明しません。計算に用いるCPUコア数の設定このチュートリアルは計算負荷が高いらしくCPUコア6個のマシンでの計算を前提にしており、このコア数より少ないマシンでは計算できません。ここでは大抵の家庭用パソコンで対応できるであろうCPUコア2個に設定します。system/decomposeParDictファイルをテキストエディタで開き、以下のように編集します。numberOfSubdomains 2;hierarchicalCoeffs{ n (21 1); // Nx * Ny = number of CPU cores order xyz;}風洞モデルの調整system/blockMeshDictファイルをテキストエディタで開き、以下のように編集します。チュートリアルの評価対象であったバイクより翼が大きいため、風洞の寸法とメッシュの数をそれぞれ２倍しています。vertices( (-10 -8 0) ( 30 -8 0) ( 30 8 0) (-10 8 0) (-10 -8 16) ( 30 -8 16) ( 30 8 16) (-10 8 16));blocks( hex (0 1 2 3 4 5 6 7) (40 16 16) simpleGrading (1 1 1));system/snappyHexMeshDictファイルを開き、以下のように編集します。翼モデルを包むようにメッシュを細分化します。翼の風下側は気流の乱れが尾を引くため、20mと長めに距離を取っています。refinementBox{ type searchableBox; min (-4.0 -8.0 0.0); max (20.0 8.0 5.0);}力係数を精度よく計算するため、refinementSurfaceの細分化レベルをもう一段引き上げます。motorBike{ level (6 7); ... 中略 ...}ついでにaddLayersControlsで翼まわりに３層のメッシュを挿入します。このとき対気速度と同じ速度の移動壁となる地面も、同様にメッシュが細かくなります。これは、翼に働く地面効果を観察する際に役立つものと思われます。layers{ "(lowerWall|motorBike).*" { nSurfaceLayers 3; }}0/include/initialConditionsファイルを開き、以下のように編集します。対気速度60ktは30.86m/sに相当します。flowVelocity (30.860 0);system/forceCoeffsファイルを開き、以下のように編集します。CofRは回転モーメントの中心を指定するものです。翼弦長25%位置を翼のX軸とY軸の原点としており、かつZ位置を高さ2mと定義しています。CofR (00 2); // Centre of rotation for moment calculationspitchAxis (0 1 0);magUInf 30.86;lRef 1.49; // Mean Aerodynamic ChordAref 14.38; // 1.49 * cos10deg * 9.8計算開始caseファイルの編集ができたら、次のコマンドで計算を開始します。$ cd $FOAM_RUN/motorBikeSteady$ ./Allrun解析結果の確認解析結果はparaFormで観察することができます。$ paraFoamメッシュ高さ1.8mに流線を表示させると、翼の後方に翼端渦が生じている様子が見えます。本来なら翼の上面を高速で通過し負圧を発生するはずの気流が下面の外向きに流れており、圧力差が小さくなっていることが分かります。これが誘導抗力として働き、揚力を減少させたり余分にエネルギーを消費したりする原因になります。こういう現象を目の当たりにすると、実際に翼を制作する時は何かしらの対策――欲を言えば見栄えのいいウィングレット――が欲しくなってきます。得られた結果を正確に評価するには、まだまだ考慮すべきこと、見直すべきパラメータが多いので、解析結果については評価しません。
07Oct
- OpenFOAMで任意の翼型を評価する
  OpenFOAM11で任意の翼型を評価します。必要なソフトは次の通りです。 OpenFOAM11 ・・・公式手順でインストールする FreeCAD 0.18 ・・・バージョンがこれより古くても新しくてもバグ有り Anaconda3 ・・・ Python3が実行できれば何でも良い代表的な軽飛行機セスナ172の主翼に使われているNACA2412を評価対象とします。翼型データの取得ここからNACA翼型プロファイルジェネレータを入手します。Anaconda Powershell Promptを立ち上げ、以下のようにnaca.pyを実行すると、NACA2412の座標データがNACA2412.datに書き出されます。python naca.py -p 2412 -n 500 | Out-File -Encoding utf8 NACA2412.dat生成されたNACA2412.datをメモ帳などのテキストエディタで開き、以下のように編集します。1行目はヘッダーで、翼型の名称を記述します。2行目以降が座標データです。座標は後縁部から記述し、翼の上面を前縁部まで移動し、下面を通って後縁部に戻るという反時計周りの順で並べることになっています。ジェネレータで生成した後縁部のY座標は限りなくゼロに近い(が、ゼロではない)ため、0.0に書き換えます。最終行の1002行目は改行しません。1行目：NACA 24122行目：1.0 0.03行目：0.9980192478918325 0.0004227951019013875（中略）1001行目：0.9979807521081675 -0.000156572879679165021002行目：1.0 0.0編集し終えたら、エンコードをUTF-8(BOM無し)として上書き保存します。これで、翼弦長1の翼型データができました。翼モデルの作成FreeCADを立ち上げ、[File]-[Open...]でNACA2412.datを選択すると以下のメニューが現れるのでCommon airfoil dataを選択します。以降はワークベンチを切り替えながら作業します。翼型データを読み込んだ直後はTOP視点にサーフェスを展開されていますので、サーフェスの向きをFront側に倒します。ワークベンチ：Draft[Base]-[Placement]の右側にある三点リーダーをクリックし、[Rotation]-[Axis]をXに切り替えます。[Rotation]-[Angle]に90を入力してOKします。モデル右上のキューブを操作してFrontを表示すると、NACA2412の翼型が見えるようになります。ここで、[View]-[Toggle axis cross]をクリックしてXYZ軸を表示させておきます。ワークベンチ：Part[Part]-[Extrude...]を実行します。Custom direction：0, -1, 0[Length]-[Along]：1Symmetric：オンCreate solid：オン以下のようになればOKです。翼弦長1m、原点を中心として翼幅1mとなります。解析対象の翼弦長はレイノルズの相似則を利用し、流速の設定を以てそのスケールを調整します。ワークベンチ：Draft画面左のModelビューでExtrudeを選択し、[Draft]-[Downgrade]を２回実行します。この操作でモデルがFaceの集合体に分解されます。ModelビューでFaceが全選択されていることを確認して、[File]-[Export...]を実行してWavefront OBJ形式のNACA2412.objを出力します。チュートリアルファイルのコピーOpenFOAMをインストールしたときに、乱流モデルにk-ωSSTを使用したチュートリアルが用意されています。この中に定常状態を解析するaerfoilNACA0012Steadyというケースファイルがありますので、これを作業フォルダにコピーします。また、フォルダ名がそのままだと翼型が混乱するため、フォルダ名を変更します。$ cd /opt/openfoam11/tutorials/fluid$ cp -r aerofoilNACA0012Steady $FOAM_RUN$ cd $FOAM_RUN$ mv aerofoilNACA0012Steady aerofoilNACA2412Steady翼モデルの組み込み先ほど作成したNACA2412.objをaerofoilNACA2412Steady/constant/geometry/フォルダにコピー＆ペーストします。Allrunファイルを開き、下記の行を#でコメントアウトして保存します。この行は解析実行時に都度翼モデルをコピーしてくるコマンドですが、さきほど解析したい翼モデルをgeometry配下にコピーしたので不要になります。cp $FOAM_TUTORIALS/resources/geometry/NACA0012.obj.gz constant/geometry/↓# cp $FOAM_TUTORIALS/resources/geometry/NACA0012.obj.gz constant/geometry/Allcleanファイルを開き、下記の行を#でコメントアウトして保存します。この行は初期化に伴って翼モデルも削除するコマンドですが、この操作も不要です。rm -f constant/geometry/NACA0012.obj.gz↓# rm -f constant/geometry/NACA0012.obj.gzaerofoilNACA2412Steady/system/blockMeshDictファイルを開き、翼モデルを指定している55行目を以下のように変更し、保存します。file "NACA0012.obj";↓file "NACA2412.obj";解析条件の設定aerofoilNACA2412Steady/0/Uファイルを開き、流速と迎角を設定します。空力平均翼弦1.49mのセスナ172が高度0km、対気速度60kt（30.86m/s）の離陸体制にある場合を例に取ります。実機ではこの対気速度を超えたあたりで機首の引き上げ操作を行います。解析したい流速を項目speedに設定します。単位は[m/s]です。この時のレイノルズ数はおよそ3.1e+6であり、翼弦長1m換算で代表流速は45.28m/sとなります。NACAがこのレイノルズ数のデータを取っていることから、彼らが軽飛行機の離陸性能をとくに気にしていたことが窺えます。解析したい迎角を項目angleOfAttackに設定します。単位は[deg]です。上向き10度の時は10となります。aerofoilNACA2412Steady/system/controlDictファイルを開き、26行目のendTime（上限イテレーション数）を設定します。ここでは10000としています。収束が早かった場合は、10000に到達する前に計算を終了する場合があります。途中経過はおおざっぱに確認できれば良いので、32行目のwriteIntervalは200とします。途中で計算が停止した場合は、収束したか発散した可能性がありますので、心配なときはlog.foamRunファイルのログを確認してください。解析開始以下のコマンドで解析が始まります。$ cd $FOAM_RUN/aerofoilNACA2412Steady$ ./Allrun解析をやり直す際は、再度Allrunを実行する前に以下のコマンドでcaseフォルダを初期化します。$ ./Allclean解析結果の確認解析結果はparaFormで観察することができます。以下は迎角10度での画像です。$ paraFoamメッシュ流速分布静圧分布流線表示揚力と抗力はaerofoilNACA2412Steady/postProcessing/forceCoeffsCompressible/0/forceCoeffs.datに書き出されています。最後の行に書き出された数値が解となります。なお、上記の条件で解析を行い、Time 1108で収束した時の数値はCL：1.2338432702e+00CD：1.6032227531e-02となっています。これより、迎角10度における揚抗比(CL／CD)は76.96となります。揚力Lの計算式は以下の通りです。機体の全幅11mから空力に寄与しないキャビン幅1.2mを引いて翼幅9.8mとします。L = 1／2・空気密度ρ[kg/m3]・(代表流速U[m/s])2・代表面積A[m2]・揚力係数CL= 1／2・1.225・(30.86)2・1.49・cos10°・9.8・1.233843= 10,350[N] = 1055.40[kgf]セスナ172の最大離陸重量は1157kgです。対気速度60ktでは乗員、貨物、燃料を満した状態で離陸できないという結果になりました。とはいえ揚力に対して代表流速が二乗で効くことが幸いして、あと少し、63ktまで加速すれば離陸できる計算です。迎角を変えながらOctaveでプロットすると以下のようになります。ネットで手に入る各所の実測結果に照らして、それらに近い数値が得られたように思います。この手法は、チュートリアルに使われているNACA0012モデルを差し替えるだけです。したがって、どんな翼型でも翼弦長1で作成できれば高い精度で揚力と抗力を解析することができます。

前ページ
次ページ