やりたいこと
- OpenFOAMのTjunctionチュートリアルで流出流量を均等に分流する
- outlet1,outlet2ともに2.5 [m/s]を目指す（inlet 5 [m/s]で流入）
- outlet2の方が圧力が低く流れやすいので傾斜をつけてoutlet1側の流量を増やす
\begin{align}
{\rm min} \: \: \:& (2.5 -U_{outlet2})^2 \\
{\rm subject \: to} \: \: \:& -0.025\leq h \leq 0.025
\end{align}
各種バージョン情報
はじめに今回の最適化を実行する環境は以下の通り
- OpenFOAM==4.1
- Python==3.6.0
- PyFoam==0.6.6
- OpenMDAO==1.7.3
下準備
-
TJunctionチュートリアルのコピー
bashcp -r $FOAM_TUTORIALS/incompressible/pimpleFoam/TJunction . -
境界条件の変更
- 実はチュートリアルそのままでは都合が悪い(最適解を見つけるより先にメッシュが歪む)
- 都合よい問題になるように境界条件を変更する
- 非定常解析→定常解析
- 境界条件の変更
- 圧力流入→強制流入
- outlet*圧力の変更
変更箇所
-
非定常解析→定常解析(systemディレクトリ内(blockMeshDict除く))
diffdiff --exclude blockMeshDict -u ../TJunction.org/system/controlDict system/controlDict --- ../TJunction.org/system/controlDict 2017-06-02 20:13:23.764228447 +0900 +++ system/controlDict 2017-06-03 10:19:41.519505802 +0900 @@ -10,26 +10,25 @@ version 2.0; format ascii; class dictionary; - location "system"; object controlDict; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -application pimpleFoam; +application simpleFoam; -startFrom startTime; +startFrom latestTime; startTime 0; stopAt endTime; -endTime 1.5; +endTime 500; -deltaT 0.001; +deltaT 1; -writeControl adjustableRunTime; +writeControl timeStep; -writeInterval 0.1; +writeInterval 500; purgeWrite 0; @@ -37,7 +36,7 @@ writePrecision 6; -writeCompression uncompressed; +writeCompression compressed; timeFormat general; @@ -45,40 +44,9 @@ runTimeModifiable true; -adjustTimeStep yes; - -maxCo 5; - functions { - probes - { - // Where to load it from - libs ( "libsampling.so" ); - - type probes; - - // Name of the directory for probe data - name probes; - - // Write at same frequency as fields - writeControl writeTime; - writeInterval 1; - - // Fields to be probed - fields - ( - p U - ); - - probeLocations - ( - ( 1e-06 0 0.01 ) // at inlet - (0.21 -0.20999 0.01) // at outlet1 - (0.21 0.20999 0.01) // at outlet2 - (0.21 0 0.01) // at central block - ); - } } + // ************************************************************************* // diff --exclude blockMeshDict -u ../TJunction.org/system/fvSchemes system/fvSchemes --- ../TJunction.org/system/fvSchemes 2017-06-02 20:13:23.760228447 +0900 +++ system/fvSchemes 2017-06-03 10:19:41.519505802 +0900 @@ -17,7 +17,7 @@ ddtSchemes { - default Euler; + default steadyState; } gradSchemes @@ -28,13 +28,13 @@ divSchemes { default none; - div(phi,U) Gauss limitedLinearV 1; - div(phi,k) Gauss limitedLinear 1; - div(phi,epsilon) Gauss limitedLinear 1; - div(phi,R) Gauss limitedLinear 1; - div(R) Gauss linear; - div(phi,nuTilda) Gauss limitedLinear 1; + div(phi,U) bounded Gauss linearUpwind grad(U); + div(phi,k) bounded Gauss limitedLinear 1; + div(phi,epsilon) bounded Gauss limitedLinear 1; + div(phi,omega) bounded Gauss limitedLinear 1; + div(phi,v2) bounded Gauss limitedLinear 1; div((nuEff*dev2(T(grad(U))))) Gauss linear; + div(nonlinearStress) Gauss linear; } laplacianSchemes @@ -52,5 +52,10 @@ default corrected; } +wallDist +{ + method meshWave; +} + // ************************************************************************* // diff --exclude blockMeshDict -u ../TJunction.org/system/fvSolution system/fvSolution --- ../TJunction.org/system/fvSolution 2017-06-02 20:13:23.760228447 +0900 +++ system/fvSolution 2017-06-03 10:19:41.519505802 +0900 @@ -21,50 +21,38 @@ { solver GAMG; tolerance 1e-06; - relTol 0.01; - smoother GaussSeidel; - } - - pFinal - { - solver GAMG; - tolerance 1e-06; - relTol 0; + relTol 0.1; smoother GaussSeidel; } - "(U|k|epsilon)" + "(U|k|epsilon|omega|f|v2)" { solver smoothSolver; smoother symGaussSeidel; tolerance 1e-05; relTol 0.1; } - - "(U|k|epsilon)Final" - { - $U; - tolerance 1e-05; - relTol 0; - } } -PIMPLE +SIMPLE { - nOuterCorrectors 1; - nCorrectors 2; nNonOrthogonalCorrectors 0; - pRefCell 0; - pRefValue 0; + consistent yes; + + residualControl + { + p 1e-2; + U 1e-3; + "(k|epsilon|omega|f|v2)" 1e-3; + } } relaxationFactors { equations { - "U.*" 1; - "k.*" 1; - "epsilon.*" 1; + U 0.9; // 0.9 is more stable but 0.95 more convergent + ".*" 0.9; // 0.9 is more stable but 0.95 more convergent } } -
境界条件の変更(0ディレクトリ内)
diffdiff -u ../TJunction.org/0/p 0/p --- ../TJunction.org/0/p 2017-06-02 20:13:23.744228447 +0900 +++ 0/p 2017-06-03 10:15:21.253437790 +0900 @@ -22,18 +22,13 @@ { inlet { - type uniformTotalPressure; - p0 table - ( - (0 10) - (1 40) - ); + type zeroGradient; } outlet1 { type fixedValue; - value uniform 10; + value uniform 2; } outlet2 diff -u ../TJunction.org/0/U 0/U --- ../TJunction.org/0/U 2017-06-02 20:13:23.752228447 +0900 +++ 0/U 2017-06-03 10:15:21.253437790 +0900 @@ -22,7 +22,9 @@ { inlet { - type pressureInletOutletVelocity; + + type flowRateInletVelocity; + volumetricFlowRate constant 0.002; // [m3/s] value uniform (0 0 0); }
最適化計算
OpenFOAMの解析をPyFoamで自動化し、OpenMDAOの最適化エンジンで最適化する
まずはクラスを作成する
Python
import numpy as np
import shlex
from PyFoam.RunDictionary.SolutionDirectory import SolutionDirectory
from PyFoam.RunDictionary.ParsedParameterFile import ParsedParameterFile
from PyFoam.Basics.DataStructures import Vector
from PyFoam.Execution.BasicRunner import BasicRunner
from PyFoam.LogAnalysis.UtilityAnalyzer import UtilityAnalyzer
from PyFoam.LogAnalysis.LogAnalyzerApplication import LogAnalyzerApplication
from openmdao.api import Component
class TJunctionFlow(Component):
""" discharged flow velocity in upper pipe of TJunction """
def __init__(self):
super(TJunctionFlow, self).__init__()
self.add_param('h', val=0.0)
self.add_output('dVel', shape=1)
def solve_nonlinear(self, params, unknowns, resids):
"""
Make the flow velocities of the two ports
with different pressures the same
target flow velocities are 2.5 [m/s].
inlet 5.0 [m/s]
outlet1 2.5 [m/s]
outlet2 2.5 [m/s]
"""
h = params['h']
#下記に定義してあるcalc_TJunctionは速度ベクトルを返す
#今回はUyが必要
vel = self.calc_TJunction(h)[1]
unknowns['dVel'] = np.square(2.5 - vel)
def calc_TJunction(self, h):
"""
Calc TJunction Flown of OpenFOAM.
input h is inlet hight.
Static pressure at outlet1 is 2.
Static pressure at outlet2 is 0.
"""
TJunction = SolutionDirectory("./TJunction")
#結果の削除
TJunction.clearResults()
#PyFoamの出力を削除
TJunction.clearPattern("PyFoam*")
blockMeshDict = ParsedParameterFile(TJunction.blockMesh())
blockMeshDict.content["vertices"][0] = Vector(0, -0.01+h, 0)
blockMeshDict.content["vertices"][3] = Vector(0, 0.01+h, 0)
blockMeshDict.content["vertices"][10] = Vector(0, -0.01+h, 0.02)
blockMeshDict.content["vertices"][13] = Vector(0, 0.01+h, 0.02)
blockMeshDict.writeFile()
#blockMesh -case path_to_case"コマンドを分割
blockMeshCMD = shlex.split('blockMesh -case %s' % TJunction.name )
#blockMeshを実行を行うインスタンスの生成
blockMeshRunner = BasicRunner(blockMeshCMD, silent=True)
#blockMeshの実行。実行結果情報を返す
blockMeshState = blockMeshRunner.start()
#"simpleFoam -case path_to_case"コマンドを分割
simpleFoamCMD = shlex.split('simpleFoam -case %s' % TJunction.name )
#simpleFoamを実行を行うインスタンスの生成
simpleFoamRunner = BasicRunner(simpleFoamCMD, silent=True)
#blockMeshの実行。実行結果情報を返す
simpleFoamState = simpleFoamRunner.start()
#"postProcess -func "patchAverage(name=outlet2,U)" -latestTime -case path_to_case"コマンドを分割
postProcessCMD = shlex.split(
'postProcess -func \"patchAverage(name=outlet2,U)\" -latestTime -case %s' % TJunction.name, posix=False )
#postProcessを実行を行うインスタンスの生成
postProcessRunner = BasicRunner(postProcessCMD, silent=True)
#postProcessの実行。実行結果情報を返す
postProcessState = postProcessRunner.start()
# logの分析
ua = UtilityAnalyzer()
#logファイル内のピックアップする正規表現
ua.addExpression("U_outlet2",
"average\(outlet2\) of U = \((%f%) (%f%) (%f%)\)")
logApp = LogAnalyzerApplication(ua)
logApp.run(postProcessState["logfile"])
#output2のU(ベクトル量)を返す
return np.array(ua.getData("U_outlet2"),dtype=float)
OpenMDAOのComponentの派生クラスを作成している。
勾配が計算できる場合linearizeメソッドでヤコビ行列を計算させるが、
今回は変数が少ないので差分法でOpenMDAOに頑張ってもらう。
クラスの定義ができたら実行を行う
Python
from __future__ import print_function
from openmdao.api import IndepVarComp, Component, Problem, Group, SqliteRecorder
from openmdao.api import ScipyOptimizer
top = Problem()
root = top.root = Group()
root.add('p1', IndepVarComp('h', 0.0))
root.add('TJunc', TJunctionFlow())
root.connect('p1.h', 'TJunc.h')
#差分法による勾配計算のための設定
root.TJunc.deriv_options["type"] = "fd"
root.TJunc.deriv_options["form"] = "central"
root.TJunc.deriv_options["step_size"] = 1.0e-4
#実際の最適化エンジンはscipyのSLSQP
top.driver = ScipyOptimizer()
top.driver.options['optimizer'] = 'SLSQP'
top.driver.add_desvar('p1.h', lower=-0.025, upper=0.025)
top.driver.add_objective('TJunc.dVel')
#計算過程のレコーダーの設定
recorder = SqliteRecorder('TJunc')
recorder.options['record_params'] = True
recorder.options['record_metadata'] = True
top.driver.add_recorder(recorder)
top.setup()
top.run()
top.cleanup()
CAEの最適化は、差分法を用いてわざわざ勾配を計算して勾配法を用いるより
勾配を必要としないGAで最適化する方がいい気がするが気にしない
最適化結果
以下のスクリプトにより結果を可視化
Python
from matplotlib import pyplot as plt
import pandas as pd
import sqlitedict
db =sqlitedict.SqliteDict("TJunc","iterations")
df = pd.DataFrame()
for k in db.keys():
a = pd.DataFrame.from_dict(db[k]["Unknowns"],orient="index")
df = df.append(a.T)
df = df.reset_index(drop=True)
df.plot(grid=True)
plt.savefig("opt.png")
赤色が今回の目的関数で青色が変数。無事、目的関数が最小化されている。
尚、圧力コンターと流速コンターは下記の通りである
変形量も小さくインパクトのない結果になってしまった。。。