LoginSignup
4
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

OpenFOAMAdvent Calendar 2017

Day 8
@TatsuyaKatayama

OpenMDAOでOpenFOAMのTjunctionの流出流量を等しくする

Last updated at Posted at 2017-12-08
1 / 9

#やりたいこと

  • 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

下準備

  1. TJunctionチュートリアルのコピー

    bash
     cp -r $FOAM_TUTORIALS/incompressible/pimpleFoam/TJunction .

  2. 境界条件の変更

    • 実はチュートリアルそのままでは都合が悪い(最適解を見つけるより先にメッシュが歪む)
    • 都合よい問題になるように境界条件を変更する
      • 非定常解析→定常解析
      • 境界条件の変更
        • 圧力流入→強制流入　
        • outlet*圧力の変更

##変更箇所

  • 非定常解析→定常解析(systemディレクトリ内(blockMeshDict除く))

    diff
    diff --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ディレクトリ内)

    diff
    diff -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で最適化する方がいい気がするが気にしない

※OpenMDAOのドキュメント

最適化結果

以下のスクリプトにより結果を可視化

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")

opt.png

赤色が今回の目的関数で青色が変数。無事、目的関数が最小化されている。

尚、圧力コンターと流速コンターは下記の通りである

変形量も小さくインパクトのない結果になってしまった。。。

4
2
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
4
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?