CMPLXFOILを使って2次元翼型の最適化を行う
はじめに
CMPLXFOILを使って2次元翼型の最適化を行う
↓公式ドキュメント
Optimization Tutorial — CMPLXFOIL documentation
mdolab-cmplxfoil.readthedocs-hosted.com
↓インストール
↓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,
)基本的にいじるのはlowerとupper
ベースとなる翼型の形状を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次元翼型の最適化を行った
最適化の設計点とか形状の制約パラメータをうまくいじればいい感じに使えるんじゃなかろうか
↓関連記事
コメント
連日のコメント失礼します.
もしも間違った指摘ならごめんなさい,
冒頭のディレクトリ作成についてですが
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まで再現出来ました.
コメントありがとうございます！
ご指摘の通り，naca4412birdman/Desktop/CMPLXFOILがコピペミスでしたので，正しいPATHに修正しました．
翼型の最適化をお楽しみください！