こんにちは（@t_kun_kamakiri)
メモに残しておきます。
たまに追加していくので困ったときに知っておくと役に立つかもしれません。
気になる部分は目次からジャンプしてください_(._.)_
OpenFOAMv2212
流体解析で計算がうまくいかないとき
OpenFOAMをメインに流体解析でうまくいかないときの原因や対策などのメモを残しています。
ご参考ください_(._.)_
bashのプロンプトを変更する
ターミナルのパスが長すぎたりすると改行したり見づらかったりしますよね。
今回は以下のように時刻と改行を入れた設定に変えるようにします。
以下をターミナルで打って、
1 | vi ~/.bashrc |
一番最後の行に下記を追加します。
1 | PS1='\ek\e\\\][\e[32m\u@\h \e[35m\t \e[34m\w\e[0m]\n\$ ' |
「:wq」と打てば保存してvimから抜けることができます。
↓こちらにするとパスの文字が黄色になります。
1 | PS1='\ek\e\\\][\e[32m\u@\h \e[36m\t \e[33m\w\e[0m]\n\$ ' |
好きなカラーで設定してください。
OpenFOAMグループ
こちらに質問と回答が蓄積されているので同じ悩みの方に出会えるかもしれません。
ただ良い回答を得るためには良い質問の仕方があります。
snappyHexMeshでメッシュが思った通りに生成されていないとき
- 「locationInMesh (2.1 0.09 -1.1);」の座標がおかしい
- メッシュサイズが大きくて形状変化の解像度をとらえられていない
- stlファイルをよく見ると欠けている
だいたいどれか。
並列計算しているなら結合しないとParaViewでは見れない。
並列分割してフォルダ内にprocessorがあるなら、ParaviewのCase TypeをDecomposed Caseに変えると並列分割した状態でも読み込めます。
流量を知りたい
①に残差の例もあります。
残差と連続式の誤差を出力
system/controlDictには残差や連続式の誤差を出力する設定を記述します。
system/controlDict
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 | application simpleFoam; startFrom startTime; startTime 0; stopAt endTime; endTime 500; deltaT 1.0; adjustTimeStep no; maxCo 0.9; writeControl timeStep; writeInterval 50; purgeWrite 0; writeFormat ascii; writePrecision 6; writeCompression off; timeFormat general; timePrecision 6; runTimeModifiable yes; functions { continuityError1 { type continuityError; libs ( fieldFunctionObjects ); // Mandatory entries (unmodifiable)// Optional entries (runtime modifiable) phi phi; // Optional (inherited) entries writePrecision 8; writeToFile yes; useUserTime yes; region region0; enabled yes; log yes; timeStart 0; timeEnd 6000; executeControl timeStep; executeInterval 1; writeControl timeStep; writeInterval 50; } residuals { type solverInfo; libs ("libutilityFunctionObjects.so"); fields (U p k omega); } //#includeFunc residuals(U,p,k,epsilon); #includeFunc yPlus; #includeFunc CourantNo; } |
残差と連続の式による誤差は最低確認するようにしましょう。
連続の式の誤差
- sum local : 誤差絶対値の格子体積重み付け平均
- global : 誤差(符号あり)の格子体積重み付け平均
- cumulative : globalの累積
グラフを出力するPythonスクリプトはこちらです。
連続の式の誤差
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 | import numpy as np import matplotlib.pyplot as plt import pandas as pd import os def graph_func(df_continuityError, imin, imax, iposi,icol): plt.subplot(imin,imax,iposi) plt.plot(df_continuityError['time'],df_continuityError[icol],linestyle='solid',color='b',marker="o",label=icol) plt.xlabel('time') plt.ylabel(icol) plt.legend() plt.grid() def df_continuityError_func(dir_): data_continuityError = np.loadtxt(f'./postProcessing/continuityError1/{dir_}/continuityError.dat') df_continuityError = pd.DataFrame(data_continuityError, columns=['time','Local','Global','Cumulative']) return df_continuityError def df_concat(dir_List): df = pd.DataFrame() for dir_ in dir_List: try: df_ = df_continuityError_func(dir_) df = pd.concat([df, df_]) except Exception: print("Empty") return df dir_List = os.listdir('./postProcessing/continuityError1') df_continuityError = df_concat(dir_List) plt.figure(figsize=(20,3)) graph_func(df_continuityError, 1,3,1,'Local') graph_func(df_continuityError, 1,3,2,'Global') graph_func(df_continuityError, 1,3,3,'Cumulative') print('time : ',df_continuityError['time'].iloc[-1], 'continuityError : ', df_continuityError['Local'].iloc[-1]) df_continuityError.head() |
残差
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 | import pandas as pd import matplotlib.pyplot as plt import os def graph_layout(): plt.grid() plt.legend(loc='best',bbox_to_anchor=(1, 1)) plt.yscale('log') plt.xlabel('time', fontsize=16) plt.ylabel('residual', fontsize=16) def df_residual_func(dir_): df_residual = pd.read_table(f'./postProcessing/residuals/{dir_}/solverInfo.dat',skiprows=1) df_residual = pd.DataFrame(df_residual) return df_residual def df_concat(dir_List): df = pd.DataFrame() for dir_ in dir_List: try: df_ = df_residual_func(dir_) df = pd.concat([df, df_]) except Exception: print("Empty") return df dir_List = os.listdir('./postProcessing/residuals') df_residual = df_concat(dir_List) initial_residial = [data for data in df_residual.columns if "initial" in data] final_residial = [data for data in df_residual.columns if "final" in data] df_residual.plot(x="# Time ",y=initial_residial, figsize=(6, 4)) graph_layout() #df_residual.plot(x="# Time ",y=final_residial) |
残差の分布を出力す量にカスタマイズした有益な議論があります。
実験データと比較したい
だいたいここを見ます。
境界条件を楽にする
共通の設定の場合
1 2 3 4 5 | "(inlet|outlet)" { type fixedValue; value uniform (0 0 0); } |
デフォルト(全て壁条件)
1 2 3 4 | ".*" { type noSlip; } |
設定の使いまわし
1 2 3 4 5 6 7 8 9 10 | inlet1 { type fixedValue; value uniform (1 0 0); } inlet2 { $inlet1; } |
並列計算
1 2 3 | decomposePar reconstructParMesh -constant -mergeTol 1e-6 mpirun -np 4 renumberMesh -overwrite -parallel |
特定のフォルダの容量調べる
1 2 3 4 | du -sh /mnt/d/OpenFOAM/ #出力結果 50G /mnt/d/OpenFOAM/ |
境界条件の候補を調べる
1 | foamHelp boundary -field U |
以下が出力される。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 | Create mesh for time = 0 Available boundary conditions for vector field: U SRFFreestreamVelocity SRFVelocity SRFWallVelocity acousticWaveTransmissive activeBaffleVelocity activePressureForceBaffleVelocity advective calculated codedFixedValue codedMixed cyclic cyclicACMI cyclicAMI cyclicSlip cylindricalInletVelocity directionMixed empty exprFixedValue exprMixed extrapolatedCalculated fixedGradient fixedInternalValue fixedJump fixedJumpAMI fixedMean fixedMeanOutletInlet fixedNormalInletOutletVelocity fixedNormalSlip fixedProfile fixedShearStress fixedValue flowRateInletVelocity flowRateOutletVelocity fluxCorrectedVelocity freestream freestreamVelocity inletOutlet interstitialInletVelocity kqRWallFunction mapped mappedField mappedFixedInternalValue mappedFixedPushedInternalValue mappedFlowRate mappedMixed mappedMixedField mappedVelocityFlux matchedFlowRateOutletVelocity mixed movingWallVelocity noSlip nonuniformTransformCyclic outletInlet outletMappedUniformInlet outletPhaseMeanVelocity partialSlip pressureDirectedInletOutletVelocity pressureDirectedInletVelocity pressureInletOutletParSlipVelocity pressureInletOutletVelocity pressureInletUniformVelocity pressureInletVelocity pressureNormalInletOutletVelocity pressurePIDControlInletVelocity processor processorCyclic rotatingPressureInletOutletVelocity rotatingWallVelocity scaledFixedValue sliced slip supersonicFreestream surfaceNormalFixedValue swirlFanVelocity swirlFlowRateInletVelocity swirlInletVelocity symmetry symmetryPlane timeVaryingMappedFixedValue translatingWallVelocity turbulentDFSEMInlet turbulentDigitalFilterInlet turbulentInlet uniformFixedGradient uniformFixedValue uniformInletOutlet uniformJump uniformJumpAMI uniformNormalFixedValue variableHeightFlowRateInletVelocity waveTransmissive wedge zeroGradient End |
ソースコードを探す
全ての$FOAM_APPを出力しておいて、laplacianFoam.Cの検索に引っかかったものだけを抽出
1 | find $FOAM_APP |grep "laplacianFoam.C" |
$FOAM_APPの中で拡張子.Cだけを出力しておいて、laplacianFoam検索に引っかかったものだけを抽出
1 | find $FOAM_APP -name *.C | grep "laplacianFoam" |
検索したい文字列のファイルを探す
1 2 | # find ファイルパス -type f | xargs grep -n '検索したい文字列' find $FOAM_TUTORIALS -type f| xargs grep -ln "noSlip" |
- l：ファイル名のみ出力
- n：行番号を出力
メッシュ品質の基準を検索
checkMeshコマンドでメッシュ情報を確認すると、メッシュ数の他にメッシュ品質に関わる内容も出力されます。
しかし、「メッシュ品質って何なの？」となることもあるので、以下に基準値が記載してあるファイルのありかを調べる方法を記します。
1 | find $FOAM_SRC -name "primitiveMeshCheck" |
以下にあることがわかる。
1 | /opt/OpenFOAM/OpenFOAM-v2012/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck |
viコマンドで中身を確認。
1 | vi /opt/OpenFOAM/OpenFOAM-v2012/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C |
品質基準は以下を参照。
1 2 3 4 5 | Foam::scalar Foam::primitiveMesh::closedThreshold_ = 1.0e-6; Foam::scalar Foam::primitiveMesh::aspectThreshold_ = 1000; Foam::scalar Foam::primitiveMesh::nonOrthThreshold_ = 70; // deg Foam::scalar Foam::primitiveMesh::skewThreshold_ = 4; Foam::scalar Foam::primitiveMesh::planarCosAngle_ = 1.0e-6; |
メッシュ品質は過去にオープンCAEのサマースクールで発表した際の資料にも書いてあります。
メッシュ品質の悪い部分を確認
メッシュ品質が悪い箇所のメッシュを修正したい場合、ParaViewでその場所を特定することができます。
まず、以下でメッシュ品質項目のvtkファイルを作成します。
1 2 | foamToVTK -faceSet skewFaces -time 0 foamToVTK -faceSet nonOrthoFaces -time 0 |
ls constant/polyMesh/setsで品質の悪い項目のファイルを確認。
constant/polyMesh/sets/ファイル名
ファイル名 = skewFaces
次に、paraviewでvtkを読み込むとメッシュ品質の悪い箇所を確認できます。
- 高アスペクト比（Aspect Ratio）
非常に微細な境界層で現れる。
ソルバーの安定性にとって致命的なものではありませんが、収束速度を著しく低下させる可能性があり - 非直交性（Orthogonality）
3つの値の範囲を定義する
・n0 < 70 – 安全な値
70 < n0 < 90 – fvSolutionのnonOrthoCorrectorやfvSchemesの数値スキームに対して特別な処理を行う必要あり
・n0 > 90 – シミュレーションに使用できない悪いメッシュ - 歪度(skewness)
値が大きいと結果の品質（正確さ）が損なわれる可能性がありますが、適度な大きさの歪度であれば、シミュレーションに使用することができる - Smoothness
- メッシュ間隔の最大変化は20%未満が理想
viコマンド
例えばメッシュ品質基準を調べたい場合。
1 | vi /opt/OpenFOAM/OpenFOAM-v2012/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C |
#行数を表示
:set number
/primitiveMesh
n：下へ検索
N(Shift + n)：上へ検索
- :q （上書きせずに閉じる）
- :q! （編集した場合でも上書きせずに閉じる）
- :wq （上書き保存）
- iでその場所から編集
- Aで最終列に移動して編集
- uで戻す
圧縮ファイルの解凍
tar.gz
1 2 3 4 | 圧縮 tar -zcvf xxxx.tar.gz directory 解凍 tar -zxvf xxxx.tar.gz |
zip
1 2 3 4 | 圧縮 zip -r xxxx.zip directory 解凍 unzip xxxx.zip |
treeコマンドが使えないとき
1 | pwd;find . | sort | sed '1d;s/^\.//;s/\/\([^/]*\)$/|--\1/;s/\/[^/|]*/| /g' |
OpenFOAMの環境変数の確認
1 2 3 | echo $FOAM_TUTORIALS # 出力結果 /opt/OpenFOAM/OpenFOAM-v2012/tutorials |
ケースファイルの初期化
初期化したい内容に応じてコマンドを使い分ける。
1 2 3 4 5 6 7 8 | #設定の初期化 foamCleanTutorials #メッシュを初期化 foamCleanPolyMesh #計算結果の初期化 foamListTimes -rm |
postProcess
controlDictで設定。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 | Q1 { // Mandatory entries (unmodifiable) type Q; libs (fieldFunctionObjects); // Optional (inherited) entries field <inpField>; result <fieldResult>; region region0; enabled true; log true; timeStart 0; timeEnd 1000; executeControl timeStep; executeInterval 1; writeControl timeStep; writeInterval 1; } |
コマンド実行
1 | postProcess -func Q |
とすると時刻のフォルダにQが出力
sampling
samplingによりフィールドデータをサンプリングするためのサンプリング関数オブジェクトが用意されています。グラフにプロットするための1Dライン、画像として表示するための2D平面や3Dサーフェスのいずれかを使用できます。各サンプリングツールは、controlDictファイルのメイン関数辞書、またはcaseシステムディレクトリの別ファイルのいずれかで辞書に指定されています。データは、Grace/xmgr、gnuplot、jPlotのような有名なグラフ作成パッケージを含む様々なフォーマットで書き出すことができます。
以下はsystem/sampleDict内に記述する例を示しています。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 | sampleDict { type sets; libs (sampling); setFormat raw; interpolationScheme cell; fields (U p wallShearStress); writeControl writeTime; sets ( z_0.001m { type face; axis xyz; start (-0.01 0.01 0.001); end ( 0.01 0.01 0.001); } z_0.05m { type face; axis xyz; start (-0.01 0.01 0.05); end ( 0.01 0.01 0.05); } z_0.1m { type face; axis xyz; start (-0.01 0.01 0.1); end ( 0.01 0.01 0.1); } z_0.15m { type face; axis xyz; start (-0.01 0.01 0.15); end ( 0.01 0.01 0.15); } ); } |
以下のコマンドで結果をpostProcessingに出力。
1 | postProcess -latestTime -func sampleDict | tee log.sample |
壁面摩擦応力(wallShearStress）の出力
1 2 3 | simpleFoam simpleFoam -postProcess -func wallShearStress -noZero -noFunctionObjects postProcess -func sampleDict -noZero -latestTime |
- simpleFoam
計算実行 - simpleFoam -postProcess -func wallShearStress -noZero -noFunctionObjects
計算実行後に実行する。各時刻にwallShearStressを出力される。
-noZero：0ディレクトリは使わないし出力もしない
-noFunctionObjects：controlDictのfunctionObjectは使わない。このオプションを付けないとfunctionObjectの記述に従って出力されてしまう。
各時刻のディレクトリ内にwallShearStressがあればParaViewでも表示できる。 - postProcess -func sampleDict -noZero -latestTime
postProcessingにsampleDictができて上記のsampleDictのfieldsを出力してくれる。wallShearStressを書いておくと出力される。
※ただしwallShearStressはinterpolationSchemeをcellPointにしないと出力してくれませんでした。
1 | interpolationScheme cellPoint;//cell; |
interpolationScheme
- cell：Cell-centre value assumed constant over cell
- cellPoint：Linear weighted interpolation using cell values
yPlusの出力
上記同様のことをyPlusで行えばよい。
- 各時間でyPlusを出力
system/controlDict
1 2 3 4 | functions { #includeFunc yPlus } |
もしくは
1 | simpleFoam -postProcess -func yPlus -noZero -latestTime |
-latestTime：最後の結果のみ出力するオプション
これでyPlusをParaViewで確認できます。
- system/sampleDict
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 | sampleDict { type sets; libs (sampling); setFormat raw; interpolationScheme cell; fields (U p wallShearStress yPlus); writeControl writeTime; sets ( z_0.0 { type face; axis y; start (0 -0.005 0.001); end (0 0.005 0.001); } ); } |
「fields (U p wallShearStress yPlus);」にyPlusを追加。
1 | postProcess -func sampleDict -noZero -latestTime |
これでpostProcessing/postProcessing/最終時刻/z_0.0m_yPlus.xyに結果が出力されます。
1 2 3 4 5 6 | 0.0450069 23.1191 0.0452398 0 0.0454727 0 0.0457056 0 ... （以下省略） |
特定の場所のy+の具体的な値をを知りたい場合は便利です。
特定値の流速のみを表示(ParaView)
ParaViewで設定を変えます。
「Edit」>[Settings]から以下の項目にちゅえっくを入れます。
流速のコンターも描くことができます。
特定の流速の範囲に絞って表示することができます。
WSLでメモリ不足になったとき
WSLは何も設定しない場合、ホスト(Windows)メモリの 50%または8GBのどちらか少ない方になるように調整します。
Windows側が16GBのメモリの場合は以下のように、
1 | free -h |
で確認するとトータル7.7BGB使えるようになっており、あふれた分はスワップで2GBまで保持できるようになっています。
シミュレーションをしているとためにメモリが足らなくなるのでもう少し使えるように設定を変えます。
使えるメモリを80%までに引き上げて、残り足らなくなった分は少々遅くなっても良いのでスワップで保持するようにします。
.wslconfig は以下のパスに作成することで、WSL2 起動時に設定が読み込まれる仕組みになっています。
1 | C:\Users\[ユーザ名]\.wslconfig |
.wslconfigを以下とします。
.wslconfig
1 2 3 | [wsl2] memory=13GB swap=10GB |
PowerShellを管理者権限でシャットダウンします。
1 | wsl --shutdown |
環境変数が展開されない
例えば以下のように$FOAM_TUTORIALSの後にTABを押して、補完機能を使いながらパスをたどっていきたい場合があるとします。
1 | vi $FOAM_TUTORIALS/[TAB] |
このようにすると、バックスラッシュが先頭に入って面倒だったりします。
「shopt」で、シェルオプションの設定状況（on／off）を一覧表示します。「shopt シェルオプション名」で表示対象を指定できます。
1 | shopt -s direxpand |
これで環境変数が自動で展開されます。
1 2 3 4 5 | vi $FOAM_TUTORIALS/[TAB] ↓以下に展開される vi /opt/OpenFOAM/OpenFOAM-v2012/tutorials/ |
境界条件の参考
非圧縮流れの場合(simpleFoma、pimpleFoamなど）
| patch type | U | p | k | epsilon | nut | omega | |
| 流速指定 | patch | Region0 { type fixedValue; value uniform (1 0 0); } | Region0 { type zeroGradient; } | Region0 { type fixedValue; value uniform 0.001; } | Region0 { type fixedValue; value uniform 0.01; } | Region0 { type calculated; value uniform 0.001; } | Region1 { type fixedValue; value $internalField; } |
| 質量流 | patch | Region { type flowRateInletVelocity; massFlowRate constant 1; //[kg/s] rhoInlet 1000.0; //[kg/m3] } | Region { type zeroGradient; } | Region0 { type fixedValue; value uniform 0.001; } | Region { type fixedValue; value uniform 0.01; } | Region { type calculated; value uniform 0.001; } | Region1 { type fixedValue; value $internalField; } |
| 体積流量指定 | patch | Region0 { type flowRateInletVelocity; volumetricFlowRate constant 1; //[m3/s] } | Region0 { type zeroGradient; } | Region0 { type fixedValue; value uniform 0.001; } | Region0 { type fixedValue; value uniform 0.01; } | Region0 { type calculated; value uniform 0.001; } | Region1 { type fixedValue; value $internalField; } |
| 法線方向流速 | patch | Region { type surfaceNormalFixedValue; refValue uniform -1.0; } | Region { type zeroGradient; } | Region { type fixedValue; value uniform 0.001; } | Region { type fixedValue; value uniform 0.01; } | Region { type zeroGradient; } | Region1 { type zeroGradient; } |
| 静止圧指定 | patch | Region1 { type zeroGradient; } | Region1 { type fixedValue; value uniform 0.0; } | Region1 { type zeroGradient; } | Region1 { type zeroGradient; } | Region1 { type calculated; value uniform 0.001; } | Region1 { type fixedValue; value $internalField; } |
| 全圧指定 | patch | Region1 { type zeroGradient; } | Region1 { type totalPressure; p0 uniform 0.0; value uniform 0.0; } | Region1 { type zeroGradient; } | Region1 { type zeroGradient; } | Region1 { type calculated; value uniform 0.001; } | Region1 { type fixedValue; value $internalField; } |
| 自然流出流入 | patch | Region2 { type pressureInletOutletVelocity; value uniform (0 0 0); } | Region2 { type totalPressure; p0 uniform 0; value uniform 0; } | Region2 { type inletOutlet; inletValue uniform 0.001; value uniform 0.001; } | Region2 { type inletOutlet; inletValue uniform 0.001; value uniform 0.001; } | Region2 { type calculated; value uniform 0.001; } | Region1 { type zeroGradient; } |
| 静止壁 | wall | Region2 { type fixedValue; value uniform (0 0 0); } | Region2 { type zeroGradient; } | Region2 { type kqRWallFunction; value uniform 0.0; } | Region2 { type epsilonWallFunction; value uniform 0.0; } | Region2 { type nutkWallFunction; value uniform 0.0; } | Region1 { type omegaWallFunction; } |
| 回転壁 | type wall; inGroups 1(wall); | Region { type rotatingWallVelocity; origin (0 0 0); axis (0 0 1); omega 0.0;//[rad/s] } | Region { type zeroGradient; } | Region { type kqRWallFunction; value uniform 0.0; } | Region { type epsilonWallFunction; value uniform 0.0; } | Region { type nutkWallFunction; value uniform 0.0; } | |
| 対称面 | type symmetry; inGroups 1(symmetry); | Region1 { type symmetry; } | Region1 { type symmetry; } | Region1 { type symmetry; } | Region1 { type symmetry; } | Region1 { type symmetry; } | |
| 滑り条件 | patch |
熱流体の場合
| U | T | P_rgh | P | K | epsilon | nut | alphat | |
| 流速 | type fixedValue; value uniform (1 0 0); | type fixedValue; value uniform 273.15; | type fixedFluxPressure; value uniform 101325; | type calculated; value uniform 101325.0; | type fixedValue; value uniform 0.001; | type fixedValue; value uniform 0.01; | type calculated; value uniform 0.001; | type calculated; value uniform 0.0; |
| 質量流量 | type flowRateInletVelocity; massFlowRate constant 1; //[kg/s] rhoInlet -1.0; //[kg/m3] | type fixedValue; value uniform 273.15; | type fixedFluxPressure; value uniform 101325; | type calculated; value uniform 101325.0; | type fixedValue; value uniform 0.001; | type fixedValue; value uniform 0.01; | type calculated; value uniform 0.001; | type calculated; value uniform 0.0; |
| 体積流量 | type flowRateInletVelocity; volumetricFlowRate constant 1; //[m3/s] | type fixedValue; value uniform 273.15; | type fixedFluxPressure; value uniform 101325; | type calculated; value uniform 101325.0; | type fixedValue; value uniform 0.001; | type fixedValue; value uniform 0.01; | type calculated; value uniform 0.001; | type calculated; value uniform 0.0; |
| 法線方向流速 | type surfaceNormalFixedValue; refValue uniform -1.0; | type fixedValue; value uniform 273.15; | type fixedFluxPressure; value uniform 101325; | type calculated; value uniform 101325.0; | type fixedValue; value uniform 0.001; | type fixedValue; value uniform 0.01; | type calculated; value uniform 0.001; | type calculated; value uniform 0.0; |
| 静圧 | type zeroGradient; | type zeroGradient; | type fixedValue; value uniform 101325.0; | type calculated; value uniform 101325.0; | type zeroGradient; | type zeroGradient; | type calculated; value uniform 0.001; | type compressible::alphatWallFunction; value uniform 0.0; |
| 自然流入出 | type pressureInletOutletVelocity; value uniform (0 0 0); | type fixedValue; value uniform 273.15; | type totalPressure; p0 uniform 101325; value uniform 101325; | type calculated; value uniform 101325.0; | type inletOutlet; inletValue uniform 0.001; value uniform 0.001; | type inletOutlet; inletValue uniform 0.001; value uniform 0.001; | type calculated; value uniform 0.001; | type calculated; value uniform 0.0; |
| すべり壁（温度） | type slip; | type fixedValue; value uniform 273.15; | type fixedFluxPressure; value uniform 101325; | type calculated; value uniform 101325.0; | type zeroGradient; | type zeroGradient; | type zeroGradient; | type compressible::alphatWallFunction; value uniform 0.0; |
| すべり壁（断熱） | type slip; | type zeroGradient; | type fixedFluxPressure; value uniform 101325; | type calculated; value uniform 101325.0; | type zeroGradient; | type zeroGradient; | type zeroGradient; | type compressible::alphatWallFunction; value uniform 0.0; |
| 回転壁（温度） | type rotatingWallVelocity; origin (0 0 0); axis (0 0 1); omega 0.0;//[rad/s] | type fixedValue; value uniform 273.15; | type fixedFluxPressure; value uniform 101325; | type calculated; value uniform 101325.0; | type kqRWallFunction; value uniform 0.0; | type epsilonWallFunction; value uniform 0.0; | type calculated; value uniform 0.001; | type compressible::alphatWallFunction; value uniform 0.0; |
| 全圧 | type zeroGradient; | type zeroGradient; | type totalPressure; p0 uniform 101325.0; value uniform 101325.0; | type calculated; value uniform 101325.0; | type zeroGradient; | type zeroGradient; | type calculated; value uniform 0.001; | type compressible::alphatWallFunction; value uniform 0.0; |
モデルの形状線を表示
例えばFreeCADなどでmm単位でモデル作成したとします。
それをmodelというフォルダにmodel.stlとして保存した場合を考えます。
FreeCADは数値しか拾っていないので、SI単位系で流体解析を行う場合には0.001倍にしてm単位として計算する必要があります。
- model/model.stl：mm単位
- model/model_m.stl：m単位
このようにしたいとします。
surfaceConvertでスケール変換し、それを特徴線として出力(eMesh)、さらにParaViewで表示するためにobjファイルにするコマンドは以下です。
1 2 3 4 | surfaceConvert -scale 0.001 model/model.stl model/model_m.stl cp -r model/model_m.stl constant/triSurface/ surfaceFeatureExtract surfaceFeatureConvert constant/triSurface/model_m.eMesh constant/triSurface/model_m.obj |
特徴線は以下のように設定しておきます。
system/surfaceFeatureExtractDict
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 | model_m.stl { // How to obtain raw features (extractFromFile || extractFromSurface) extractionMethod extractFromSurface; // Mark edges whose adjacent surface normals are at an angle less // than includedAngle as features // - 0 : selects no edges // - 180: selects all edges includedAngle 150; subsetFeatures { // Keep nonManifold edges (edges with >2 connected faces) nonManifoldEdges no; // Keep open edges (edges with 1 connected face) openEdges yes; } // Write options // Write features to obj format for postprocessing writeObj yes; } |
model_m.objはParaViewで読み込めば特徴線も読み込めます。
モデルのスケール変換
OpenFOAMのメッシュを作成した後に節点座標を「移動・回転・スケール」させることができます。
以下の例は、スケールを1/1000倍にした例です。
1 2 | # 節点の場合 transformPoints -scale "(0.001 0.0001 0.001)" |
また、stlを変換することもできます。
1 2 | # stl変換の場合 surfaceConvert -scale 0.001 constant/triSurface/model.stl constant/triSurface/model_m.stl |
クーラン数の確認
クーラン数は$Co = \frac{u\Delta t}{\Delta x}$で定義されます。
陰解法の場合は、安定性の観点からクーラン数を必ずしも1以下にする必要はなく、$\Delta t$があまりにも小さくなりすぎて計算が終わらない場合は、精度を多少犠牲にしてクーラン数を大きくとり$Delta t$が小さくなりすぎないようにします。
しかし、クーラン数の最大値を指定できますが、最小値は指定できません（たぶん）
なので、1度過去に計算したもののクーラン数を確認したい場合があります。
以下のコマンドで各時刻（ステップ）ごとのクーラン数を出力します。
1 | postProcess -func CourantNo |
これでParaViewでのクーラン数の分布を見ることもできます。
クーラン数の最大、最小を見たい場合は以下のようにsystem/controlDictのfunctionsに記述して、fieldの最大最小を出力しておきます。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 | fieldMinMax { type fieldMinMax; libs ( "libfieldFunctionObjects.so" ); region region0; writeToFile yes; log yes; location yes; mode magnitude; fields (U p Co ); } |
system/controlDictのfunctionsを書きます。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 | CourantNo1 { // Mandatory entries (unmodifiable) type CourantNo; libs (fieldFunctionObjects); // Optional entries (runtime modifiable) rho rho; // Optional (inherited) entries field phi; result Co;//CourantNumberField; region region0; enabled true; log true; timeStart 0; timeEnd 1000; executeControl timeStep; executeInterval 1; writeControl writeTime;;//timeStep; writeInterval 1; } |
記述後にpostProcess のオプション付きで計算を実行します（出力計算のみ）。
1 | interFoam -postProcess |
postProcessingに出力されるので、グラフにするなりすればOK。
例えば、以下のPythonプログラムなど。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 | import pandas as pd import matplotlib.pyplot as plt dir_ = 20 df_Co = pd.read_table(f'./postProcessing/fieldMinMax/{dir_}/fieldMinMax.dat',skiprows=6, sep="\s+", header=None) df_MaxCo = df_Co[[0,1, 6]] df_MaxCo.columns = ["time", "filed", "MaxCo"] df_MaxCo = df_MaxCo[df_MaxCo["filed"] == "Co"] plt.plot(df_MaxCo["time"], df_MaxCo["MaxCo"]) plt.ylim([0,12]) plt.xlabel("time[sec]") plt.ylabel("MaxCo") plt.legend(loc='best',bbox_to_anchor=(1, 1)) plt.title(f'MaxCo') plt.grid() |
参考書
PENGUINさんサイトを体系的に学べる書籍となっています。ネット記事でも十分勉強できるのですが、OpenFOAMの初学者でOpenFOAMをインストール済みであれば一冊持って置き、体系的に学ぶのが良いでしょう。
あとは初心者向けに丁寧に解説がされているこちらの書籍もお勧めです。最後の章にはオーバーセットメッシュ（重合メッシュ）の機能を使った解析を最後まで丁寧に解説しているので挫折することはないでしょう。
辞書的な参考書としては、こちらが良いです。
Foundation版ですが、大変参考になります。
Follow @t_kun_kamakiri