PR

CMPLXFOILで翼型の最適化を行う

CMPLXFOILを使って2次元翼型の最適化を行う

スポンサーリンク

はじめに

CMPLXFOILを使って2次元翼型の最適化を行う

↓公式ドキュメント

Optimization Tutorial — CMPLXFOIL documentation

↓インストール

↓XFLRの使い方について

前回デスクトップに作成したCMPLXFOILというフォルダの中に「naca4412」というフォルダを作り,その中に必要なデータを入れていく

:~$ mkdir -p /mnt/c/Users/mtkbirdman/Desktop/CMPLXFOIL/naca4412

NACA4412の準備

XFLR5を開き「Module>Direct Foil Design」をクリックする

上のメニューから「Foil>Naca Foil」をクリックする

「4 or 5 digits: 4412」「Number of Panels: 100」と入力して「OK」をクリックする

NACA4412ができた

上のメニューから「Foil>Export」で翼を出力する

「C:\Users\mtkbirdman\Desktop\CMPLXFOIL\naca4412」に保存する

最後に1行目の翼型名を削除すればいい

 1.0000000     0.0012600
 0.9928006     0.0032198
 0.9798895     0.0066835
 0.9635184     0.0109823
 0.9445469     0.0158354
 0.9234998     0.0210604
 0.9007495     0.0265227
 0.8765830     0.0321168
 0.8512328     0.0377572
 0.8248944     0.0433730
 0.7977362     0.0489045
 0.7699063     0.0543008
 0.7415373     0.0595184
 0.7127487     0.0645195
 0.6836500     0.0692712
 0.6543419     0.0737448
 0.6249182     0.0779148
 0.5954661     0.0817592
 0.5660677     0.0852584
 0.5368005     0.0883955
 0.5077378     0.0911557
 0.4789491     0.0935267
 0.4505011     0.0954983
 0.4224571     0.0970626
 0.3948782     0.0982103
 0.3678229     0.0988053
 0.3413476     0.0987893
 0.3155065     0.0981781
 0.2903524     0.0969911
 0.2659361     0.0952510
 0.2423068     0.0929845
 0.2195125     0.0902210
 0.1975998     0.0869933
 0.1766138     0.0833369
 0.1565988     0.0792897
 0.1375977     0.0748918
 0.1196525     0.0701852
 0.1028043     0.0652129
 0.0870931     0.0600189
 0.0725582     0.0546478
 0.0592379     0.0491435
 0.0471700     0.0435495
 0.0363913     0.0379076
 0.0269381     0.0322577
 0.0188458     0.0266369
 0.0121494     0.0210790
 0.0068832     0.0156139
 0.0030809     0.0102667
 0.0007756     0.0050573
 0.0000000     0.0000000
 0.0007756    -0.0047474
 0.0030809    -0.0090391
 0.0068832    -0.0128843
 0.0121494    -0.0162931
 0.0188458    -0.0192762
 0.0269381    -0.0218453
 0.0363913    -0.0240133
 0.0471700    -0.0257940
 0.0592379    -0.0272029
 0.0725582    -0.0282569
 0.0870931    -0.0289743
 0.1028043    -0.0293755
 0.1196525    -0.0294825
 0.1375977    -0.0293193
 0.1565988    -0.0289118
 0.1766138    -0.0282876
 0.1975998    -0.0274762
 0.2195125    -0.0265089
 0.2423068    -0.0254180
 0.2659361    -0.0242376
 0.2903524    -0.0230023
 0.3155065    -0.0217477
 0.3413476    -0.0205094
 0.3678229    -0.0193230
 0.3948782    -0.0182234
 0.4224571    -0.0171746
 0.4505011    -0.0160650
 0.4789491    -0.0149118
 0.5077378    -0.0137351
 0.5368005    -0.0125542
 0.5660677    -0.0113870
 0.5954661    -0.0102497
 0.6249182    -0.0091567
 0.6543419    -0.0081203
 0.6836500    -0.0071506
 0.7127487    -0.0062555
 0.7415373    -0.0054401
 0.7699063    -0.0047076
 0.7977362    -0.0040587
 0.8248944    -0.0034919
 0.8512328    -0.0030041
 0.8765830    -0.0025904
 0.9007495    -0.0022449
 0.9234998    -0.0019608
 0.9445469    -0.0017313
 0.9635184    -0.0015496
 0.9798895    -0.0014106
 0.9928006    -0.0013115
 1.0000000    -0.0012600

CMPLXFOILのプログラムの準備

Ubuntuを起動する

以下のコマンドを実行する

:~$ . /mnt/c/Users/mtkbirdman/Desktop/CMPLXFOIL/.py310venv/bin/activate
(.py310venv) :~$ cd /mnt/c/Users/mtkbirdman/Desktop/CMPLXFOIL/naca4412
(.py310venv) :/mnt/c/Users/mtkbirdman/Desktop/CMPLXFOIL/naca4412$

ソースコード

ソースコードはこれ

# ======================================================================
#         Import modules
# ======================================================================
# rst imports (beg)
import os
import numpy as np
from mpi4py import MPI
from baseclasses import AeroProblem
from pygeo import DVConstraints, DVGeometryCST
from pyoptsparse import Optimization, OPT
from multipoint import multiPointSparse
from cmplxfoil import CMPLXFOIL, AnimateAirfoilOpt

# rst imports (end)

# ======================================================================
#         Specify parameters for optimization
# ======================================================================
# rst params (beg)
mycl = 1.0  # lift coefficient constraint
alpha = 5.0 # initial angle of attack (zero if the target cl is zero)
mach = 0.0  # Mach number
Re = 5e5  # Reynolds number
T = 288.15  # 1976 US Standard Atmosphere temperature @ sea level (K)
# rst params (end)

airfoilname = "naca4412.dat"

# ======================================================================
#         Create multipoint communication object
# ======================================================================
# rst procs (beg)
MP = multiPointSparse(MPI.COMM_WORLD)
MP.addProcessorSet("cruise", nMembers=1, memberSizes=MPI.COMM_WORLD.size)
MP.createCommunicators()
# rst procs (end)

# ======================================================================
#         Create output directory
# ======================================================================
# rst dir (beg)
curDir = os.path.abspath(os.path.dirname(__file__))
outputDir = os.path.join(curDir, "output")

if not os.path.exists(outputDir):
    os.mkdir(outputDir)
# rst dir (end)

# ======================================================================
#         CFD solver set-up
# ======================================================================
# rst solver (beg)
aeroOptions = {
    "writeSolution": True,
    "writeSliceFile": True,
    "writeCoordinates": True,
    "plotAirfoil": True,
    "outputDirectory": outputDir,
}

# Create solver
CFDSolver = CMPLXFOIL(os.path.join(curDir, airfoilname), options=aeroOptions)
# rst solver (end)

# ======================================================================
#         Set up flow conditions with AeroProblem
# ======================================================================
# rst ap (beg)
ap = AeroProblem(
    name="fc",
    alpha=alpha if mycl != 0.0 else 0.0,
    mach=mach,
    reynolds=Re,
    reynoldsLength=1.0,
    T=T,
    areaRef=1.0,
    chordRef=1.0,
    evalFuncs=["cl", "cd"],
)
# Add angle of attack variable
if mycl != 0.0:
    ap.addDV("alpha", value=ap.alpha, lower=-10.0, upper=10.0, scale=1.0)
# rst ap (end)

# ======================================================================
#         Geometric Design Variable Set-up
# ======================================================================
# rst geom (beg)
nCoeff = 4  # number of CST coefficients on each surface
DVGeo = DVGeometryCST(os.path.join(curDir, airfoilname), numCST=nCoeff)

DVGeo.addDV("upper_shape", dvType="upper", lowerBound=-0.1, upperBound=0.5)
DVGeo.addDV("lower_shape", dvType="lower", lowerBound=-0.5, upperBound=0.1)

# Add DVGeo object to CFD solver
CFDSolver.setDVGeo(DVGeo)
# rst geom (end)

# ======================================================================
#         DVConstraint Setup
# ======================================================================
# rst cons (beg)
DVCon = DVConstraints()
DVCon.setDVGeo(DVGeo)
DVCon.setSurface(CFDSolver.getTriangulatedMeshSurface())

# Thickness, volume, and leading edge radius constraints
le = 0.0001
wingtipSpacing = 0.1
leList = [[le, 0, wingtipSpacing], [le, 0, 1.0 - wingtipSpacing]]
teList = [[1.0 - le, 0, wingtipSpacing], [1.0 - le, 0, 1.0 - wingtipSpacing]]

DVCon.addVolumeConstraint(
        leList=leList,
        teList=teList,
        nSpan=2,
        nChord=100,
        lower=1.00,
        upper=3.00,
        scaled=True,
)

DVCon.addThicknessConstraints2D(
        leList=leList,
        teList=teList,
        nSpan=2,
        nChord=100,
        lower=1.00,
        upper=3.00,
        scaled=True,
)

le = 0.01
leList = [[le, 0, wingtipSpacing], [le, 0, 1.0 - wingtipSpacing]]

DVCon.addLERadiusConstraints(
        leList=leList,
        nSpan=2,
        axis=[0,1,0],
        chordDir=[-1,0,0],
        lower=1.00,
        upper=3.00,
        scaled=True,
)

fileName = os.path.join(outputDir, "constraints.dat")
DVCon.writeTecplot(fileName)
# rst cons (end)

# ======================================================================
#         Functions:
# ======================================================================
# rst funcs (beg)
def cruiseFuncs(x):
    print(x)
    # Set design vars
    DVGeo.setDesignVars(x)
    ap.setDesignVars(x)
    # Run CFD
    CFDSolver(ap)
    # Evaluate functions
    funcs = {}
    DVCon.evalFunctions(funcs)
    CFDSolver.evalFunctions(ap, funcs)
    CFDSolver.checkSolutionFailure(ap, funcs)
    if MPI.COMM_WORLD.rank == 0:
        print("functions:")
        for key, val in funcs.items():
            if key == "DVCon1_thickness_constraints_0":
                continue
            print(f"    {key}: {val}")
    return funcs


def cruiseFuncsSens(x, funcs):
    funcsSens = {}
    DVCon.evalFunctionsSens(funcsSens)
    CFDSolver.evalFunctionsSens(ap, funcsSens)
    CFDSolver.checkAdjointFailure(ap, funcsSens)
    print("function sensitivities:")
    evalFunc = ["fc_cd", "fc_cl", "fail"]
    for var in evalFunc:
        print(f"    {var}: {funcsSens[var]}")
    return funcsSens


def objCon(funcs, printOK):
    # Assemble the objective and any additional constraints:
    funcs["obj"] = funcs[ap["cd"]]
    funcs["cl_con_" + ap.name] = funcs[ap["cl"]] - mycl
    if printOK:
        print("funcs in obj:", funcs)
    return funcs


# rst funcs (end)

# ======================================================================
#         Optimization Problem Set-up
# ======================================================================
# rst optprob (beg)
# Create optimization problem
optProb = Optimization("opt", MP.obj)

# Add objective
optProb.addObj("obj", scale=1e4)

# Add variables from the AeroProblem
ap.addVariablesPyOpt(optProb)

# Add DVGeo variables
DVGeo.addVariablesPyOpt(optProb)

# Add constraints
DVCon.addConstraintsPyOpt(optProb)

# Add cl constraint
optProb.addCon("cl_con_" + ap.name, lower=0.0, upper=0.0, scale=1.0)

# Enforce first upper and lower CST coefficients to add to zero
# to maintain continuity at the leading edge
jac = np.zeros((1, nCoeff), dtype=float)
jac[0, 0] = 1.0
optProb.addCon(
    "first_cst_coeff_match",
    lower=0.0,
    upper=0.0,
    linear=True,
    wrt=["upper_shape", "lower_shape"],
    jac={"upper_shape": jac, "lower_shape": jac},
)

# The MP object needs the 'obj' and 'sens' function for each proc set,
# the optimization problem and what the objcon function is:
MP.setProcSetObjFunc("cruise", cruiseFuncs)
MP.setProcSetSensFunc("cruise", cruiseFuncsSens)
MP.setObjCon(objCon)
MP.setOptProb(optProb)
optProb.printSparsity()
optProb.getDVConIndex()
# rst optprob (end)

# rst opt (beg)
# Run optimization
optOptions = {"IFILE": os.path.join(outputDir, "SLSQP.out")}
opt = OPT("SLSQP", options=optOptions)
sol = opt(optProb, MP.sens, storeHistory=os.path.join(outputDir, "opt.hst"))
if MPI.COMM_WORLD.rank == 0:
    print(sol)
# rst opt (end)

# ======================================================================
#         Postprocessing
# ======================================================================
# rst postprocessing (beg)
# Save the final figure
CFDSolver.airfoilAxs[1].legend(["Original", "Optimized"], labelcolor="linecolor")
CFDSolver.airfoilFig.savefig(os.path.join(outputDir, "OptFoil.pdf"))

# Animate the optimization
AnimateAirfoilOpt(outputDir, "fc").animate(
    outputFileName=os.path.join(outputDir, "OptFoil"), fps=10, dpi=300, extra_args=["-vcodec", "libx264"]
)
# rst postprocessing (end)
# ======================================================================
#         Import modules
# ======================================================================
# rst imports (beg)
import os
import numpy as np
import argparse
import ast
from mpi4py import MPI
from baseclasses import AeroProblem
from pygeo import DVConstraints, DVGeometryCST
from pyoptsparse import Optimization, OPT
from multipoint import multiPointSparse
from cmplxfoil import CMPLXFOIL, AnimateAirfoilOpt

# rst imports (end)
parser = argparse.ArgumentParser()
parser.add_argument("--output", type=str, default="output")
parser.add_argument("--opt", type=str, default="SLSQP", choices=["SLSQP", "SNOPT"])
parser.add_argument("--gridFile", type=str, default="n0012.cgns")
parser.add_argument("--optOptions", type=ast.literal_eval, default={}, help="additional optimizer options to be added")
args = parser.parse_args()
# ======================================================================
#         Specify parameters for optimization
# ======================================================================
# rst params (beg)
mycl = [0.6,1.0]  # lift coefficient constraint
alpha = [2.0,5.0] # initial angle of attack (zero if the target cl is zero)
mach = [0.0,0.0]  # Mach number
Re = [5e5,5e5]  # Reynolds number
T = [288.15,288.15]  # 1976 US Standard Atmosphere temperature @ sea level (K)
nFlowCases = len(Re)
# rst params (end)

airfoilname = "naca4412.dat"

# ======================================================================
#         Create multipoint communication object
# ======================================================================
# rst procs (beg)
# assign number of processors
nGroup = 1
nProcPerGroup = MPI.COMM_WORLD.size
MP = multiPointSparse(MPI.COMM_WORLD)
MP.addProcessorSet("cruise", nMembers=nGroup, memberSizes=nProcPerGroup)
comm, setComm, setFlags, groupFlags, ptID = MP.createCommunicators()
if not os.path.exists(args.output):
    if comm.rank == 0:
        os.mkdir(args.output)
# rst procs (end)

# ======================================================================
#         Create output directory
# ======================================================================
# rst dir (beg)
curDir = os.path.abspath(os.path.dirname(__file__))
outputDir = os.path.join(curDir, "output")

if not os.path.exists(outputDir):
    os.mkdir(outputDir)
# rst dir (end)

# ======================================================================
#         CFD solver set-up
# ======================================================================
# rst solver (beg)
aeroOptions = {
    "writeSolution": True,
    "writeSliceFile": True,
    "writeCoordinates": True,
    "plotAirfoil": True,
    "outputDirectory": outputDir,
}

# Create solver
CFDSolver = CMPLXFOIL(os.path.join(curDir, airfoilname), options=aeroOptions)
# rst solver (end)

# ======================================================================
#         Set up flow conditions with AeroProblem
# ======================================================================
# rst ap (beg)
aeroProblems = []
for i in range(nFlowCases):
    ap = AeroProblem(
        name="fc%d" % i,
        alpha=alpha[i],
        mach=mach[i],
        reynolds=Re[i],
        reynoldsLength=1.0,
        T=T[i],
        areaRef=1.0,
        chordRef=1.0,
        evalFuncs=["cl", "cd"],
    )
    # Add angle of attack variable
    ap.addDV("alpha", value=alpha[i], lower=0, upper=10.0, scale=1.0)
    aeroProblems.append(ap)
# rst ap (end)

# ======================================================================
#         Geometric Design Variable Set-up
# ======================================================================
# rst geom (beg)
nCoeff = 4  # number of CST coefficients on each surface
DVGeo = DVGeometryCST(os.path.join(curDir, airfoilname), numCST=nCoeff)

DVGeo.addDV("upper_shape", dvType="upper", lowerBound=-0.1, upperBound=0.5)
DVGeo.addDV("lower_shape", dvType="lower", lowerBound=-0.5, upperBound=0.1)

# Add DVGeo object to CFD solver
CFDSolver.setDVGeo(DVGeo)
# rst geom (end)

# ======================================================================
#         DVConstraint Setup
# ======================================================================
# rst cons (beg)
DVCon = DVConstraints()
DVCon.setDVGeo(DVGeo)
DVCon.setSurface(CFDSolver.getTriangulatedMeshSurface())

# Thickness, volume, and leading edge radius constraints
le = 0.0001
wingtipSpacing = 0.1
leList = [[le, 0, wingtipSpacing], [le, 0, 1.0 - wingtipSpacing]]
teList = [[1.0 - le, 0, wingtipSpacing], [1.0 - le, 0, 1.0 - wingtipSpacing]]

DVCon.addVolumeConstraint(
        leList=leList,
        teList=teList,
        nSpan=2,
        nChord=100,
        lower=1.00,
        upper=3.0,
        scaled=True,
)

DVCon.addThicknessConstraints2D(
        leList=leList,
        teList=teList,
        nSpan=2,
        nChord=100,
        lower=1.00,
        upper=3.0,
        scaled=True,
)

le = 0.01
leList = [[le, 0, wingtipSpacing], [le, 0, 1.0 - wingtipSpacing]]

DVCon.addLERadiusConstraints(
        leList=leList,
        nSpan=2,
        axis=[0,1,0],
        chordDir=[-1,0,0],
        lower=1.00,
        upper=3.0,
        scaled=True,
)

fileName = os.path.join(outputDir, "constraints.dat")
DVCon.writeTecplot(fileName)
# rst cons (end)

# ======================================================================
#         Functions:
# ======================================================================
# rst funcs (beg)
def cruiseFuncs(x):
    if MPI.COMM_WORLD.rank == 0:
        print(x)
    # Set design vars
    DVGeo.setDesignVars(x)
    # Evaluate functions
    funcs = {}
    DVCon.evalFunctions(funcs)

    for i in range(nFlowCases):
        if i % nGroup == ptID:
            aeroProblems[i].setDesignVars(x)
            CFDSolver(aeroProblems[i])
            CFDSolver.evalFunctions(aeroProblems[i], funcs)
            CFDSolver.checkSolutionFailure(aeroProblems[i], funcs)
    if MPI.COMM_WORLD.rank == 0:
        print(funcs)
    return funcs


def cruiseFuncsSens(x, funcs):
    funcsSens = {}
    DVCon.evalFunctionsSens(funcsSens)
    for i in range(nFlowCases):
        if i % nGroup == ptID:
            CFDSolver.evalFunctionsSens(aeroProblems[i], funcsSens)
            CFDSolver.checkAdjointFailure(aeroProblems[i], funcsSens)
    if MPI.COMM_WORLD.rank == 0:
        print(funcsSens)
    return funcsSens


def objCon(funcs, printOK):
    # Assemble the objective and any additional constraints:
    funcs["obj"] = 0.0
    for i in range(nFlowCases):
        ap = aeroProblems[i]
        funcs["obj"] += funcs[ap["cd"]] / nFlowCases
        funcs["cl_con_" + ap.name] = funcs[ap["cl"]] - mycl[i]
    if printOK:
        print("funcs in obj:", funcs)
    return funcs
# rst funcs (end)

# ======================================================================
#         Optimization Problem Set-up
# ======================================================================
# rst optprob (beg)
# Create optimization problem
optProb = Optimization("opt", MP.obj, comm=MPI.COMM_WORLD)

# Add objective
optProb.addObj("obj", scale=1e4)

# Add variables from the AeroProblem
for ap in aeroProblems:
    ap.addVariablesPyOpt(optProb)

# Add DVGeo variables
DVGeo.addVariablesPyOpt(optProb)

# Add constraints
DVCon.addConstraintsPyOpt(optProb)
for ap in aeroProblems:
    optProb.addCon("cl_con_" + ap.name, lower=0.0, upper=0.0, scale=1.0)

# Enforce first upper and lower CST coefficients to add to zero
# to maintain continuity at the leading edge
jac = np.zeros((1, nCoeff), dtype=float)
jac[0, 0] = 1.0
optProb.addCon(
    "first_cst_coeff_match",
    lower=0.0,
    upper=0.0,
    linear=True,
    wrt=["upper_shape", "lower_shape"],
    jac={"upper_shape": jac, "lower_shape": jac},
)

# The MP object needs the 'obj' and 'sens' function for each proc set,
# the optimization problem and what the objcon function is:
MP.setProcSetObjFunc("cruise", cruiseFuncs)
MP.setProcSetSensFunc("cruise", cruiseFuncsSens)
MP.setObjCon(objCon)
MP.setOptProb(optProb)
optProb.printSparsity()
optProb.getDVConIndex()
# rst optprob (end)

# rst optimizer
# Set up optimizer
if args.opt == "SLSQP":
    optOptions = {"IFILE": os.path.join(args.output, "SLSQP.out")}
elif args.opt == "SNOPT":
    optOptions = {
        "Major feasibility tolerance": 1e-4,
        "Major optimality tolerance": 1e-4,
        "Hessian full memory": None,
        "Function precision": 1e-8,
        "Print file": os.path.join(args.output, "SNOPT_print.out"),
        "Summary file": os.path.join(args.output, "SNOPT_summary.out"),
    }
optOptions.update(args.optOptions)
opt = OPT(args.opt, options=optOptions)

# Run Optimization
sol = opt(optProb, MP.sens, storeHistory=os.path.join(args.output, "opt.hst"))
if MPI.COMM_WORLD.rank == 0:
    print(sol)

# ======================================================================
#         Postprocessing
# ======================================================================
# rst postprocessing (beg)
# Save the final figure
CFDSolver.airfoilAxs[1].legend(["Original", "Optimized"], labelcolor="linecolor")
CFDSolver.airfoilFig.savefig(os.path.join(outputDir, "OptFoil.pdf"))

# Animate the optimization

for i in range(nFlowCases):
    AnimateAirfoilOpt(outputDir, "fc%d" % i).animate(
        outputFileName=os.path.join(outputDir, "OptFoil"), fps=10, dpi=300, extra_args=["-vcodec", "libx264"]
    )
# rst postprocessing (end)

パラメータ設定

パラメータ設定は以下の部分

DVCon.addVolumeConstraint(
        leList=leList,
        teList=teList,
        nSpan=2,
        nChord=100,
        lower=1.00,
        upper=3.0,
        scaled=True,
)
DVCon.addThicknessConstraints2D(
        leList=leList,
        teList=teList,
        nSpan=2,
        nChord=100,
        lower=1.00,
        upper=3.0,
        scaled=True,
)
DVCon.addLERadiusConstraints(
        leList=leList,
        nSpan=2,
        axis=[0,1,0],
        chordDir=[-1,0,0],
        lower=1.00,
        upper=3.0,
        scaled=True,
)

基本的にいじるのはlowerupper

ベースとなる翼型の形状を1.0とし,新たに作る翼型の形状の下限と上限を倍率で指定する

指定できるのは翼型の面積(Volume),翼厚(Thickness),前縁半径(LERadius)の3つ

パラメータ設定はよくわからん

CL=1.0で最適化する

single_point.pyを使い,CL=1.0,Re=500,000でNACA4412を最適化する

# ======================================================================
#         Specify parameters for optimization
# ======================================================================
# rst params (beg)
mycl = 1.0  # lift coefficient constraint
alpha = 5.0 # initial angle of attack (zero if the target cl is zero)
mach = 0.0  # Mach number
Re = 5e5  # Reynolds number
T = 288.15  # 1976 US Standard Atmosphere temperature @ sea level (K)
# rst params (end)

airfoilname = "naca4412.dat"

次のコマンドでoutputディレクトリを削除しつつsingle_point.pyを実行する

rm -r output | python single_point.py
(.py310venv) :/mnt/c/Users/mtkbirdman/Desktop/CMPLXFOIL/naca4412$ rm -r output | python single_point.py
/mnt/c/Users/mtkbirdman/Desktop/CMPLXFOIL/.py310venv/lib/python3.10/site-packages/baseclasses/problems/pyAero_problem.py:928: RuntimeWarning: divide by zero encountered in double_scalars
  self.__dict__["rho"] = self.re * self.mu / self.V
/mnt/c/Users/mtkbirdman/Desktop/CMPLXFOIL/.py310venv/lib/python3.10/site-packages/baseclasses/problems/pyAero_problem.py:937: RuntimeWarning: invalid value encountered in double_scalars
  self.q = 0.5 * self.rho * self.V**2
######## Fitting CST coefficients to coordinates in /mnt/c/Users/mtkbirdman/Desktop/CMPLXFOIL/naca4412/naca4412.dat ########
Upper surface
    L2 norm of coordinates in dat file versus fit coordinates: 0.0007718469932335163
    Fit CST coefficients: [0.19840079 0.31210138 0.21393084 0.30791771]
Lower surface
    L2 norm of coordinates in dat file versus fit coordinates: 0.0008112070833303378
    Fit CST coefficients: [-0.14295354  0.00227122 -0.0548071  -0.00507369]
pyGeo Initialization Type is: liftingSurface
Interpolating missing sections ...
Symmetrizing Knot Vectors ...
Computing TE surfaces ...
Computing Tip surfaces ...
+-------------------------------------------------------------------------------+
|                   Sparsity structure of constraint Jacobian                   |
+-------------------------------------------------------------------------------+
                                          alpha_fc (1)   upper_shape (4)   lower_shape (4)
                                        +--------------------------------------------------+
         DVCon1_volume_constraint_0 (1) |              |        X        |        X        |
                                        ----------------------------------------------------
   DVCon1_thickness_constraints_0 (200) |              |        X        |        X        |
                                        ----------------------------------------------------
      DVCon1_leradius_constraints_0 (2) |              |        X        |        X        |
                                        ----------------------------------------------------
                          cl_con_fc (1) |       X      |        X        |        X        |
                                        ----------------------------------------------------
           first_cst_coeff_match(L) (1) |              |        X        |        X        |
                                        +--------------------------------------------------+
### DESIGN VARIABLES ###
alpha_fc [1]
upper_shape [2, 5]
lower_shape [6, 9]
### CONSTRAINTS ###
DVCon1_volume_constraint_0 [1]
DVCon1_thickness_constraints_0 [2, 201]
DVCon1_leradius_constraints_0 [202, 203]
cl_con_fc [204]
first_cst_coeff_match [205]
OrderedDict([('alpha_fc', array(5.)), ('upper_shape', array([0.19840079, 0.31210138, 0.21393084, 0.30791771])), ('lower_shape', array([-0.14295354,  0.00227122, -0.0548071 , -0.00507369]))])
+----------------------------------------------------------------------+
|  Switching to Aero Problem: fc                                       |
+----------------------------------------------------------------------+
functions:
    DVCon1_volume_constraint_0: 1.01573699552304
    DVCon1_leradius_constraints_0: [0.97764811 0.97764811]
    fc_cd: 0.00944166942668497
    fc_cl: 1.0145047127668452
    fail: False
funcs in obj: {'fc_cd': 0.00944166942668497, 'fc_cl': 1.0145047127668452, 'obj': 0.00944166942668497, 'cl_con_fc': 0.014504712766845174}
function sensitivities:
    fc_cd: {'upper_shape': array([ 0.01076249,  0.00732753, -0.00246365,  0.00299908]), 'lower_shape': array([-0.00229123,  0.00026035,  0.00104259,  0.00110376]), 'alpha_fc': array([0.00053889])}
    fc_cl: {'upper_shape': array([0.25253025, 0.626825  , 0.60554456, 0.41360192]), 'lower_shape': array([0.29386779, 0.3734998 , 0.54493649, 0.74788055]), 'alpha_fc': array([0.10007496])}
    fail: False
OrderedDict([('alpha_fc', array(5.30475417)), ('upper_shape', array([ 0.17304847, -0.1       ,  0.5       ,  0.38229858])), ('lower_shape', array([-0.17304847, -0.36821269,  0.1       ,  0.1       ]))])
	:
OrderedDict([('alpha_fc', array(3.40753513)), ('upper_shape', array([0.17242851, 0.24545086, 0.25624909, 0.4376909 ])), ('lower_shape', array([-0.17242851, -0.06437285,  0.00197058,  0.1       ]))])
functions:
    DVCon1_volume_constraint_0: 1.018828032255285
    DVCon1_leradius_constraints_0: [1.00453248 1.00453248]
    fc_cd: 0.007317203636262957
    fc_cl: 0.9999993797917858
    fail: True
funcs in obj: {'fc_cd': 0.007317203636262957, 'fc_cl': 0.9999993797917858, 'obj': 0.007317203636262957, 'cl_con_fc': -6.202082142303667e-07}


Optimization Problem -- opt
================================================================================
    Objective Function: obj

    Solution:
--------------------------------------------------------------------------------
    Total Time:                  518.1149
       User Objective Time :      46.9870
       User Sensitivity Time :   469.5620
       Interface Time :            1.5383
       Opt Solver Time:            0.0276
    Calls to Objective Function :     142
    Calls to Sens Function :           48


   Objectives
      Index  Name            Value
          0  obj     7.317204E-03

   Variables (c - continuous, i - integer, d - discrete)
      Index  Name            Type      Lower Bound            Value      Upper Bound     Status
          0  alpha_fc_0         c    -1.000000E+01     3.407535E+00     1.000000E+01
          1  upper_shape_0      c    -1.000000E-01     1.724285E-01     5.000000E-01
          2  upper_shape_1      c    -1.000000E-01     2.454509E-01     5.000000E-01
          3  upper_shape_2      c    -1.000000E-01     2.562491E-01     5.000000E-01
          4  upper_shape_3      c    -1.000000E-01     4.376909E-01     5.000000E-01
          5  lower_shape_0      c    -5.000000E-01    -1.724285E-01     1.000000E-01
          6  lower_shape_1      c    -5.000000E-01    -6.437285E-02     1.000000E-01
          7  lower_shape_2      c    -5.000000E-01     1.970583E-03     1.000000E-01
          8  lower_shape_3      c    -5.000000E-01     1.000000E-01     1.000000E-01          u

   Constraints (i - inequality, e - equality)
      Index  Name                           Type          Lower           Value           Upper    Status  Lagrange Multiplier (N/A)
          0  DVCon1_volume_constraint_0        i   1.000000E+00    1.018828E+00    3.000000E+00              9.00000E+100
          1  DVCon1_thickness_constraints_0    i   1.000000E+00    1.377185E+00    3.000000E+00              9.00000E+100
          2  DVCon1_thickness_constraints_0    i   1.000000E+00    1.000000E+00    3.000000E+00         l    9.00000E+100
          3  DVCon1_thickness_constraints_0    i   1.000000E+00    1.007753E+00    3.000000E+00              9.00000E+100
		:
        202  DVCon1_leradius_constraints_0     i   1.000000E+00    1.004532E+00    3.000000E+00              9.00000E+100
        203  cl_con_fc                         e   0.000000E+00   -6.202082E-07    0.000000E+00              9.00000E+100
        204  first_cst_coeff_match             e   0.000000E+00    0.000000E+00    0.000000E+00              9.00000E+100

--------------------------------------------------------------------------------

Found dat and pkl files for 142 iterations

できた翼型がこれ

1.0000000		0.00126
0.9928006		0.004362755437
0.9798895		0.009734052047
0.9635184		0.016200612922
0.9445469		0.023234834523
0.9234998		0.030492140422
0.9007495		0.037728539638
0.8765830		0.044767730463
0.8512328		0.051483433617
0.8248944		0.05778756071
0.7977362		0.063621676264
0.7699063		0.068950151107
0.7415373		0.073754604784
0.7127487		0.078029529869
0.6836500		0.081778464579
0.6543419		0.085011109734
0.6249182		0.087740925379
0.5954661		0.089983433741
0.5660677		0.091754864337
0.5368005		0.093071316971
0.5077378		0.093948262972
0.4789491		0.094400323555
0.4505011		0.094441276135
0.4224571		0.094084251652
0.3948782		0.093342029775
0.3678229		0.09222740093
0.3413476		0.090753550993
0.3155065		0.088934407252
0.2903524		0.086784993668
0.2659361		0.084321618797
0.2423068		0.081562027961
0.2195125		0.078525477924
0.1975998		0.075232638545
0.1766138		0.071705376803
0.1565988		0.067966578035
0.1375977		0.064039671772
0.1196525		0.059948271117
0.1028043		0.055715688009
0.0870931		0.051364360218
0.0725582		0.046915428646
0.0592379		0.04238815356
0.0471700		0.037799627837
0.0363913		0.033164232658
0.0269381		0.028493545349
0.0188458		0.023795992046
0.0121494		0.019077060088
0.0068832		0.014339263142
0.0030809		0.009582436222
0.0007756		0.004804042134
0.0000000		-0.0
0.0007756		-0.004792322516
0.0030809		-0.009489988345
0.0068832		-0.014032385463
0.0121494		-0.018363373255
0.0188458		-0.022431689919
0.0269381		-0.026191622558
0.0363913		-0.029603434007
0.0471700		-0.032633796885
0.0592379		-0.035255828258
0.0725582		-0.037449315474
0.0870931		-0.039200604034
0.1028043		-0.040502570849
0.1196525		-0.041354353259
0.1375977		-0.041761098266
0.1565988		-0.041733595549
0.1766138		-0.041287885366
0.1975998		-0.040444799679
0.2195125		-0.039229501814
0.2423068		-0.037670960694
0.2659361		-0.035801476279
0.2903524		-0.033656172902
0.3155065		-0.031272477679
0.3413476		-0.028689666096
0.3678229		-0.02594843786
0.3948782		-0.023090409287
0.4224571		-0.020157769435
0.4505011		-0.017192864918
0.4789491		-0.014237861094
0.5077378		-0.011334317063
0.5368005		-0.00852292468
0.5660677		-0.005843050872
0.5954661		-0.003332415119
0.6249182		-0.00102666035
0.6543419		0.001041131964
0.6836500		0.002840987261
0.7127487		0.004346616477
0.7415373		0.005536187004
0.7699063		0.00639313423
0.7977362		0.006907162828
0.8248944		0.007075369626
0.8512328		0.006903591675
0.876583		0.006407949942
0.9007495		0.00561670371
0.9234998		0.00457251103
0.9445469		0.003335411539
0.9635184		0.001987167731
0.9798895		0.000639081629
0.9928006		-0.000548737883
1.0000000		-0.00126

XFLR5で解析してみるとこうなる

CL=1.0で猛烈にとがっている

CL=0.6, 1.0で最適化する

multi_point.pyを使い,CL=[0.6, 1.0],Re=500,000でNACA4412を最適化する

# ======================================================================
#         Specify parameters for optimization
# ======================================================================
# rst params (beg)
mycl = [0.6,1.0]  # lift coefficient constraint
alpha = [2.0,5.0] # initial angle of attack (zero if the target cl is zero)
mach = [0.0,0.0]  # Mach number
Re = [5e5,5e5]  # Reynolds number
T = [288.15,288.15]  # 1976 US Standard Atmosphere temperature @ sea level (K)
nFlowCases = len(Re)
# rst params (end)

airfoilname = "naca4412.dat"

次のコマンドでoutputディレクトリを削除しつつmulti_point.pyを実行する

rm -r output | python multi_point.py

(.py310venv) :/mnt/c/Users/mtkbirdman/Desktop/CMPLXFOIL/naca4412$ rm -r output | python multi_point.py
/mnt/c/Users/mtkbirdman/Desktop/CMPLXFOIL/.py310venv/lib/python3.10/site-packages/baseclasses/problems/pyAero_problem.py:928: RuntimeWarning: divide by zero encountered in double_scalars
  self.__dict__["rho"] = self.re * self.mu / self.V
/mnt/c/Users/mtkbirdman/Desktop/CMPLXFOIL/.py310venv/lib/python3.10/site-packages/baseclasses/problems/pyAero_problem.py:937: RuntimeWarning: invalid value encountered in double_scalars
  self.q = 0.5 * self.rho * self.V**2
######## Fitting CST coefficients to coordinates in /mnt/c/Users/mtkbirdman/Desktop/CMPLXFOIL/naca4412/naca4412.dat ########
Upper surface
    L2 norm of coordinates in dat file versus fit coordinates: 0.0007718469932335163
    Fit CST coefficients: [0.19840079 0.31210138 0.21393084 0.30791771]
Lower surface
    L2 norm of coordinates in dat file versus fit coordinates: 0.0008112070833303378
    Fit CST coefficients: [-0.14295354  0.00227122 -0.0548071  -0.00507369]
pyGeo Initialization Type is: liftingSurface
Interpolating missing sections ...
Symmetrizing Knot Vectors ...
Computing TE surfaces ...
Computing Tip surfaces ...
+-------------------------------------------------------------------------------+
|                   Sparsity structure of constraint Jacobian                   |
+-------------------------------------------------------------------------------+
                                          alpha_fc0 (1)   alpha_fc1 (1)   upper_shape (4)   lower_shape (4)
                                        +-------------------------------------------------------------------+
         DVCon1_volume_constraint_0 (1) |               |               |        X        |        X        |
                                        ---------------------------------------------------------------------
   DVCon1_thickness_constraints_0 (200) |               |               |        X        |        X        |
                                        ---------------------------------------------------------------------
      DVCon1_leradius_constraints_0 (2) |               |               |        X        |        X        |
                                        ---------------------------------------------------------------------
                         cl_con_fc0 (1) |       X       |       X       |        X        |        X        |
                                        ---------------------------------------------------------------------
                         cl_con_fc1 (1) |       X       |       X       |        X        |        X        |
                                        ---------------------------------------------------------------------
           first_cst_coeff_match(L) (1) |               |               |        X        |        X        |
                                        +-------------------------------------------------------------------+
### DESIGN VARIABLES ###
alpha_fc0 [1]
alpha_fc1 [2]
upper_shape [3, 6]
lower_shape [7, 10]
### CONSTRAINTS ###
DVCon1_volume_constraint_0 [1]
DVCon1_thickness_constraints_0 [2, 201]
DVCon1_leradius_constraints_0 [202, 203]
cl_con_fc0 [204]
cl_con_fc1 [205]
first_cst_coeff_match [206]
OrderedDict([('alpha_fc0', array(2.)), ('alpha_fc1', array(5.)), ('upper_shape', array([0.19840079, 0.31210138, 0.21393084, 0.30791771])), ('lower_shape', array([-0.14295354,  0.00227122, -0.0548071 , -0.00507369]))])
+----------------------------------------------------------------------+
|  Switching to Aero Problem: fc0                                      |
+----------------------------------------------------------------------+
+----------------------------------------------------------------------+
|  Switching to Aero Problem: fc1                                      |
+----------------------------------------------------------------------+
{'DVCon1_volume_constraint_0': 1.01573699552304, 'DVCon1_thickness_constraints_0': array([1.36320157, 0.99014066, 0.99812991, 1.00273913, 1.00562219,
       1.00752997, 1.0089155 , 1.0098028 , 1.01033061, 1.01059078,
       1.01078522, 1.01078866, 1.01068377, 1.01052361, 1.01027269,
       1.01000679, 1.00970039, 1.00940716, 1.00907138, 1.00876555,
		:
       1.05979155, 1.06195505, 1.06403079, 1.06612797, 1.06830362,
       1.07035994, 1.07225836, 1.07391039, 1.07525838, 1.07603075,
       1.07595793, 1.07413975, 1.06874689, 1.05492332, 1.00120263]), 'DVCon1_leradius_constraints_0': array([0.97764811, 0.97764811]), 'fc0_cd': 0.007742382075777451, 'fc0_cl': 0.7040047225796933, 'fail': False, 'fc1_cd': 0.00944166942668497, 'fc1_cl': 1.0145047127668452}
funcs in obj: {'fc0_cd': 0.007742382075777451, 'fc0_cl': 0.7040047225796933, 'fc1_cd': 0.00944166942668497, 'fc1_cl': 1.0145047127668452, 'obj': 0.00859202575123121, 'cl_con_fc0': 0.10400472257969329, 'cl_con_fc1': 0.014504712766845174}
+----------------------------------------------------------------------+
|  Switching to Aero Problem: fc0                                      |
+----------------------------------------------------------------------+
+----------------------------------------------------------------------+
|  Switching to Aero Problem: fc1                                      |
+----------------------------------------------------------------------+
{'DVCon1_volume_constraint_0': {'upper_shape': array([0.89743927, 1.0118917 , 0.84325189, 0.4917935 ]), 'lower_shape': array([-0.89743927, -1.0118917 , -0.84325189, -0.4917935 ])}, 'DVCon1_thickness_constraints_0': {'upper_shape': array([[3.99212685e+00, 1.19775783e-03, 1.19787762e-07, 3.99332473e-12],
       [2.81864126e+00, 8.71305246e-02, 8.97799993e-04, 3.08366799e-06],
       [2.76222033e+00, 1.71687384e-01, 3.55710919e-03, 2.45660173e-05],
       [2.69703695e+00, 2.53655613e-01, 7.95208116e-03, 8.30990149e-05],
		:

Optimization Problem -- opt
================================================================================
    Objective Function: obj

    Solution:
--------------------------------------------------------------------------------
    Total Time:                  430.7806
       User Objective Time :      29.3784
       User Sensitivity Time :   400.6035
       Interface Time :            0.7834
       Opt Solver Time:            0.0154
    Calls to Objective Function :      60
    Calls to Sens Function :           36


   Objectives
      Index  Name            Value
          0  obj     7.617882E-03

   Variables (c - continuous, i - integer, d - discrete)
      Index  Name            Type      Lower Bound            Value      Upper Bound     Status
          0  alpha_fc0_0        c     0.000000E+00     9.283732E-01     1.000000E+01
          1  alpha_fc1_0        c     0.000000E+00     4.693599E+00     1.000000E+01
          2  upper_shape_0      c    -1.000000E-01     1.728902E-01     5.000000E-01
          3  upper_shape_1      c    -1.000000E-01     2.611330E-01     5.000000E-01
          4  upper_shape_2      c    -1.000000E-01     2.810529E-01     5.000000E-01
          5  upper_shape_3      c    -1.000000E-01     3.359472E-01     5.000000E-01
          6  lower_shape_0      c    -5.000000E-01    -1.728902E-01     1.000000E-01
          7  lower_shape_1      c    -5.000000E-01    -1.793646E-02     1.000000E-01
          8  lower_shape_2      c    -5.000000E-01    -5.914019E-02     1.000000E-01
          9  lower_shape_3      c    -5.000000E-01     6.752422E-03     1.000000E-01

   Constraints (i - inequality, e - equality)
      Index  Name                           Type          Lower           Value           Upper    Status  Lagrange Multiplier (N/A)
          0  DVCon1_volume_constraint_0        i   1.000000E+00    1.056806E+00    3.000000E+00              9.00000E+100
          1  DVCon1_thickness_constraints_0    i   1.000000E+00    1.380834E+00    3.000000E+00              9.00000E+100
          2  DVCon1_thickness_constraints_0    i   1.000000E+00    1.000000E+00    3.000000E+00         l    9.00000E+100
          3  DVCon1_thickness_constraints_0    i   1.000000E+00    1.005329E+00    3.000000E+00              9.00000E+100
		:
        203  cl_con_fc0                        e   0.000000E+00    1.322324E-08    0.000000E+00              9.00000E+100
        204  cl_con_fc1                        e   0.000000E+00   -1.768477E-09    0.000000E+00              9.00000E+100
        205  first_cst_coeff_match             e   0.000000E+00    0.000000E+00    0.000000E+00              9.00000E+100

--------------------------------------------------------------------------------

Found dat and pkl files for 60 iterations
Found dat and pkl files for 60 iterations

できた翼型がこれ

1.0000000		0.00126
0.9928006		0.003652358783
0.9798895		0.007857346024
0.9635184		0.013033968669
0.9445469		0.01881982222
0.9234998		0.024976621602
0.9007495		0.031328206848
0.876583		0.037737037963
0.8512328		0.044092675997
0.8248944		0.050304825869
0.7977362		0.056298942422
0.7699063		0.062013028785
0.7415373		0.067395272715
0.7127487		0.072402474371
0.6836500		0.076998594775
0.6543419		0.081153883606
0.6249182		0.084844081953
0.5954661		0.088050014429
0.5660677		0.090757145444
0.5368005		0.092955376459
0.5077378		0.094638908646
0.4789491		0.095806139989
0.4505011		0.096459572556
0.4224571		0.096605781854
0.3948782		0.096255309466
0.3678229		0.095422557663
0.3413476		0.094125620615
0.3155065		0.092386045312
0.2903524		0.090228584662
0.2659361		0.087680797746
0.2423068		0.084772635874
0.2195125		0.081535999209
0.1975998		0.078004171146
0.1766138		0.074211194446
0.1565988		0.070191362195
0.1375977		0.065978466416
0.1196525		0.061605241531
0.1028043		0.05710274205
0.0870931		0.052499687432
0.0725582		0.047822019446
0.0592379		0.043092332213
0.0471700		0.038329648186
0.0363913		0.033548949424
0.0269381		0.028761195986
0.0188458		0.023973089421
0.0121494		0.019187414004
0.0068832		0.014403122667
0.0030809		0.009615756874
0.0007756		0.004817874515
0.0000000		-0.0
0.0007756		-0.004802140772
0.0030809		-0.009491791695
0.0068832		-0.013992433801
0.0121494		-0.018234920724
0.0188458		-0.022158693246
0.0269381		-0.025713038444
0.0363913		-0.028857876165
0.0471700		-0.031564305179
0.0592379		-0.033814541243
0.0725582		-0.035601801259
0.0870931		-0.036929675372
0.1028043		-0.037811378743
0.1196525		-0.038268643665
0.1375977		-0.038330498378
0.1565988		-0.038031872854
0.1766138		-0.037412130731
0.1975998		-0.036513547851
0.2195125		-0.035379833601
0.2423068		-0.034054667781
0.2659361		-0.032580407255
0.2903524		-0.030996908021
0.3155065		-0.029340520924
0.3413476		-0.027643342467
0.3678229		-0.02593270239
0.3948782		-0.02423083289
0.4224571		-0.022554879997
0.4505011		-0.020917099517
0.4789491		-0.019325316228
0.5077378		-0.017783532467
0.5368005		-0.016292791071
0.5660677		-0.014852060437
0.5954661		-0.013459274335
0.6249182		-0.012112348392
0.6543419		-0.010810144861
0.6836500		-0.009553323614
0.7127487		-0.008345048153
0.7415373		-0.007191385224
0.7699063		-0.006101455284
0.7977362		-0.005087180928
0.8248944		-0.00416268814
0.8512328		-0.003343226461
0.876583		-0.002643702381
0.9007495		-0.002076794111
0.9234998		-0.001650686382
0.9445469		-0.001366531547
0.9635184		-0.001215691159
0.9798895		-0.00117684023
0.9928006		-0.001212580992
1.0000000		-0.00126

XFLR5で解析してみるとこうなる

CL=0.6とCL=1.0の両方の点でNACA4412を上回る翼型ができた

おわりに

CMPLXFOILを使って2次元翼型の最適化を行った

最適化の設計点とか形状の制約パラメータをうまくいじればいい感じに使えるんじゃなかろうか

↓関連記事

コメント

  1. SnipeLover より:

    連日のコメント失礼します.
    もしも間違った指摘ならごめんなさい,

    冒頭のディレクトリ作成についてですが
    mkdir -p /mnt/c/Users/mtkbirdman/Desktop/CMPLXFOIL/naca4412birdman/Desktop/CMPLXFOIL/naca4412

    に対して、
    パイソンスクリプト実行前に移動したディレクトリが

    (.py310venv) :~$ cd /mnt/c/Users/mtkbirdman/Desktop/CMPLXFOIL/naca4412

    であるため

    naca4412birdman/Desktop/CMPLXFOIL の部分が誤植かもしれません.

    それ以外はブログ通りで,multi_point.pyまで再現出来ました.

  2. @mtk_birdman @mtk_birdman より:

    コメントありがとうございます!
    ご指摘の通り,naca4412birdman/Desktop/CMPLXFOILがコピペミスでしたので,正しいPATHに修正しました.

    翼型の最適化をお楽しみください!