From 17fac0d402ab6259d4535da43ee8e7be58f47bbe Mon Sep 17 00:00:00 2001 From: nmadison11 Date: Wed, 19 Nov 2025 13:18:50 -0600 Subject: [PATCH 1/2] Added RAE2822_Airfoil_HISA to validations --- RAE2822_Airfoil_HISA/0_orig/T | 49 ++++ RAE2822_Airfoil_HISA/0_orig/U | 51 ++++ RAE2822_Airfoil_HISA/0_orig/alphat | 45 +++ .../0_orig/include/freestreamConditions | 14 + RAE2822_Airfoil_HISA/0_orig/nuTilda | 46 +++ RAE2822_Airfoil_HISA/0_orig/nut | 46 +++ RAE2822_Airfoil_HISA/0_orig/p | 51 ++++ RAE2822_Airfoil_HISA/Allclean.sh | 24 ++ RAE2822_Airfoil_HISA/FFD/genFFD.py | 99 +++++++ RAE2822_Airfoil_HISA/FFD/wingFFD.xyz | 5 + .../constant/thermophysicalProperties | 50 ++++ .../constant/turbulenceProperties | 31 ++ RAE2822_Airfoil_HISA/genAirFoilMesh.py | 269 ++++++++++++++++++ RAE2822_Airfoil_HISA/paraview.foam | 0 RAE2822_Airfoil_HISA/preProcessing.sh | 20 ++ .../profiles/RAE2822PS.profile | 63 ++++ .../profiles/RAE2822SS.profile | 63 ++++ RAE2822_Airfoil_HISA/runScript.py | 268 +++++++++++++++++ RAE2822_Airfoil_HISA/system/controlDict | 39 +++ RAE2822_Airfoil_HISA/system/createPatchDict | 82 ++++++ RAE2822_Airfoil_HISA/system/decomposeParDict | 32 +++ RAE2822_Airfoil_HISA/system/fvSchemes | 69 +++++ RAE2822_Airfoil_HISA/system/fvSolution | 74 +++++ 23 files changed, 1490 insertions(+) create mode 100644 RAE2822_Airfoil_HISA/0_orig/T create mode 100644 RAE2822_Airfoil_HISA/0_orig/U create mode 100644 RAE2822_Airfoil_HISA/0_orig/alphat create mode 100644 RAE2822_Airfoil_HISA/0_orig/include/freestreamConditions create mode 100644 RAE2822_Airfoil_HISA/0_orig/nuTilda create mode 100644 RAE2822_Airfoil_HISA/0_orig/nut create mode 100644 RAE2822_Airfoil_HISA/0_orig/p create mode 100644 RAE2822_Airfoil_HISA/Allclean.sh create mode 100644 RAE2822_Airfoil_HISA/FFD/genFFD.py create mode 100644 RAE2822_Airfoil_HISA/FFD/wingFFD.xyz create mode 100644 RAE2822_Airfoil_HISA/constant/thermophysicalProperties create mode 100644 RAE2822_Airfoil_HISA/constant/turbulenceProperties create mode 100644 RAE2822_Airfoil_HISA/genAirFoilMesh.py create mode 100644 RAE2822_Airfoil_HISA/paraview.foam create mode 100644 RAE2822_Airfoil_HISA/preProcessing.sh create mode 100644 RAE2822_Airfoil_HISA/profiles/RAE2822PS.profile create mode 100644 RAE2822_Airfoil_HISA/profiles/RAE2822SS.profile create mode 100644 RAE2822_Airfoil_HISA/runScript.py create mode 100644 RAE2822_Airfoil_HISA/system/controlDict create mode 100644 RAE2822_Airfoil_HISA/system/createPatchDict create mode 100644 RAE2822_Airfoil_HISA/system/decomposeParDict create mode 100644 RAE2822_Airfoil_HISA/system/fvSchemes create mode 100644 RAE2822_Airfoil_HISA/system/fvSolution diff --git a/RAE2822_Airfoil_HISA/0_orig/T b/RAE2822_Airfoil_HISA/0_orig/T new file mode 100644 index 0000000..ae3a909 --- /dev/null +++ b/RAE2822_Airfoil_HISA/0_orig/T @@ -0,0 +1,49 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: plus | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + location "0"; + object T; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +#include "include/freestreamConditions" +dimensions [0 0 0 1 0 0 0]; + +internalField uniform $T; + +boundaryField +{ + "wing.*" + { + type characteristicWallTemperature; + value $internalField; + } + symmetry1 + { + type symmetry; + } + symmetry2 + { + type symmetry; + } + inout + { + type characteristicFarfieldTemperature; + #include "include/freestreamConditions" + value $internalField; + //type inletOutlet; + //inletValue $internalField; + //value $internalField; + } +} + + +// ************************************************************************* // diff --git a/RAE2822_Airfoil_HISA/0_orig/U b/RAE2822_Airfoil_HISA/0_orig/U new file mode 100644 index 0000000..5122bab --- /dev/null +++ b/RAE2822_Airfoil_HISA/0_orig/U @@ -0,0 +1,51 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: plus | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volVectorField; + location "0"; + object U; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#include "include/freestreamConditions" + +dimensions [0 1 -1 0 0 0 0]; + +internalField uniform $U; + +boundaryField +{ + "wing.*" + { + type boundaryCorrectedFixedValue; + value uniform (0 0 0); + } + symmetry1 + { + type symmetry; + } + symmetry2 + { + type symmetry; + } + inout + { + type characteristicFarfieldVelocity; + #include "include/freestreamConditions" + value uniform (0 0 0); + //type inletOutlet; + //inletValue $internalField; + //value $internalField; + } +} + + +// ************************************************************************* // diff --git a/RAE2822_Airfoil_HISA/0_orig/alphat b/RAE2822_Airfoil_HISA/0_orig/alphat new file mode 100644 index 0000000..fdddca5 --- /dev/null +++ b/RAE2822_Airfoil_HISA/0_orig/alphat @@ -0,0 +1,45 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: plus | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + location "0"; + object alphat; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [1 -1 -1 0 0 0 0]; + +internalField uniform 0; + +boundaryField +{ + "wing.*" + { + type boundaryCorrectedFixedValue; + value uniform 0; + } + symmetry1 + { + type symmetry; + } + symmetry2 + { + type symmetry; + } + inout + { + type calculated; + value $internalField; + } +} + + +// ************************************************************************* // diff --git a/RAE2822_Airfoil_HISA/0_orig/include/freestreamConditions b/RAE2822_Airfoil_HISA/0_orig/include/freestreamConditions new file mode 100644 index 0000000..98ce389 --- /dev/null +++ b/RAE2822_Airfoil_HISA/0_orig/include/freestreamConditions @@ -0,0 +1,14 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| | +| HiSA: High Speed Aerodynamic solver | +| Copyright (C) 2014-2017 Johan Heyns - CSIR, South Africa | +| Copyright (C) 2014-2017 Oliver Oxtoby - CSIR, South Africa | +| | +\*---------------------------------------------------------------------------*/ + +U (233.603 0 0); +p 108988.0; +T 255.56; +#inputMode merge + +// ************************************************************************* // diff --git a/RAE2822_Airfoil_HISA/0_orig/nuTilda b/RAE2822_Airfoil_HISA/0_orig/nuTilda new file mode 100644 index 0000000..b6e139b --- /dev/null +++ b/RAE2822_Airfoil_HISA/0_orig/nuTilda @@ -0,0 +1,46 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: plus | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + location "0"; + object nuTilda; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 2 -1 0 0 0 0]; + +internalField uniform 4.5e-05; + +boundaryField +{ + "wing.*" + { + type boundaryCorrectedFixedValue; + value uniform 0.0; + } + symmetry1 + { + type symmetry; + } + symmetry2 + { + type symmetry; + } + inout + { + type inletOutlet; + inletValue $internalField; + value $internalField; + } +} + + +// ************************************************************************* // diff --git a/RAE2822_Airfoil_HISA/0_orig/nut b/RAE2822_Airfoil_HISA/0_orig/nut new file mode 100644 index 0000000..9590a53 --- /dev/null +++ b/RAE2822_Airfoil_HISA/0_orig/nut @@ -0,0 +1,46 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: plus | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + location "0"; + object nut; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 2 -1 0 0 0 0]; + +internalField uniform 0; + +boundaryField +{ + "wing.*" + { + type boundaryCorrectedNutLowReWallFunction; + value uniform 0; + } + symmetry1 + { + type symmetry; + } + symmetry2 + { + type symmetry; + } + inout + { + type calculated; + value $internalField; + } + +} + + +// ************************************************************************* // diff --git a/RAE2822_Airfoil_HISA/0_orig/p b/RAE2822_Airfoil_HISA/0_orig/p new file mode 100644 index 0000000..0a2fb93 --- /dev/null +++ b/RAE2822_Airfoil_HISA/0_orig/p @@ -0,0 +1,51 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: plus | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + location "0"; + object p; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +#include "include/freestreamConditions" + +dimensions [1 -1 -2 0 0 0 0]; + +internalField uniform $p; + +boundaryField +{ + "wing.*" + { + type characteristicWallPressure; + value $internalField; + } + symmetry1 + { + type symmetry; + } + symmetry2 + { + type symmetry; + } + inout + { + type characteristicFarfieldPressure; + #include "include/freestreamConditions" + value $internalField; + + //type outletInlet; + //outletValue $internalField; + //value $internalField; + } +} + + +// ************************************************************************* // diff --git a/RAE2822_Airfoil_HISA/Allclean.sh b/RAE2822_Airfoil_HISA/Allclean.sh new file mode 100644 index 0000000..8bfb3cd --- /dev/null +++ b/RAE2822_Airfoil_HISA/Allclean.sh @@ -0,0 +1,24 @@ +#!/bin/bash + +while true +do + read -p "Delete everything and resume to the default setup (y/n)?" yn + case $yn in + [Yy]* ) + # clean everyting + echo "Cleaning..." + rm -rf 0 + rm -rf postProcessing + rm -rf constant/extendedFeatureEdgeMesh + rm -rf constant/triSurface + rm -rf constant/polyMesh/ + rm -rf *.bin *.info *.dat *.xyz *.stl *.txt + rm -rf processor* 0.0000* + rm -rf {1..9}* + exit + ;; + [Nn]* ) exit;; + * ) echo "Please answer yes or no.";; + esac +done + diff --git a/RAE2822_Airfoil_HISA/FFD/genFFD.py b/RAE2822_Airfoil_HISA/FFD/genFFD.py new file mode 100644 index 0000000..b35e34d --- /dev/null +++ b/RAE2822_Airfoil_HISA/FFD/genFFD.py @@ -0,0 +1,99 @@ +import numpy as np +import sys + +def writeFFDFile(fileName,nBlocks,nx,ny,nz,points): + ''' + Take in a set of points and write the plot 3dFile + ''' + + f = open(fileName,'w') + + f.write('%d\n'%nBlocks) + for i in range(nBlocks): + f.write('%d %d %d '%(nx[i],ny[i],nz[i])) + # end + f.write('\n') + for block in range(nBlocks): + for k in range(nz[block]): + for j in range(ny[block]): + for i in range(nx[block]): + f.write('%f '%points[block][i,j,k,0]) + # end + # end + # end + f.write('\n') + + for k in range(nz[block]): + for j in range(ny[block]): + for i in range(nx[block]): + f.write('%f '%points[block][i,j,k,1]) + # end + # end + # end + f.write('\n') + + for k in range(nz[block]): + for j in range(ny[block]): + for i in range(nx[block]): + f.write('%f '%points[block][i,j,k,2]) + # end + # end + # end + # end + f.close() + return + +def returnBlockPoints(corners,nx,ny,nz): + ''' + corners needs to be 8 x 3 + ''' + points = np.zeros([nx,ny,nz,3]) + + # points 1 - 4 are the iMin face + # points 5 - 8 are the iMax face + + for idim in range(3): + edge1 = np.linspace(corners[0][idim],corners[4][idim],nx) + edge2 = np.linspace(corners[1][idim],corners[5][idim],nx) + edge3 = np.linspace(corners[2][idim],corners[6][idim],nx) + edge4 = np.linspace(corners[3][idim],corners[7][idim],nx) + + for i in range(nx): + edge5 = np.linspace(edge1[i],edge3[i],ny) + edge6 = np.linspace(edge2[i],edge4[i],ny) + for j in range(ny): + edge7 = np.linspace(edge5[j],edge6[j],nz) + points[i,j,:,idim] = edge7 + # end + # end + # end + + return points + +################ FFD ############## +nBlocks = 1 + +nx = [8] +ny = [2] +nz = [2] + +corners = np.zeros([nBlocks,8,3]) + +corners[0,0,:] = [-0.010,-.0700,0.0] +corners[0,1,:] = [-0.010,-0.0700,0.01] +corners[0,2,:] = [-0.010,0.0700,0.0] +corners[0,3,:] = [-0.010,0.0700,0.01] +corners[0,4,:] = [ 1.01,-.07,0.0] +corners[0,5,:] = [ 1.01,-.07,0.010] +corners[0,6,:] = [ 1.01,0.07,0.0] +corners[0,7,:] = [ 1.01,0.07,0.01] + + +points = [] +for block in range(nBlocks): + points.append(returnBlockPoints(corners[block],nx[block],ny[block],nz[block])) + +#print points +fileName = 'wingFFD.xyz' +writeFFDFile(fileName,nBlocks,nx,ny,nz,points) + diff --git a/RAE2822_Airfoil_HISA/FFD/wingFFD.xyz b/RAE2822_Airfoil_HISA/FFD/wingFFD.xyz new file mode 100644 index 0000000..26f0558 --- /dev/null +++ b/RAE2822_Airfoil_HISA/FFD/wingFFD.xyz @@ -0,0 +1,5 @@ +1 +8 2 2 +-0.010000 0.135714 0.281429 0.427143 0.572857 0.718571 0.864286 1.010000 -0.010000 0.135714 0.281429 0.427143 0.572857 0.718571 0.864286 1.010000 -0.010000 0.135714 0.281429 0.427143 0.572857 0.718571 0.864286 1.010000 -0.010000 0.135714 0.281429 0.427143 0.572857 0.718571 0.864286 1.010000 +-0.070000 -0.070000 -0.070000 -0.070000 -0.070000 -0.070000 -0.070000 -0.070000 0.070000 0.070000 0.070000 0.070000 0.070000 0.070000 0.070000 0.070000 -0.070000 -0.070000 -0.070000 -0.070000 -0.070000 -0.070000 -0.070000 -0.070000 0.070000 0.070000 0.070000 0.070000 0.070000 0.070000 0.070000 0.070000 +0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000 \ No newline at end of file diff --git a/RAE2822_Airfoil_HISA/constant/thermophysicalProperties b/RAE2822_Airfoil_HISA/constant/thermophysicalProperties new file mode 100644 index 0000000..c4a3e4c --- /dev/null +++ b/RAE2822_Airfoil_HISA/constant/thermophysicalProperties @@ -0,0 +1,50 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| | +| HiSA: High Speed Aerodynamic solver | +| Copyright (C) 2014-2017 Johan Heyns - CSIR, South Africa | +| Copyright (C) 2014-2017 Oliver Oxtoby - CSIR, South Africa | +| | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "constant"; + object thermophysicalProperties; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +thermoType +{ + type hePsiThermo; + mixture pureMixture; + transport const; + thermo hConst; + equationOfState perfectGas; + specie specie; + energy sensibleInternalEnergy; +} + + +mixture +{ + specie + { + nMoles 1; + molWeight 28.966; + } + thermodynamics + { + Cp 1005; + Hf 0; + } + transport + { + mu 1.8e-5; + Pr 0.700000; + TRef 300.000000; + } +} + +// ************************************************************************* // diff --git a/RAE2822_Airfoil_HISA/constant/turbulenceProperties b/RAE2822_Airfoil_HISA/constant/turbulenceProperties new file mode 100644 index 0000000..9aab251 --- /dev/null +++ b/RAE2822_Airfoil_HISA/constant/turbulenceProperties @@ -0,0 +1,31 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| | +| HiSA: High Speed Aerodynamic solver | +| Copyright (C) 2014-2017 Johan Heyns - CSIR, South Africa | +| Copyright (C) 2014-2017 Oliver Oxtoby - CSIR, South Africa | +| | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "constant"; + object turbulenceProperties; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +simulationType RAS; + +RAS +{ + RASModel SpalartAllmaras; + + turbulence on; + + printCoeffs on; +Prt 0.85; +} + + +// ************************************************************************* // diff --git a/RAE2822_Airfoil_HISA/genAirFoilMesh.py b/RAE2822_Airfoil_HISA/genAirFoilMesh.py new file mode 100644 index 0000000..4b9bae1 --- /dev/null +++ b/RAE2822_Airfoil_HISA/genAirFoilMesh.py @@ -0,0 +1,269 @@ +#!/usr/bin/env python +""" +This script reads a coarse airfoil profile, refine the profile using spline function, +and outputs it as surfaceMesh.xyz. Then it generate a 3D volume mesh with nSpan layers +in the z direction using pyHyp, available at https://github.com/mdolab/pyhyp +Note: the airfoil data should be seperated into PS and SS surfaces, they should start from the +LE and ends at TE. We use blunt TE so truncate the PS and SS data at about 99.8% of the chord. +""" + +from pyhyp import pyHyp +import numpy +from pyspline import * + +########## user input ################ +# 2D +prefix = "./profiles/" +airfoilProfilePS = prefix + "RAE2822PS.profile" +airfoilProfileSS = prefix + "RAE2822SS.profile" +ZSpan = 0.01 # width in the z direction +nSpan = 2 # how many points in z +# PS parameters +dX1PS = 0.0025 # first dx from the LE +Alpha1PS = 1.2 # clustering from the LE +dX2PS = 5e-4 # first dx from the TE +Alpha2PS = 1.2 # clustering from the TE +dXMaxPS = 0.015 # max dx for PS +# SS parameters +dX1SS = 0.0025 # first dx from the LE +Alpha1SS = 1.2 # clustering from the LE +dX2SS = 5e-4 # first dx from the TE +Alpha2SS = 1.2 # clustering from the TE +dXMaxSS = 0.015 # max dx for SS +# TE parameters +NpTE = 5 # number of points for blunt TE +# 3D +NpExtrude = 65 # how many points to extrude for the 3D volume mesh in y +yWall = 1e-5 # first layer mesh length +marchDist = 30.0 # march distance for extruding +########## user input ################ + +# read profiles +fPS = open(airfoilProfilePS, "r") +linesPS = fPS.readlines() +fPS.close() +xPS = [] +yPS = [] +zPS = [] +for line in linesPS: + cols = line.split() + xPS.append(float(cols[0])) + yPS.append(float(cols[1])) + +for i in range(len(xPS)): + zPS.append(0.0) + +fSS = open(airfoilProfileSS, "r") +linesSS = fSS.readlines() +fSS.close() +xSS = [] +ySS = [] +zSS = [] +for line in linesSS: + cols = line.split() + xSS.append(float(cols[0])) + ySS.append(float(cols[1])) +for i in range(len(xSS)): + zSS.append(0.0) + +# ------------PS +# first compute how many stretching points do we need for both ends +tmp = dX1PS +for i in range(1000): + if tmp > dXMaxPS: + nStretch1PS = i + break + else: + tmp = tmp * Alpha1PS +# print (nStretch1PS) +tmp = dX2PS +for i in range(1000): + if tmp > dXMaxPS: + nStretch2PS = i + break + else: + tmp = tmp * Alpha2PS +# print (nStretch2PS) + +# now compute how much length does these two stretching end and the constant portion have +xLPSConst = xPS[-1] +xLPS1 = 0 +xLPS2 = 0 +for i in range(nStretch1PS): + xLPS1 += dX1PS * (Alpha1PS ** i) + xLPSConst -= dX1PS * (Alpha1PS ** i) +for i in range(nStretch2PS): + xLPS2 += dX2PS * (Alpha2PS ** i) + xLPSConst -= dX2PS * (Alpha2PS ** i) +# print xLPS1,xLPS2,xLPSConst +# update dXMax, we just want to make sure these three portion add up +nXConstPS = int(xLPSConst / dXMaxPS) +dXMaxPS = xLPSConst / nXConstPS +# print(dXMaxPS) + +# now we can add these three portions together +xInterpPS = [0] +tmp = dX1PS +for i in range(nStretch1PS): + xInterpPS.append(xInterpPS[-1] + tmp) + tmp = tmp * Alpha1PS +for i in range(nXConstPS): + xInterpPS.append(xInterpPS[-1] + dXMaxPS) +tmp = dX2PS * (Alpha2PS ** (nStretch2PS - 1)) +for i in range(nStretch2PS): + xInterpPS.append(xInterpPS[-1] + tmp) + tmp /= Alpha2PS +# print xInterpPS +# Finally, we interpolate the refined stretch stuff +c1PS = Curve(x=xPS, y=yPS, z=zPS, k=3) +XPS = c1PS(xInterpPS) +c2PS = Curve(X=XPS, k=3) +x1PS = c2PS.X[:, 0] +y1PS = c2PS.X[:, 1] + + +# ------------SS +# first compute how many stretching points do we need for both ends +tmp = dX1SS +for i in range(1000): + if tmp > dXMaxSS: + nStretch1SS = i + break + else: + tmp = tmp * Alpha1SS +# print (nStretch1SS) +tmp = dX2SS +for i in range(1000): + if tmp > dXMaxSS: + nStretch2SS = i + break + else: + tmp = tmp * Alpha2SS +# print (nStretch2SS) + +# now compute how much length does these two stretching end and the constant portion have +xLSSConst = xSS[-1] +xLSS1 = 0 +xLSS2 = 0 +for i in range(nStretch1SS): + xLSS1 += dX1SS * (Alpha1SS ** i) + xLSSConst -= dX1SS * (Alpha1SS ** i) +for i in range(nStretch2SS): + xLSS2 += dX2SS * (Alpha2SS ** i) + xLSSConst -= dX2SS * (Alpha2SS ** i) +# print xLSS1,xLSS2,xLSSConst +# update dXMax, we just want to make sure these three portion add up +nXConstSS = int(xLSSConst / dXMaxSS) +dXMaxSS = xLSSConst / nXConstSS +# print(dXMaxSS) + +# now we can add these three portions together +xInterpSS = [0] +tmp = dX1SS +for i in range(nStretch1SS): + xInterpSS.append(xInterpSS[-1] + tmp) + tmp = tmp * Alpha1SS +for i in range(nXConstSS): + xInterpSS.append(xInterpSS[-1] + dXMaxSS) +tmp = dX2SS * (Alpha2SS ** (nStretch2SS - 1)) +for i in range(nStretch2SS): + xInterpSS.append(xInterpSS[-1] + tmp) + tmp /= Alpha2SS +# print xInterpSS +# Finally, we interpolate the refined stretch stuff +c1SS = Curve(x=xSS, y=ySS, z=zSS, k=3) +XSS = c1SS(xInterpSS) +c2SS = Curve(X=XSS, k=3) +x1SS = c2SS.X[:, 0] +y1SS = c2SS.X[:, 1] + +# Since the TE is open we need to close it. Close it multiple linear segments. +delta_y = numpy.linspace(y1PS[-1], y1SS[-1], NpTE, "d") +delta_y = delta_y[1:] +delta_x = numpy.ones_like(delta_y, "d") +for i in range(len(delta_x)): + delta_x[i] = x1SS[-1] + +x1SS_Flip = x1SS[::-1] # reverse the array #numpy.flip(x1SS,axis=0) +xAll = numpy.append(x1SS_Flip, x1PS[1:]) +xAll = numpy.append(xAll, delta_x) + +y1SS_Flip = y1SS[::-1] # reverse the array # numpy.flip(y1SS,axis=0) +yAll = numpy.append(y1SS_Flip, y1PS[1:]) +yAll = numpy.append(yAll, delta_y) + + +# print mesh statistics +print("nPoints for PS: ", nStretch1PS + nStretch2PS + nXConstPS) +print("nPoints for SS: ", nStretch1SS + nStretch2SS + nXConstSS) +print("nPoints for TE: ", NpTE) +print("nPoints Total: ", nStretch1PS + nStretch2PS + nXConstPS + nStretch1SS + nStretch2SS + nXConstSS + NpTE) +print( + "Mesh cells: ", + (nStretch1PS + nStretch2PS + nXConstPS + nStretch1SS + nStretch2SS + nXConstSS + NpTE - 1) + * (NpExtrude - 1) + * (nSpan - 1), +) + +# plt.plot(xPS,yPS,'-k',linewidth=1) +# plt.plot(xAll,yAll,'ro',markersize=2) +# plt.plot(xSS,ySS,'-k',linewidth=1) +# plt.gca().set_aspect('equal', adjustable='box') +# plt.show() +# plt.savefig('figure.png',bbox_inches='tight') # save the figure to file +# plt.close() # close the figure + + +# Write the plot3d input file: +f = open("surfaceMesh.xyz", "w") +f.write("1\n") +f.write("%d %d %d\n" % (len(xAll), nSpan, 1)) +for iDim in range(3): + for z in numpy.linspace(0.0, ZSpan, nSpan): + for i in range(len(xAll)): + if iDim == 0: + f.write("%20.16f\n" % xAll[i]) + elif iDim == 1: + f.write("%20.16f\n" % yAll[i]) + else: + f.write("%20.16f\n" % z) +f.close() + +options = { + # --------------------------- + # Input Parameters + # --------------------------- + "inputFile": "surfaceMesh.xyz", + "unattachedEdgesAreSymmetry": False, + "outerFaceBC": "farfield", + "autoConnect": True, + "BC": {1: {"jLow": "zSymm", "jHigh": "zSymm"}}, + "families": "wall", + # --------------------------- + # Grid Parameters + # --------------------------- + "N": NpExtrude, + "s0": yWall, + "marchDist": marchDist, + # --------------------------- + # Pseudo Grid Parameters + # --------------------------- + "ps0": -1.0, + "pGridRatio": -1.0, + "cMax": 1.0, + # --------------------------- + # Smoothing parameters + # --------------------------- + "epsE": 2.0, + "epsI": 4.0, + "theta": 2.0, + "volCoef": 0.20, + "volBlend": 0.0005, + "volSmoothIter": 20, +} + + +hyp = pyHyp(options=options) +hyp.run() +# hyp.writeCGNS('volumeMesh.cgns') +hyp.writePlot3D("volumeMesh.xyz") diff --git a/RAE2822_Airfoil_HISA/paraview.foam b/RAE2822_Airfoil_HISA/paraview.foam new file mode 100644 index 0000000..e69de29 diff --git a/RAE2822_Airfoil_HISA/preProcessing.sh b/RAE2822_Airfoil_HISA/preProcessing.sh new file mode 100644 index 0000000..2d4744c --- /dev/null +++ b/RAE2822_Airfoil_HISA/preProcessing.sh @@ -0,0 +1,20 @@ +#!/bin/bash + +if [ -z "$WM_PROJECT" ]; then + echo "OpenFOAM environment not found, forgot to source the OpenFOAM bashrc?" + exit +fi + +# pre-processing + +# generate mesh +echo "Generating mesh.." +python genAirFoilMesh.py &> logMeshGeneration.txt +plot3dToFoam -noBlank volumeMesh.xyz >> logMeshGeneration.txt +autoPatch 45 -overwrite >> logMeshGeneration.txt +createPatch -overwrite >> logMeshGeneration.txt +renumberMesh -overwrite >> logMeshGeneration.txt +echo "Generating mesh.. Done!" + +# copy initial and boundary condition files +cp -r 0_orig 0 diff --git a/RAE2822_Airfoil_HISA/profiles/RAE2822PS.profile b/RAE2822_Airfoil_HISA/profiles/RAE2822PS.profile new file mode 100644 index 0000000..bddfd79 --- /dev/null +++ b/RAE2822_Airfoil_HISA/profiles/RAE2822PS.profile @@ -0,0 +1,63 @@ + 0.000000 0.000000 + 0.000602 -.003160 + 0.002408 -.006308 + 0.005412 -.009443 + 0.009607 -.012559 + 0.014984 -.015649 + 0.021530 -.018707 + 0.029228 -.021722 + 0.038060 -.024685 + 0.048005 -.027586 + 0.059039 -.030416 + 0.071136 -.033170 + 0.084265 -.035843 + 0.098396 -.038431 + 0.113495 -.040929 + 0.129524 -.043326 + 0.146447 -.045610 + 0.164221 -.047773 + 0.182803 -.049805 + 0.202150 -.051694 + 0.222215 -.053427 + 0.242949 -.054994 + 0.264302 -.056376 + 0.286222 -.057547 + 0.308658 -.058459 + 0.331555 -.059046 + 0.354858 -.059236 + 0.378510 -.058974 + 0.402455 -.058224 + 0.426635 -.056979 + 0.450991 -.055257 + 0.475466 -.053099 + 0.500000 -.050563 + 0.524534 -.047719 + 0.549009 -.044642 + 0.573365 -.041397 + 0.597545 -.038043 + 0.621490 -.034631 + 0.645142 -.031207 + 0.668445 -.027814 + 0.691342 -.024495 + 0.713778 -.021289 + 0.735698 -.018232 + 0.757051 -.015357 + 0.777785 -.012690 + 0.797850 -.010244 + 0.817197 -.008027 + 0.835779 -.006048 + 0.853553 -.004314 + 0.870476 -.002829 + 0.886505 -.001592 + 0.901604 -.000600 + 0.915735 0.000157 + 0.928864 0.000694 + 0.940961 0.001033 + 0.951995 0.001197 + 0.961940 0.001212 + 0.970772 0.001112 + 0.978470 0.000935 + 0.985016 0.000719 + 0.990393 0.000497 + 0.994588 0.000296 + 0.997592 0.000137 \ No newline at end of file diff --git a/RAE2822_Airfoil_HISA/profiles/RAE2822SS.profile b/RAE2822_Airfoil_HISA/profiles/RAE2822SS.profile new file mode 100644 index 0000000..bc472c9 --- /dev/null +++ b/RAE2822_Airfoil_HISA/profiles/RAE2822SS.profile @@ -0,0 +1,63 @@ + 0.000000 0.000000 + 0.000602 0.003165 + 0.002408 0.006306 + 0.005412 0.009416 + 0.009607 0.012480 + 0.014984 0.015489 + 0.021530 0.018441 + 0.029228 0.021348 + 0.038060 0.024219 + 0.048005 0.027062 + 0.059039 0.029874 + 0.071136 0.032644 + 0.084265 0.035360 + 0.098396 0.038011 + 0.113495 0.040585 + 0.129524 0.043071 + 0.146447 0.045457 + 0.164221 0.047729 + 0.182803 0.049874 + 0.202150 0.051885 + 0.222215 0.053753 + 0.242949 0.055470 + 0.264302 0.057026 + 0.286222 0.058414 + 0.308658 0.059629 + 0.331555 0.060660 + 0.354858 0.061497 + 0.378510 0.062133 + 0.402455 0.062562 + 0.426635 0.062779 + 0.450991 0.062774 + 0.475466 0.062530 + 0.500000 0.062029 + 0.524534 0.061254 + 0.549009 0.060194 + 0.573365 0.058845 + 0.597545 0.057218 + 0.621490 0.055344 + 0.645142 0.053258 + 0.668445 0.050993 + 0.691342 0.048575 + 0.713778 0.046029 + 0.735698 0.043377 + 0.757051 0.040641 + 0.777785 0.037847 + 0.797850 0.035017 + 0.817197 0.032176 + 0.835779 0.029347 + 0.853553 0.026554 + 0.870476 0.023817 + 0.886505 0.021153 + 0.901604 0.018580 + 0.915735 0.016113 + 0.928864 0.013769 + 0.940961 0.011562 + 0.951995 0.009508 + 0.961940 0.007622 + 0.970772 0.005915 + 0.978470 0.004401 + 0.985016 0.003092 + 0.990393 0.002001 + 0.994588 0.001137 + 0.997592 0.000510 \ No newline at end of file diff --git a/RAE2822_Airfoil_HISA/runScript.py b/RAE2822_Airfoil_HISA/runScript.py new file mode 100644 index 0000000..25cb1c1 --- /dev/null +++ b/RAE2822_Airfoil_HISA/runScript.py @@ -0,0 +1,268 @@ +#!/usr/bin/env python +""" +DAFoam run script for the RAE2822 airfoil at transonic conditions +""" + +# ============================================================================= +# Imports +# ============================================================================= +import os +import argparse +import numpy as np +from mpi4py import MPI +import openmdao.api as om +from mphys.multipoint import Multipoint +from dafoam.mphys import DAFoamBuilder, OptFuncs +from mphys.scenario_aerodynamic import ScenarioAerodynamic +from pygeo.mphys import OM_DVGEOCOMP +from pygeo import geo_utils + + +parser = argparse.ArgumentParser() +# which optimizer to use. Options are: IPOPT (default), SLSQP, and SNOPT +parser.add_argument("-optimizer", help="optimizer to use", type=str, default="IPOPT") +# which task to run. Options are: run_driver (default), run_model, compute_totals, check_totals +parser.add_argument("-task", help="type of run to do", type=str, default="run_model") +args = parser.parse_args() + +# ============================================================================= +# Input Parameters +# ============================================================================= + +# NOTE: Input Parameters Must Also be Specified in "0/include/freestreamconditions" +U0 = 233.603 +p0 = 108988.0 +T0 = 255.56 +nuTilda0 = 4.5e-5 +CL_target = 0.7 +twist0 = 2.52517169 +A0 = 0.01 +# rho is used for normalizing CD and CL +rho0 = p0 / T0 / 287 + +# Input parameters for DAFoam +daOptions = { + "designSurfaces": ["wing"], + "solverName": "DAHisaFoam", + # "primalInitCondition": {"U": [Ux, Uy, 0.0], "p": p0, "T": T0}, + # "primalBC": { + # "U0": {"variable": "U", "patches": ["inout"], "value": [Ux, Uy, 0.0]}, + # "p0": {"variable": "p", "patches": ["inout"], "value": [p0]}, + # "T0": {"variable": "T", "patches": ["inout"], "value": [T0]}, + # "useWallFunction": True, + # }, + "function": { + "CD": { + "type": "force", + "source": "patchToFace", + "patches": ["wing"], + "directionMode": "fixedDirection", + "direction": [1.0, 0.0, 0.0], + "scale": 1.0 / (0.5 * U0 * U0 * A0 * rho0), + }, + "CL": { + "type": "force", + "source": "patchToFace", + "patches": ["wing"], + "directionMode": "fixedDirection", + "direction": [0.0, 1.0, 0.0], + "scale": 1.0 / (0.5 * U0 * U0 * A0 * rho0), + }, + }, + "adjStateOrdering": "cell", + "adjEqnOption": { + "gmresRelTol": 1.0e-6, + "pcFillLevel": 1, + "jacMatReOrdering": "natural", + "gmresMaxIters": 2000, + "gmresRestart": 2000, + }, + "normalizeStates": { + "U": U0, + "p": p0, + "T": T0, + "nuTilda": nuTilda0 * 10.0, + }, + "checkMeshThreshold": {"maxNonOrth": 70.0, "maxSkewness": 6.0, "maxAspectRatio": 5000.0}, + "inputInfo": { + "aero_vol_coords": {"type": "volCoord", "components": ["solver", "function"]}, + }, +} + +# Mesh deformation setup +meshOptions = { + "gridFile": os.getcwd(), + "fileType": "OpenFOAM", + "useRotations": False, + # point and normal for the symmetry plane + "symmetryPlanes": [[[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], [[0.0, 0.0, 0.01], [0.0, 0.0, 1.0]]], +} + + +# Top class to setup the optimization problem +class Top(Multipoint): + def setup(self): + + # create the builder to initialize the DASolvers + dafoam_builder = DAFoamBuilder(daOptions, meshOptions, scenario="aerodynamic") + dafoam_builder.initialize(self.comm) + + # add the design variable component to keep the top level design variables + self.add_subsystem("dvs", om.IndepVarComp(), promotes=["*"]) + + # add the mesh component + self.add_subsystem("mesh", dafoam_builder.get_mesh_coordinate_subsystem()) + + # add the geometry component (FFD) + self.add_subsystem("geometry", OM_DVGEOCOMP(file="FFD/wingFFD.xyz", type="ffd")) + + # add a scenario (flow condition) for optimization, we pass the builder + # to the scenario to actually run the flow and adjoint + self.mphys_add_scenario("scenario1", ScenarioAerodynamic(aero_builder=dafoam_builder)) + + # need to manually connect the x_aero0 between the mesh and geometry components + # here x_aero0 means the surface coordinates of structurally undeformed mesh + self.connect("mesh.x_aero0", "geometry.x_aero_in") + # need to manually connect the x_aero0 between the geometry component and the scenario1 + # scenario group + self.connect("geometry.x_aero0", "scenario1.x_aero") + + def configure(self): + + # get the surface coordinates from the mesh component + points = self.mesh.mphys_get_surface_mesh() + + # add pointset to the geometry component + self.geometry.nom_add_discipline_coords("aero", points) + + # set the triangular points to the geometry component for geometric constraints + tri_points = self.mesh.mphys_get_triangulated_surface() + self.geometry.nom_setConstraintSurface(tri_points) + + # Create reference axis for the twist variable + self.geometry.nom_addRefAxis(name="wingAxis", xFraction=0.25, alignIndex="k") + + # use the shape function to define shape variables for 2D airfoil + pts = self.geometry.DVGeo.getLocalIndex(0) + dir_y = np.array([0.0, 1.0, 0.0]) + shapes = [] + for i in range(1, pts.shape[0] - 1): + for j in range(pts.shape[1]): + # k=0 and k=1 move together to ensure symmetry + shapes.append({pts[i, j, 0]: dir_y, pts[i, j, 1]: dir_y}) + # LE/TE shape, the j=0 and j=1 move in opposite directions so that + # the LE/TE are fixed + for i in [0, pts.shape[0] - 1]: + shapes.append({pts[i, 0, 0]: dir_y, pts[i, 0, 1]: dir_y, pts[i, 1, 0]: -dir_y, pts[i, 1, 1]: -dir_y}) + self.geometry.nom_addShapeFunctionDV(dvName="shape", shapes=shapes) + + # Set up global design variables. We dont change the root twist + def twist(val, geo): + for i in range(2): + geo.rot_z["wingAxis"].coef[i] = -val[0] + + self.geometry.nom_addGlobalDV(dvName="twist", value=np.ones(1) * twist0, func=twist) + + # setup the volume and thickness constraints + leList = [[1e-3, 0.0, 1e-3], [1e-3, 0.0, 0.01 - 1e-3]] + teList = [[0.961, 0.005, 1e-3], [0.961, 0.005, 0.01 - 1e-3]] + self.geometry.nom_addThicknessConstraints2D("thickcon", leList, teList, nSpan=2, nChord=10) + self.geometry.nom_addVolumeConstraint("volcon", leList, teList, nSpan=2, nChord=10) + self.geometry.nom_addLERadiusConstraints("rcon", leList, 2, [0.0, 1.0, 0.0], [-1.0, 0.0, 0.0]) + # NOTE: we no longer need to define the sym and LE/TE constraints + # because these constraints are defined in the above shape function + + # add the design variables to the dvs component's output + self.dvs.add_output("shape", val=np.array([0] * len(shapes))) + self.dvs.add_output("twist", val=np.ones(1) * twist0) + + # manually connect the dvs output to the geometry and scenario1 + + self.connect("shape", "geometry.shape") + self.connect("twist", "geometry.twist") + + # define the design variables to the top level + self.add_design_var("shape", lower=-1.0, upper=1.0, scaler=10.0) + self.add_design_var("twist", lower=-10.0, upper=10.0, scaler=0.1) + + # add objective and constraints to the top level + self.add_objective("scenario1.aero_post.CD", scaler=1.0) + self.add_constraint("scenario1.aero_post.CL", equals=CL_target, scaler=1.0) + self.add_constraint("geometry.thickcon", lower=0.5, upper=3.0, scaler=1.0) + self.add_constraint("geometry.volcon", lower=1.0, scaler=1.0) + self.add_constraint("geometry.rcon", lower=0.8, scaler=1.0) + + +# OpenMDAO setup +prob = om.Problem() +prob.model = Top() +prob.setup(mode="rev") +om.n2(prob, show_browser=False, outfile="mphys.html") + +# initialize the optimization function +optFuncs = OptFuncs(daOptions, prob) + +# use pyoptsparse to setup optimization +prob.driver = om.pyOptSparseDriver() +prob.driver.options["optimizer"] = args.optimizer +# options for optimizers +if args.optimizer == "SNOPT": + prob.driver.opt_settings = { + "Major feasibility tolerance": 1.0e-5, + "Major optimality tolerance": 1.0e-5, + "Minor feasibility tolerance": 1.0e-5, + "Verify level": -1, + "Function precision": 1.0e-5, + "Major iterations limit": 100, + "Nonderivative linesearch": None, + "Print file": "opt_SNOPT_print.txt", + "Summary file": "opt_SNOPT_summary.txt", + } +elif args.optimizer == "IPOPT": + prob.driver.opt_settings = { + "tol": 1.0e-5, + "constr_viol_tol": 1.0e-5, + "max_iter": 100, + "print_level": 5, + "output_file": "opt_IPOPT.txt", + "mu_strategy": "adaptive", + "limited_memory_max_history": 10, + "nlp_scaling_method": "none", + "alpha_for_y": "full", + "recalc_y": "yes", + } +elif args.optimizer == "SLSQP": + prob.driver.opt_settings = { + "ACC": 1.0e-5, + "MAXIT": 100, + "IFILE": "opt_SLSQP.txt", + } +else: + print("optimizer arg not valid!") + exit(1) + +prob.driver.options["debug_print"] = ["nl_cons", "objs", "desvars"] +prob.driver.options["print_opt_prob"] = True +prob.driver.hist_file = "OptView.hst" + +if args.task == "run_driver": + # solve CL + optFuncs.findFeasibleDesign(["scenario1.aero_post.CL"], ["twist"], targets=[CL_target], designVarsComp=[0], epsFD=[1e-1]) + # run the optimization + prob.run_driver() +elif args.task == "run_model": + # just run the primal once + prob.run_model() +elif args.task == "compute_totals": + # just run the primal and adjoint once + prob.run_model() + totals = prob.compute_totals() + if MPI.COMM_WORLD.rank == 0: + print(totals) +elif args.task == "check_totals": + # verify the total derivatives against the finite-difference + prob.run_model() + prob.check_totals(compact_print=False, step=1e-3, form="central", step_calc="abs") +else: + print("task arg not found!") + exit(1) diff --git a/RAE2822_Airfoil_HISA/system/controlDict b/RAE2822_Airfoil_HISA/system/controlDict new file mode 100644 index 0000000..9ee733e --- /dev/null +++ b/RAE2822_Airfoil_HISA/system/controlDict @@ -0,0 +1,39 @@ +/*--------------------------------*- C++ -*---------------------------------*\ +| ======== | | +| \ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \ / O peration | Version: v1812 | +| \ / A nd | Web: www.OpenFOAM.com | +| \/ M anipulation | | +\*--------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object controlDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + +application hisa; +startFrom latestTime; +startTime 0; +stopAt endTime; +endTime 8000; +deltaT 1; +writeControl timeStep; +writeInterval 8000; +purgeWrite 0; +writeFormat ascii; +writePrecision 16; +writeCompression on; +timeFormat general; +timePrecision 8; +runTimeModifiable false; +HiSAPrint 0; + +DebugSwitches +{ + SolverPerformance 0; +} diff --git a/RAE2822_Airfoil_HISA/system/createPatchDict b/RAE2822_Airfoil_HISA/system/createPatchDict new file mode 100644 index 0000000..5e2c407 --- /dev/null +++ b/RAE2822_Airfoil_HISA/system/createPatchDict @@ -0,0 +1,82 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: v1812 | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object createPatchDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +pointSync false; + +// Patches to create. +patches +( + { + // Name of new patch + name symmetry1; + + // Dictionary to construct new patch from + patchInfo + { + type symmetry; + } + + // How to construct: either from 'patches' or 'set' + constructFrom patches; + patches (auto0); + } + { + // Name of new patch + name symmetry2; + + // Dictionary to construct new patch from + patchInfo + { + type symmetry; + } + + // How to construct: either from 'patches' or 'set' + constructFrom patches; + patches (auto1); + } + + { + // Name of new patch + name wing; + + // Dictionary to construct new patch from + patchInfo + { + type wall; + } + + // How to construct: either from 'patches' or 'set' + constructFrom patches; + patches (auto2 auto3); + } + + { + // Name of new patch + name inout; + + // Dictionary to construct new patch from + patchInfo + { + type patch; + } + + // How to construct: either from 'patches' or 'set' + constructFrom patches; + patches (auto4); + } +); + +// ************************************************************************* // diff --git a/RAE2822_Airfoil_HISA/system/decomposeParDict b/RAE2822_Airfoil_HISA/system/decomposeParDict new file mode 100644 index 0000000..3f04af2 --- /dev/null +++ b/RAE2822_Airfoil_HISA/system/decomposeParDict @@ -0,0 +1,32 @@ +/*--------------------------------*- C++ -*---------------------------------*\ +| ======== | | +| \ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \ / O peration | Version: v1812 | +| \ / A nd | Web: www.OpenFOAM.com | +| \/ M anipulation | | +\*--------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object decomposeParDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +numberOfSubdomains 64; + +method scotch; + +simpleCoeffs +{ + n (2 2 1); + delta 0.001; +} + +distributed false; + +roots(); + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/RAE2822_Airfoil_HISA/system/fvSchemes b/RAE2822_Airfoil_HISA/system/fvSchemes new file mode 100644 index 0000000..f39c0a3 --- /dev/null +++ b/RAE2822_Airfoil_HISA/system/fvSchemes @@ -0,0 +1,69 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| | +| HiSA: High Speed Aerodynamic solver | +| Copyright (C) 2014-2017 Johan Heyns - CSIR, South Africa | +| Copyright (C) 2014-2017 Oliver Oxtoby - CSIR, South Africa | +| | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object fvSchemes; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +fluxScheme JST; // AUSMPlusUp; +jst_k2 0.5; +jst_k4 0.02; +sensorName pressure; +lowMachAusm false; + +ddtSchemes +{ + default bounded dualTime rPseudoDeltaT steadyState; +} + +gradSchemes +{ + default faceLeastSquares linear; + grad(nuTilda) Gauss linear; +} + +divSchemes +{ + default none; + div(tauMC) Gauss linear; + div(phi,nuTilda) bounded Gauss upwind; + div(pc) bounded Gauss upwind; +} + +laplacianSchemes +{ + default Gauss linear corrected; + laplacian(muEff,U) Gauss linear compact; + laplacian(alphaEff,e) Gauss linear compact; +} + +interpolate wVanLeer; +interpolationSchemes +{ + default linear; + reconstruct(rho) $interpolate; + reconstruct(U) $interpolate; + reconstruct(T) $interpolate; +} + +snGradSchemes +{ + default corrected; +} + +wallDist +{ + method meshWave; +} + +// ************************************************************************* // diff --git a/RAE2822_Airfoil_HISA/system/fvSolution b/RAE2822_Airfoil_HISA/system/fvSolution new file mode 100644 index 0000000..8b72658 --- /dev/null +++ b/RAE2822_Airfoil_HISA/system/fvSolution @@ -0,0 +1,74 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| | +| HiSA: High Speed Aerodynamic solver | +| Copyright (C) 2014-2017 Johan Heyns - CSIR, South Africa | +| Copyright (C) 2014-2017 Oliver Oxtoby - CSIR, South Africa | +| | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object fvSolution; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +solvers +{ + nuTilda + { + solver smoothSolver; + smoother symGaussSeidel; + tolerance 0; + relTol 1e-1; + minIter 1; + } + + yPsi + { + solver GAMG; + smoother GaussSeidel; + cacheAgglomeration true; + nCellsInCoarsestLevel 10; + agglomerator faceAreaPair; + mergeLevels 1; + tolerance 0; + relTol 0; + } +} + +relaxationFactors +{ + equations + { + nuTilda 0.5; + yWall 0.7; + } +} + +flowSolver +{ + solver GMRES; + GMRES + { + inviscidJacobian LaxFriedrich; + viscousJacobian laplacian; + preconditioner LUSGS; + + maxIter 20; + nKrylov 8; + solverTolRel 1e-1 (1e-1 1e-1 1e-1) 1e-1; + } +} + +pseudoTime +{ + pseudoTol 1e-8 (1e-8 1e-8 1e-8) 1e-8; + pseudoCoNum 1.0; + pseudoCoNumMax 10.0; + localTimestepping true; +} + +// ************************************************************************* // From f1117df042ca611e511823f5766967edc6001573 Mon Sep 17 00:00:00 2001 From: nmadison11 Date: Thu, 4 Dec 2025 10:13:45 -0600 Subject: [PATCH 2/2] Simplified Runscript --- RAE2822_Airfoil_HISA/FFD/genFFD.py | 99 ----------- RAE2822_Airfoil_HISA/FFD/wingFFD.xyz | 5 - RAE2822_Airfoil_HISA/runScript.py | 240 ++++----------------------- 3 files changed, 35 insertions(+), 309 deletions(-) delete mode 100644 RAE2822_Airfoil_HISA/FFD/genFFD.py delete mode 100644 RAE2822_Airfoil_HISA/FFD/wingFFD.xyz diff --git a/RAE2822_Airfoil_HISA/FFD/genFFD.py b/RAE2822_Airfoil_HISA/FFD/genFFD.py deleted file mode 100644 index b35e34d..0000000 --- a/RAE2822_Airfoil_HISA/FFD/genFFD.py +++ /dev/null @@ -1,99 +0,0 @@ -import numpy as np -import sys - -def writeFFDFile(fileName,nBlocks,nx,ny,nz,points): - ''' - Take in a set of points and write the plot 3dFile - ''' - - f = open(fileName,'w') - - f.write('%d\n'%nBlocks) - for i in range(nBlocks): - f.write('%d %d %d '%(nx[i],ny[i],nz[i])) - # end - f.write('\n') - for block in range(nBlocks): - for k in range(nz[block]): - for j in range(ny[block]): - for i in range(nx[block]): - f.write('%f '%points[block][i,j,k,0]) - # end - # end - # end - f.write('\n') - - for k in range(nz[block]): - for j in range(ny[block]): - for i in range(nx[block]): - f.write('%f '%points[block][i,j,k,1]) - # end - # end - # end - f.write('\n') - - for k in range(nz[block]): - for j in range(ny[block]): - for i in range(nx[block]): - f.write('%f '%points[block][i,j,k,2]) - # end - # end - # end - # end - f.close() - return - -def returnBlockPoints(corners,nx,ny,nz): - ''' - corners needs to be 8 x 3 - ''' - points = np.zeros([nx,ny,nz,3]) - - # points 1 - 4 are the iMin face - # points 5 - 8 are the iMax face - - for idim in range(3): - edge1 = np.linspace(corners[0][idim],corners[4][idim],nx) - edge2 = np.linspace(corners[1][idim],corners[5][idim],nx) - edge3 = np.linspace(corners[2][idim],corners[6][idim],nx) - edge4 = np.linspace(corners[3][idim],corners[7][idim],nx) - - for i in range(nx): - edge5 = np.linspace(edge1[i],edge3[i],ny) - edge6 = np.linspace(edge2[i],edge4[i],ny) - for j in range(ny): - edge7 = np.linspace(edge5[j],edge6[j],nz) - points[i,j,:,idim] = edge7 - # end - # end - # end - - return points - -################ FFD ############## -nBlocks = 1 - -nx = [8] -ny = [2] -nz = [2] - -corners = np.zeros([nBlocks,8,3]) - -corners[0,0,:] = [-0.010,-.0700,0.0] -corners[0,1,:] = [-0.010,-0.0700,0.01] -corners[0,2,:] = [-0.010,0.0700,0.0] -corners[0,3,:] = [-0.010,0.0700,0.01] -corners[0,4,:] = [ 1.01,-.07,0.0] -corners[0,5,:] = [ 1.01,-.07,0.010] -corners[0,6,:] = [ 1.01,0.07,0.0] -corners[0,7,:] = [ 1.01,0.07,0.01] - - -points = [] -for block in range(nBlocks): - points.append(returnBlockPoints(corners[block],nx[block],ny[block],nz[block])) - -#print points -fileName = 'wingFFD.xyz' -writeFFDFile(fileName,nBlocks,nx,ny,nz,points) - diff --git a/RAE2822_Airfoil_HISA/FFD/wingFFD.xyz b/RAE2822_Airfoil_HISA/FFD/wingFFD.xyz deleted file mode 100644 index 26f0558..0000000 --- a/RAE2822_Airfoil_HISA/FFD/wingFFD.xyz +++ /dev/null @@ -1,5 +0,0 @@ -1 -8 2 2 --0.010000 0.135714 0.281429 0.427143 0.572857 0.718571 0.864286 1.010000 -0.010000 0.135714 0.281429 0.427143 0.572857 0.718571 0.864286 1.010000 -0.010000 0.135714 0.281429 0.427143 0.572857 0.718571 0.864286 1.010000 -0.010000 0.135714 0.281429 0.427143 0.572857 0.718571 0.864286 1.010000 --0.070000 -0.070000 -0.070000 -0.070000 -0.070000 -0.070000 -0.070000 -0.070000 0.070000 0.070000 0.070000 0.070000 0.070000 0.070000 0.070000 0.070000 -0.070000 -0.070000 -0.070000 -0.070000 -0.070000 -0.070000 -0.070000 -0.070000 0.070000 0.070000 0.070000 0.070000 0.070000 0.070000 0.070000 0.070000 -0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000 \ No newline at end of file diff --git a/RAE2822_Airfoil_HISA/runScript.py b/RAE2822_Airfoil_HISA/runScript.py index 25cb1c1..bbebf99 100644 --- a/RAE2822_Airfoil_HISA/runScript.py +++ b/RAE2822_Airfoil_HISA/runScript.py @@ -16,20 +16,15 @@ from mphys.scenario_aerodynamic import ScenarioAerodynamic from pygeo.mphys import OM_DVGEOCOMP from pygeo import geo_utils +from dafoam import PYDAFOAM -parser = argparse.ArgumentParser() -# which optimizer to use. Options are: IPOPT (default), SLSQP, and SNOPT -parser.add_argument("-optimizer", help="optimizer to use", type=str, default="IPOPT") -# which task to run. Options are: run_driver (default), run_model, compute_totals, check_totals -parser.add_argument("-task", help="type of run to do", type=str, default="run_model") -args = parser.parse_args() - # ============================================================================= # Input Parameters # ============================================================================= # NOTE: Input Parameters Must Also be Specified in "0/include/freestreamconditions" + U0 = 233.603 p0 = 108988.0 T0 = 255.56 @@ -37,20 +32,40 @@ CL_target = 0.7 twist0 = 2.52517169 A0 = 0.01 +alpha0= 2.31 +k0= 8.5 +omega0= 5.6e4 +epsilon0= 4.3e4 # rho is used for normalizing CD and CL rho0 = p0 / T0 / 287 + +def calcUAndDir(UMag, alpha1): + dragDir = [float(np.cos(alpha1 * np.pi / 180)), float(np.sin(alpha1 * np.pi / 180)), 0.0] + liftDir = [float(-np.sin(alpha1 * np.pi / 180)), float(np.cos(alpha1 * np.pi / 180)), 0.0] + inletU = [float(UMag * np.cos(alpha1 * np.pi / 180)), float(UMag * np.sin(alpha1 * np.pi / 180)), 0.0] + return inletU, dragDir, liftDir + + +inletU, dragDir, liftDir = calcUAndDir(U0, alpha0) + + # Input parameters for DAFoam daOptions = { "designSurfaces": ["wing"], "solverName": "DAHisaFoam", - # "primalInitCondition": {"U": [Ux, Uy, 0.0], "p": p0, "T": T0}, - # "primalBC": { - # "U0": {"variable": "U", "patches": ["inout"], "value": [Ux, Uy, 0.0]}, - # "p0": {"variable": "p", "patches": ["inout"], "value": [p0]}, - # "T0": {"variable": "T", "patches": ["inout"], "value": [T0]}, - # "useWallFunction": True, - # }, +# "primalInitCondition": {"U": [Ux, Uy, 0.0], "p": p0, "T": T0}, + "primalMinResTol": 1.0e-7, + "primalBC": { + "U0": {"variable": "U", "patches": ["inout"], "value": inletU}, + "p0": {"variable": "p", "patches": ["inout"], "value": [p0]}, + "T0": {"variable": "T", "patches": ["inout"], "value": [T0]}, + "nuTilda0": {"variable": "nuTilda", "patches": ["inout"], "value": [nuTilda0]}, + "k0": {"variable": "k", "patches": ["inout"], "value": [k0]}, + "omega0": {"variable": "omega", "patches": ["inout"], "value": [omega0]}, + "epsilon0": {"variable": "epsilon", "patches": ["inout"], "value": [epsilon0]}, + "useWallFunction": False, + }, "function": { "CD": { "type": "force", @@ -69,200 +84,15 @@ "scale": 1.0 / (0.5 * U0 * U0 * A0 * rho0), }, }, - "adjStateOrdering": "cell", - "adjEqnOption": { - "gmresRelTol": 1.0e-6, - "pcFillLevel": 1, - "jacMatReOrdering": "natural", - "gmresMaxIters": 2000, - "gmresRestart": 2000, - }, - "normalizeStates": { - "U": U0, - "p": p0, - "T": T0, - "nuTilda": nuTilda0 * 10.0, - }, "checkMeshThreshold": {"maxNonOrth": 70.0, "maxSkewness": 6.0, "maxAspectRatio": 5000.0}, - "inputInfo": { - "aero_vol_coords": {"type": "volCoord", "components": ["solver", "function"]}, - }, -} -# Mesh deformation setup -meshOptions = { - "gridFile": os.getcwd(), - "fileType": "OpenFOAM", - "useRotations": False, - # point and normal for the symmetry plane - "symmetryPlanes": [[[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], [[0.0, 0.0, 0.01], [0.0, 0.0, 1.0]]], } -# Top class to setup the optimization problem -class Top(Multipoint): - def setup(self): - - # create the builder to initialize the DASolvers - dafoam_builder = DAFoamBuilder(daOptions, meshOptions, scenario="aerodynamic") - dafoam_builder.initialize(self.comm) - - # add the design variable component to keep the top level design variables - self.add_subsystem("dvs", om.IndepVarComp(), promotes=["*"]) - - # add the mesh component - self.add_subsystem("mesh", dafoam_builder.get_mesh_coordinate_subsystem()) - - # add the geometry component (FFD) - self.add_subsystem("geometry", OM_DVGEOCOMP(file="FFD/wingFFD.xyz", type="ffd")) - - # add a scenario (flow condition) for optimization, we pass the builder - # to the scenario to actually run the flow and adjoint - self.mphys_add_scenario("scenario1", ScenarioAerodynamic(aero_builder=dafoam_builder)) - - # need to manually connect the x_aero0 between the mesh and geometry components - # here x_aero0 means the surface coordinates of structurally undeformed mesh - self.connect("mesh.x_aero0", "geometry.x_aero_in") - # need to manually connect the x_aero0 between the geometry component and the scenario1 - # scenario group - self.connect("geometry.x_aero0", "scenario1.x_aero") - - def configure(self): - - # get the surface coordinates from the mesh component - points = self.mesh.mphys_get_surface_mesh() - - # add pointset to the geometry component - self.geometry.nom_add_discipline_coords("aero", points) - - # set the triangular points to the geometry component for geometric constraints - tri_points = self.mesh.mphys_get_triangulated_surface() - self.geometry.nom_setConstraintSurface(tri_points) - - # Create reference axis for the twist variable - self.geometry.nom_addRefAxis(name="wingAxis", xFraction=0.25, alignIndex="k") - - # use the shape function to define shape variables for 2D airfoil - pts = self.geometry.DVGeo.getLocalIndex(0) - dir_y = np.array([0.0, 1.0, 0.0]) - shapes = [] - for i in range(1, pts.shape[0] - 1): - for j in range(pts.shape[1]): - # k=0 and k=1 move together to ensure symmetry - shapes.append({pts[i, j, 0]: dir_y, pts[i, j, 1]: dir_y}) - # LE/TE shape, the j=0 and j=1 move in opposite directions so that - # the LE/TE are fixed - for i in [0, pts.shape[0] - 1]: - shapes.append({pts[i, 0, 0]: dir_y, pts[i, 0, 1]: dir_y, pts[i, 1, 0]: -dir_y, pts[i, 1, 1]: -dir_y}) - self.geometry.nom_addShapeFunctionDV(dvName="shape", shapes=shapes) - - # Set up global design variables. We dont change the root twist - def twist(val, geo): - for i in range(2): - geo.rot_z["wingAxis"].coef[i] = -val[0] - - self.geometry.nom_addGlobalDV(dvName="twist", value=np.ones(1) * twist0, func=twist) - - # setup the volume and thickness constraints - leList = [[1e-3, 0.0, 1e-3], [1e-3, 0.0, 0.01 - 1e-3]] - teList = [[0.961, 0.005, 1e-3], [0.961, 0.005, 0.01 - 1e-3]] - self.geometry.nom_addThicknessConstraints2D("thickcon", leList, teList, nSpan=2, nChord=10) - self.geometry.nom_addVolumeConstraint("volcon", leList, teList, nSpan=2, nChord=10) - self.geometry.nom_addLERadiusConstraints("rcon", leList, 2, [0.0, 1.0, 0.0], [-1.0, 0.0, 0.0]) - # NOTE: we no longer need to define the sym and LE/TE constraints - # because these constraints are defined in the above shape function - - # add the design variables to the dvs component's output - self.dvs.add_output("shape", val=np.array([0] * len(shapes))) - self.dvs.add_output("twist", val=np.ones(1) * twist0) - - # manually connect the dvs output to the geometry and scenario1 - - self.connect("shape", "geometry.shape") - self.connect("twist", "geometry.twist") - - # define the design variables to the top level - self.add_design_var("shape", lower=-1.0, upper=1.0, scaler=10.0) - self.add_design_var("twist", lower=-10.0, upper=10.0, scaler=0.1) - - # add objective and constraints to the top level - self.add_objective("scenario1.aero_post.CD", scaler=1.0) - self.add_constraint("scenario1.aero_post.CL", equals=CL_target, scaler=1.0) - self.add_constraint("geometry.thickcon", lower=0.5, upper=3.0, scaler=1.0) - self.add_constraint("geometry.volcon", lower=1.0, scaler=1.0) - self.add_constraint("geometry.rcon", lower=0.8, scaler=1.0) - - -# OpenMDAO setup -prob = om.Problem() -prob.model = Top() -prob.setup(mode="rev") -om.n2(prob, show_browser=False, outfile="mphys.html") - -# initialize the optimization function -optFuncs = OptFuncs(daOptions, prob) - -# use pyoptsparse to setup optimization -prob.driver = om.pyOptSparseDriver() -prob.driver.options["optimizer"] = args.optimizer -# options for optimizers -if args.optimizer == "SNOPT": - prob.driver.opt_settings = { - "Major feasibility tolerance": 1.0e-5, - "Major optimality tolerance": 1.0e-5, - "Minor feasibility tolerance": 1.0e-5, - "Verify level": -1, - "Function precision": 1.0e-5, - "Major iterations limit": 100, - "Nonderivative linesearch": None, - "Print file": "opt_SNOPT_print.txt", - "Summary file": "opt_SNOPT_summary.txt", - } -elif args.optimizer == "IPOPT": - prob.driver.opt_settings = { - "tol": 1.0e-5, - "constr_viol_tol": 1.0e-5, - "max_iter": 100, - "print_level": 5, - "output_file": "opt_IPOPT.txt", - "mu_strategy": "adaptive", - "limited_memory_max_history": 10, - "nlp_scaling_method": "none", - "alpha_for_y": "full", - "recalc_y": "yes", - } -elif args.optimizer == "SLSQP": - prob.driver.opt_settings = { - "ACC": 1.0e-5, - "MAXIT": 100, - "IFILE": "opt_SLSQP.txt", - } -else: - print("optimizer arg not valid!") - exit(1) - -prob.driver.options["debug_print"] = ["nl_cons", "objs", "desvars"] -prob.driver.options["print_opt_prob"] = True -prob.driver.hist_file = "OptView.hst" +DASolver = PYDAFOAM(options=daOptions, comm=MPI.COMM_WORLD) +DASolver() -if args.task == "run_driver": - # solve CL - optFuncs.findFeasibleDesign(["scenario1.aero_post.CL"], ["twist"], targets=[CL_target], designVarsComp=[0], epsFD=[1e-1]) - # run the optimization - prob.run_driver() -elif args.task == "run_model": - # just run the primal once - prob.run_model() -elif args.task == "compute_totals": - # just run the primal and adjoint once - prob.run_model() - totals = prob.compute_totals() - if MPI.COMM_WORLD.rank == 0: - print(totals) -elif args.task == "check_totals": - # verify the total derivatives against the finite-difference - prob.run_model() - prob.check_totals(compact_print=False, step=1e-3, form="central", step_calc="abs") -else: - print("task arg not found!") - exit(1) +# +#funcs = {} +#evalFuncs = ["CD", "CL"] +#DASolver.evalFunctions(funcs, evalFuncs)