Error with FLOWUnsteady Simulation on NREL 5MW Wind Turbine (Spline Knot Condition Issue)

[ky45k]

Hello,
I am currently trying to run a FLOWUnsteady simulation on the NREL 5MW baseline wind turbine, and I keep running into an issue regarding spline approximation. Below is the modified code that I have been using for the simulation:

# Simulation of NREL 5-MW wind turbine at various incidence angles. propeller incidence sweep code.


import FLOWUnsteady as uns
import FLOWVLM as vlm

case_name       = "baseline"    # Name of this sweep case

save_path       = case_name                 # Where to save this sweep
output_runs     = [0]                      # Saves the VTK output of these AOAs for viz


# ----------------- GEOMETRY PARAMETERS ----------------------------------------

# Rotor geometry
rotor_file      = "SUI.csv"              # Rotor geometry
data_path = "/home/kweon/.julia/packages/FLOWUnsteady/jQESG/database"         # Path to rotor database
pitch           = 0.0                       # (deg) collective pitch of blades
CW              = false                     # Clock-wise rotation
xfoil           = false                     # Whether to run XFOIL
read_polar      = vlm.ap.read_polar2        # What polar reader to use

# Discretization
n               = 20                        # Number of blade elements per blade
r               = 1/5                       # Geometric expansion of elements

# Read radius of this rotor and number of blades
R, B            = uns.read_rotor(rotor_file; data_path=data_path)[[1,3]]

# ----------------- SIMULATION PARAMETERS --------------------------------------

# Operating conditions
AOAs            = -180:1:180                    # (deg) incidence angles to evaluate
magVinf         = 11.4                      # (m/s) rated wind speed (Table 1-1 참조)
J               = 0.9                       # Advance ratio Vinf/(nD)
RPM             = 12.1                      # Rated rotor speed (Table 5-1 참조)

rho             = 1.225                     # (kg/m^3) air density (표준값)
mu              = 1.79e-5                   # (kg/ms) air dynamic viscosity (표준값)
speedofsound    = 342.35                    # (m/s) speed of sound (표준값)

ReD             = 2*pi*RPM/60*R * rho/mu * 2*R      # Diameter-based Reynolds number
Matip           = 80 / speedofsound                 # Rated tip speed to sound speed 비율

println("""
    RPM:    $(RPM)
    Vinf:   $(magVinf) m/s
    Matip:  $(round(Matip, digits=3))
    ReD:    $(round(ReD, digits=0))
""")

# ----------------- SOLVER PARAMETERS ------------------------------------------

# Aerodynamic solver
VehicleType     = uns.UVLMVehicle           # Unsteady solver

# Time parameters
nrevs           = 1                         # Number of revolutions in simulation
nsteps_per_rev  = 1                        # Time steps per revolution
nsteps          = 1      # Number of time steps
ttot            = 1       # (s) total simulation time

# VPM particle shedding
p_per_step      = 2                         # Sheds per time step
shed_starting   = true                      # Whether to shed starting vortex
shed_unsteady   = true                      # Whether to shed vorticity from unsteady loading
max_particles   = 500000 # Maximum number of particles

# Regularization
sigma_rotor_surf= R/40                      # Rotor-on-VPM smoothing radius
lambda_vpm      = 2.125                     # VPM core overlap
sigma_vpm_overwrite = lambda_vpm * 2*pi*R/(nsteps_per_rev*p_per_step)

# Rotor solver
vlm_rlx         = 0.7                       # VLM relaxation
hubtiploss_correction = ((0.75, 10, 0.5, 0.05), (1, 1, 1, 1.0))

if VehicleType == uns.QVLMVehicle           # Mute colinear warnings if quasi-steady solver
    uns.vlm.VLMSolver._mute_warning(true)
end

# ----------------- 1) VEHICLE DEFINITION --------------------------------------
println("Generating geometry...")

# Generate rotor
data_path = "/home/kweon/.julia/packages/FLOWUnsteady/jQESG/database"
rotor = uns.generate_rotor(rotor_file; pitch=pitch,
                                        n=n, CW=CW, blade_r=r,
                                        altReD=[RPM, J, mu/rho],
                                        xfoil=xfoil,
                                        read_polar=read_polar,
                                        data_path=data_path,
                                        verbose=true,
                                        verbose_xfoil=false,
                                        plot_disc=true,
                                        );

println("Generating vehicle...")

# Generate vehicle
system = vlm.WingSystem()                   # System of all FLOWVLM objects
vlm.addwing(system, "Rotor", rotor)

rotors = [rotor];                           # Defining this rotor as its own system
rotor_systems = (rotors, );                 # All systems of rotors

wake_system = vlm.WingSystem()              # System that will shed a VPM wake
if VehicleType != uns.QVLMVehicle
    vlm.addwing(wake_system, "Rotor", rotor)
end

vehicle = VehicleType(system; rotor_systems=rotor_systems, wake_system=wake_system)

# ------------- 2) MANEUVER DEFINITION -----------------------------------------
Vvehicle(t) = zeros(3)
anglevehicle(t) = zeros(3)
RPMcontrol(t) = 1.0

angles = ()                                 # Angle of each tilting system (none)
RPMs = (RPMcontrol, )                       # RPM of each rotor system

maneuver = uns.KinematicManeuver(angles, RPMs, Vvehicle, anglevehicle)

# ------------- 3) SIMULATION DEFINITION ---------------------------------------
Vref = 0.0                                  # Reference velocity to scale maneuver by
RPMref = RPM                                # Reference RPM to scale maneuver by

Vinit = Vref*Vvehicle(0)                    # Initial vehicle velocity
Winit = pi/180*(anglevehicle(1e-6) - anglevehicle(0))/(1e-6*ttot)

simulation = uns.Simulation(vehicle, maneuver, Vref, RPMref, ttot; Vinit=Vinit, Winit=Winit)

# ------------- RUN AOA SWEEP --------------------------------------------------
uns.gt.create_path(save_path, true)

for AOA in AOAs
    println("\n\nRunning AOA = $(AOA)")
    Vinf(X, t) = magVinf * [cosd(AOA), sind(AOA), 0]

    monitor_rotor = uns.generate_monitor_rotors(rotors, J, rho, RPM, nsteps;
                                                t_scale=RPM/60,
                                                t_lbl="Revolutions",
                                                save_path=save_path,
                                                run_name="AOA$(ceil(Int, AOA*100))",
                                                disp_conv=false
                                                )

    println("\tRunning simulation...")

    save_vtks = AOA in output_runs ? save_path*"-AOA$(ceil(Int, AOA*100))" : nothing

    uns.run_simulation(simulation, nsteps;
                        Vinf=Vinf,
                        rho=rho, mu=mu, sound_spd=speedofsound,
                        p_per_step=p_per_step,
                        max_particles=max_particles,
                        sigma_vlm_surf=sigma_rotor_surf,
                        sigma_rotor_surf=sigma_rotor_surf,
                        sigma_vpm_overwrite=sigma_vpm_overwrite,
                        vlm_rlx=vlm_rlx,
                        hubtiploss_correction=hubtiploss_correction,
                        shed_unsteady=shed_unsteady,
                        shed_starting=shed_starting,
                        extra_runtime_function=monitor_rotor,
                        save_path=save_vtks,
                        run_name="singlerotor",
                        v_lvl=1, verbose_nsteps=24
                        )
end

Below are the details of my setup and the error:
Setup:

I have set all airfoil CSV files with α values in ascending order.

Spline parameters are set as spl_k = 3 and spl_s = 0.001.

I have modified the incidence sweep code to include CSV files with updated airfoil characteristics (α, Cl, Cd, Cm) values.

The α values are set from -180 to 180 degrees, with differing step sizes for each section.

Parameters like AOAs and magVinf are adjusted accordingly in the simulation code.

Error Message:
I get the following error when running the simulation in Julia:

ERROR: Error on entry, no approximation returned. The following conditions
must hold:
1<=k<=5
x[1] < x[2] < ... < x[end]
w[i] > 0.0 for all i

Additionally, if spline knots are given:
length(xknots) <= length(x) + k + 1
x[1] < xknots[1] < xknots[k+2] < ... < xknots[end] < x[end]
The schoenberg-whitney conditions: there must be a subset of data points
xx[j] such that t[j] < xx[j] < t[j+k+1] for j=1,2,...,n-k-1
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] Dierckx.Spline1D(...)
    @ Dierckx ~/.julia/packages/Dierckx/TDOyl/src/Dierckx.jl:175
  [3] occb_af_from_data(...)
    @ FLOWVLM ~/.julia/packages/FLOWVLM/RWw6g/src/FLOWVLM_rotor_ccb.jl:57
  [4] FLOWVLM2OCCBlade(...)

It seems that there is an issue with the spline knot condition or ordering of data points. I have double-checked the α values, and they are sorted in ascending order from -180 to 180 degrees, with positive weights for all. Despite this, the error persists.

How can I solve this issue?
Thank you very much for your assistance. Any advice or insights would be greatly appreciated!

[Zollhafen]

Hello,

i ran into the same issue when I was trying to run the simulation with mainly negative alpha values in my polar data.
I couldnt find a direct solve for the problem but you might want to experiment or at least look into passing the turbine_flag=true argument in the generate_rotor function.
This will reverse the x_inflow of the ccblade solver and cl if i remember correctly.
This inturn might bring other problems and screw with your sign convention but it might be worth a try.

I solved my problem in manually reversing the sign of alpha in the polars and reversing theta_eff.
Though your polars have both positive and negative alpha values, this migth not work.
I hope my input was useful.

Best regards