L/D calculation in Python API

[Ishaan Khullar]

Hi I am trying to work out L/D ratios for specific wing geometries using the python API however I feel I am getting odd results. Could someone please look through if there are any errors

Results: 

Camber  Camber_Loc  Thickness        CL        CD        L/D

0       0         0.3       0.09  0.045286  0.001371  33.035347
1       0         0.3       0.12  0.045283  0.001371  33.027087
2       0         0.3       0.15  0.045288  0.001371  33.022597
3       0         0.4       0.09  0.045286  0.001371  33.035344
4       0         0.4       0.12  0.045283  0.001371  33.027089
5       0         0.4       0.15  0.045288  0.001371  33.022600
6       0         0.5       0.09  0.045286  0.001371  33.035344
7       0         0.5       0.12  0.045283  0.001371  33.027090
8       0         0.5       0.15  0.045288  0.001371  33.022599
9    0.02         0.3       0.09  0.062946  0.002028  31.034945
10   0.02         0.3       0.12  0.063024  0.002031  31.032436
11   0.02         0.3       0.15  0.063091  0.002034  31.015059
12   0.02         0.4       0.09  0.064161  0.002076  30.910062
13   0.02         0.4       0.12  0.064256  0.002080  30.893812
14   0.02         0.4       0.15  0.064323  0.002084  30.869482
15   0.02         0.5       0.09  0.065847  0.002146  30.679829
16   0.02         0.5       0.12  0.065892  0.002149  30.654525
17   0.02         0.5       0.15  0.065998  0.002154  30.633420
18   0.04         0.3       0.09  0.080733  0.002856  28.267882
19   0.04         0.3       0.12  0.080907  0.002869  28.197746
20   0.04         0.3       0.15  0.081109  0.002875  28.209264
21   0.04         0.4       0.09  0.083292  0.002986  27.894574
22   0.04         0.4       0.12  0.083372  0.002998  27.813695
23   0.04         0.4       0.15  0.083591  0.003006  27.808415
24   0.04         0.5       0.09  0.086530  0.003165  27.338019
25   0.04         0.5       0.12  0.086625  0.003175  27.283001
26   0.04         0.5       0.15  0.086808  0.003176  27.331718
27   0.06         0.3       0.09  0.098563  0.003890  25.334865
28   0.06         0.3       0.12  0.098854  0.003906  25.309681
29   0.06         0.3       0.15  0.099142  0.003926  25.250295
30   0.06         0.4       0.09  0.102177  0.004138  24.695113
31   0.06         0.4       0.12  0.102348  0.004161  24.599846
32   0.06         0.4       0.15  0.102663  0.004169  24.627554
33   0.06         0.5       0.09  0.106891  0.004452  24.007497
34   0.06         0.5       0.12  0.107040  0.004464  23.979841
35   0.06         0.5       0.15  0.107198  0.004472  23.970331

CODE:

import sys

import csv

import pandas as pd

import os

# Add OpenVSP Python module path

sys.path.append(r"C:\Users\khull\OneDrive\Documents\OpenVsp\OpenVSP-3.43.0-win64\python\openvsp")

import vsp

vsp.SetVSPAEROPath(r"C:\Users\khull\OneDrive\Documents\OpenVsp\OpenVSP-3.43.0-win64\python\openvsp\openvsp")

# Analysis conditions

aoa = 5.0

mach = 0.17

# Design parameter ranges

camber_list = [0, 0.02, 0.04, 0.06]

cloc_list = [0.3, 0.4, 0.5]

tc_list = [0.09, 0.12, 0.15]

# Wing constants

span = 5.5

root_chord = 1.625

tip_chord = 0.325

sref = 10.725

cref = 0.975

bref = span

recref = 4.4631e+6

df = pd.DataFrame(columns=["Camber", "Camber_Loc", "Thickness", "CL", "CD", "L/D"])

# Loop through param combinations

for m in camber_list:

    for p in cloc_list:

        for t in tc_list:

            print(f"\n Running: camber={m}%, cloc={p}, thickness={t}")

            try:

                vsp.ClearVSPModel()

                wing_id = vsp.AddGeom("WING")

                vsp.SetGeomName(wing_id, "ParametricWing")

                # Set wing shape

                print("no")

                vsp.SetParmVal(wing_id, "Sym_Planar_Flag", "Sym", vsp.SYM_XZ)

                print("HI1")

             

                vsp.SetParmVal(wing_id, "SectTess_U", "XSec_1", 20)

                vsp.SetParmVal(wing_id, "Span", "XSec_1", span)

                print("HI2")

                vsp.SetParmVal(wing_id, "Sweep", "XSec_1", 0.0)

                print("HI3")

                vsp.SetParmVal(wing_id, "Dihedral", "XSec_1", 1.0)

                print("HI4")

                vsp.SetParmVal(wing_id, "Twist", "XSec_1", -2.0)

                print("HI5")

                vsp.SetParmVal(wing_id, "Root_Chord", "XSec_1", root_chord)

                print("HI6")

                vsp.SetParmVal(wing_id, "Tip_Chord", "XSec_1", tip_chord)

                # Set airfoil type and parameters

                xsecsurf_id = vsp.GetXSecSurf(wing_id, 0)

                # Change root and tip shapes

                vsp.ChangeXSecShape(xsecsurf_id, 0, vsp.XS_FOUR_SERIES)

                vsp.ChangeXSecShape(xsecsurf_id, 1, vsp.XS_FOUR_SERIES)

                vsp.Update()

                # Get section IDs

                xsec_id_root = vsp.GetXSec(xsecsurf_id, 0)

                xsec_id_tip = vsp.GetXSec(xsecsurf_id, 1)

                # Set NACA parameters at root

                camber = vsp.GetXSecParm(xsec_id_root, "Camber")

                camber_loc = vsp.GetXSecParm(xsec_id_root, "CamberLoc")

                thickness = vsp.GetXSecParm(xsec_id_root, "ThickChord")

                vsp.SetParmVal(camber, m)

                vsp.SetParmVal(camber_loc, p)

                vsp.SetParmVal(thickness, t)

                # Set NACA parameters at tip

                camber = vsp.GetXSecParm(xsec_id_tip, "Camber")

                camber_loc = vsp.GetXSecParm(xsec_id_tip, "CamberLoc")

                thickness = vsp.GetXSecParm(xsec_id_tip, "ThickChord")

                vsp.SetParmVal(camber, m)

                vsp.SetParmVal(camber_loc, p)

                vsp.SetParmVal(thickness, t)

                vsp.Update()

                print("HI8")

               

               

                output_dir = r"C:\Users\khull\Results_Internship"

                os.makedirs(output_dir, exist_ok=True)

                # Save geometry for checking

                vsp_filename = os.path.join(output_dir, f"vspaero_geom_c{m}_cloc{p}_t{t}.vsp3")

                vsp.WriteVSPFile(vsp_filename, vsp.SET_ALL)

                print(f" Saved geometry to {vsp_filename}")

                # Save VSPAERO analysis run file

                output_dir = r"C:\Users\khull\OpenVSP-3.43.0-win64\python\openvsp\openvsp"

                os.makedirs(output_dir, exist_ok=True)

               

               

                vsp_run_filename = os.path.join(output_dir, "vspaero_run.vsp3")

                vsp.SetVSP3FileName(vsp_run_filename)

                vsp.WriteVSPFile(vsp_run_filename, vsp.SET_ALL)

                print(f" Saved VSPAERO run file to {vsp_run_filename}")

                # === VSPAERO Compute Geometry ===

                try:

                    compgeom_name = "VSPAEROComputeGeometry"

                    vsp.SetAnalysisInputDefaults(compgeom_name)

                    method = list(vsp.GetIntAnalysisInput(compgeom_name, "AnalysisMethod"))

                    method[0] = vsp.VORTEX_LATTICE

                    vsp.SetIntAnalysisInput(compgeom_name, "AnalysisMethod", method)

                    vsp.ExecAnalysis(compgeom_name)

                except:

                    print("HI10")

                # === VSPAERO Sweep ===

                analysis_name = "VSPAEROSweep"

                vsp.SetAnalysisInputDefaults(analysis_name)

                vsp.SetIntAnalysisInput(analysis_name, "GeomSet", [0])

                vsp.SetIntAnalysisInput(analysis_name, "RefFlag", [1]) #set int is for integers set double is for decimal points

                vsp.SetDoubleAnalysisInput(analysis_name, "AlphaStart", [aoa])

                vsp.SetDoubleAnalysisInput(analysis_name, "MachStart", [mach])

                vsp.SetDoubleAnalysisInput(analysis_name, "AlphaEnd", [aoa])

                vsp.SetIntAnalysisInput(analysis_name, "AlphaNpts", [1])

                vsp.SetDoubleAnalysisInput(analysis_name,"ReCrefStart", [4.4631e+6])

               

                vsp.Update()

                result_id = vsp.ExecAnalysis(analysis_name)

                # Get and save results

                try:

                    print(" Available results:")

                    history_res = vsp.FindLatestResultsID("VSPAERO_History")

                    print(history_res)

                    cl_vals = vsp.GetDoubleResults(history_res, "CL")

                    cd_vals = vsp.GetDoubleResults(history_res, "CDi")

               

                    cdo_vals = vsp.GetDoubleResults(history_res,"CDo")

                    print("jkasdhkhfakh",cdo_vals[-1])

                   

                    if cl_vals and cd_vals:

                        CL = cl_vals[-1]

                        CD = cd_vals[0] + cdo_vals[-1]

                        L_D = CL / CD if CD != 0 else 0

                       

                        print(f" Done: CL={CL:.4f}, CD={CD:.4f}, L/D={L_D:.2f}")

                        new_row = {"Camber": m, "Camber_Loc": p, "Thickness": t, "CL": CL, "CD": CD, "L/D": L_D}

                        df = pd.concat([df, pd.DataFrame([new_row])], ignore_index=True)    

                   

                       

                    else:

                        raise ValueError("Missing CL/CD values in VSPAERO results")

                except Exception as e:

                    print(f" Error for camber={m}, cloc={p}, thickness={t}: {e}")

            except Exception as e:

                print(f" Error for camber={m}, cloc={p}, thickness={t}: {e}")

print(df)