Propeller wing interaction (continued)

[EdoAlvarezR]

(starting a new thread since the previous one has become quite long)

@Turbulence1116 @inse0918 @Dhruvin108, I have now added a tutorial reproducing the PROWIM validation case capturing prop-on-wing interactions. I have tested two different levels of fidelity:

I have also added settings for mid-low fidelity, but I'm not confident that mid-low would produce accurate results (it's very fast though). Could you guys try both the mid-low and mid-high fidelity? If the mid-low fidelity is not matching the experiment, perhaps p_per_step=5 would help.

• Mid-low fidelity run time: 13 minutes a Dell Precision 7760 laptop.
• Mid-high fidelity run time: 70 minutes a Dell Precision 7760 laptop.
• High fidelity runtime: ~2 days on a 16-core AMD EPYC 7302 processor.

TUTORIAL: https://flow.byu.edu/FLOWUnsteady/examples/prowim-aero/

[inse0918]

Hi, Edo. Thank you for sharing your 'prowim' script !

I ran this script by changing the fidelity using 'Low' and 'Mid' options.

This is the results about sectional normal force coefficient along the semi spanwise.

The, results obtained by computing the 'Low' fidelity option are not fit with experimental data by Veldhuis.

But, when I using the 'Mid' fidelity option, the results were in good agreement with experimental data.

System Performance

• CPU : Intel Xeon Gold 5222, 3.8GHz, 8 cores
• RAM : 256 GB

Low Fidelity Analysis : 24 min.
Mid Fidelity Analysis : 157 min.

All the analysis were ran individually.

Thank you.

[EdoAlvarezR]

From your runs it seems that either p_per_step=4+SFS or p_per_step=5 without SFS are the coarsest we can go while still getting accurate results. How fast are both cases?

Also vpm_integration = vpm.euler seems to be accurate enough for the temporal scheme (as opposed to using RK3).

[inse0918]

I agree about it.

[p_per_step=4+SFS] was almost required 2 hours and [p_per_step=5 no SFS] was required under 1 hours.

So, the single propeller with wing validation cases are sufficiently validated.

Next step is multi propeller with wing case.

Reference paper : "Numerical Investigation of Wing–Multiple Propeller Aerodynamic Interaction Using Actuator Disk Method"

[EdoAlvarezR]

Great! @inse0918, is P5CONC a RANS simulation?

[EdoAlvarezR]

@inse0918, any progress comparing to that paper? I'm looking forward to seeing your results!

(something to watch out for: actuator disk models often don't capture any swirl velocity, which can significantly change the results of a blown wing. The actuator line model in the VPM captures the swirl velocity, so it might be more physically accurate and won't necessarily match the actuator disk)

[inse0918]

Hi @EdoAlvarezR , yes these results are ran in RANS simulation.

In physically, the Actuator Line Model could be more accurate than the Actuator Disk Model as you mentioned.

I'll compare it soon(maybe, next week :) ).

[Turbulence1116]

Hi Edo, thank for your scripts, it works well at AoA=4, with a pitch at 2 degrees.

However, I have checked my script, apart from the "wake treatment", everything else are the same as yours, my sectional lift is over-predicted quite a lot. Strangely, at AoA=0, it works well. I am thinking that this over-prediction may be related to the wake? If I understand correctly, the wake treatment removes less relevant particles from the wake (the particles shedding from the middle of the wing for example), thus reducing the computational time (40mins to 30mins).

The Euler time integration without SFS seems to work well enough in this case, under my workstation, i7-13700 32G RAM, it takes 30mins for 10 revolution, for Mid fidelity.

#=#############################################################################
# DESCRIPTION
    Simulation of prowim, Veldhuis 2005
    Strong interaction of propeller wake and wing wake, higher level of discretisation should be applied.
    Convergence for reference, https://flow.byu.edu/FLOWUnsteady/theory/convergence/#Immersed-Vorticity
    180 sheds per revo, 50 wing elements and 20 rotor elements should be sufficient to have converged results.
    Full wing span 1.28m is simulated to remove the reflection plate in the experiment.
    
    The current setup with reference to: https://flow.byu.edu/FLOWUnsteady/examples/prowim-aero/   
    
    Euler, no SFS, 50 elements on wing, 20 elements on blade, 36 time steps per revo (10 degress each step), 
    5 sheds per step, total 180 sheds per revo
    10 revos, 40mins
    
    Total 0.0407s, 360 steps and dt = 0.000113, hub correction on, denotes "" 
    
    Freesteam: 49.5m/s
    PropRadi: 0.118m
    J: 0.85
    RPM: 14743 
    TipMach: 0.534
    ReD: 2.86e6
    ReC: 784985 (<8e5)
=###############################################################################


import FLOWUnsteady as uns
import FLOWUnsteady: vlm, vpm, gt, Im
import FLOWVLM as vlm
import FLOWVPM as vpm

run_name        = "prowim4"                  # Name of this simulation
save_path       = "prowim4"                  # Where to save this simulation
paraview        = false                           # Whether to visualize with Paraview

# ----------------- GEOMETRY PARAMETERS ----------------------------------------
# Wing geometry
b               = 1.28                     # (m) span length
ar              = 5.33                       # Aspect ratio b/c_tip
tr              = 1.0                       # Taper ratio c_tip/c_root
twist_root      = 0.0                       # (deg) twist at root
twist_tip       = 0.0                       # (deg) twist at tip
lambda          = 0.0                      # (deg) sweep
gamma           = 0.0                       # (deg) dihedral

thickness_w     = 0.15                      # Wing airfoil thickness t/c (ASM activation)

# Rotor geometry
rotor_file      = "beaver.csv"              # Rotor geometry
data_path       = uns.def_data_path         # Path to rotor database
read_polar      = vlm.ap.read_polar2        # What polar reader to use
pitch           = 2.0                       # (deg) collective pitch of blades, verified that C_t is around 0.12 at N_shed = 360
CW              = false                     # Clock-wise rotation
xfoil           = false                      # Whether to run XFOIL
#ncrit           = 9                         # Turbulence criterion for XFOIL

# Wing Discretization
n_wing          = 50                        # Number of spanwise elements per side
r_wing          = 2.0                       # Geometric expansion of elements
# Rotor Discretization
n_rotor         = 20                        # Number of blade elements per blade
r_rotor         = 1/10                       # Geometric expansion of elements

# Vehicle assembly
AOAwing         = 0.0                        # (deg) wing angle of attack
spanpos         = [0.47,-0.47]               # Semi-span position of each rotor, 2*y/b
xpos            = [-0.84,-0.84]              # x-position of rotors relative to LE, x/c
zpos            = [0.0,0.0]                  # z-position of rotors relative to LE, z/c
CWs             = [false,true]               # Clockwise rotation for each rotor
nrotors         = length(spanpos)            # Number of rotors

# Check that we declared all the inputs that we need for each rotor
@assert length(spanpos)==length(xpos)==length(zpos)==length(CWs) ""*
    "Invalid rotor inputs! Check that spanpos, xpos, zpos, and CWs have the same length"

# ----------------- SIMULATION PARAMETERS --------------------------------------
# Vehicle motion
magVvehicle     = 0.0                      # (m/s) vehicle velocity
AOA             = 4.0                       # (deg) freestream angle of attack

# Freestream
#magVinf         = 1e-8                      # (m/s) freestream velocity
magVinf         = 49.5 
rho             = 1.225                       # (kg/m^3) air density
mu              = 1.79e-5
speedofsound    = 342.35                    # (m/s) speed of sound
magVref         = sqrt(magVinf^2 + magVvehicle^2) # (m/s) reference velocity
qinf            = 0.5*rho*magVref^2         # (Pa) reference static pressure

#Vinf(X, t)      = t==0 ? magVvehicle*[1,0,0] : magVinf*[1,0,0] # Freestream function
Vinf(X, t)      = magVinf*[cosd(AOA), 0.0, sind(AOA)]  # Freestream function

# Rotor operation
# Read radius of this rotor and number of blades
R, B            = uns.read_rotor(rotor_file; data_path=data_path)[[1,3]]
J               = 0.85                       # Advance ratio Vref/(nD)
RPM             = 60*magVref/(J*2*R)        # RPM

Rec             = rho * magVref * (b/ar) / mu       # Chord-based wing Reynolds number
ReD             = 2*pi*RPM/60*R * rho/mu * 2*R      # Diameter-based rotor Reynolds number
Matip           = 2*pi*RPM/60 * R / speedofsound    # Tip Mach number

println("""
    Vref:   $(round(magVref, digits=1)) m/s
    RPM:    $(RPM)
    Matip:  $(round(Matip, digits=3))
    ReD:    $(round(ReD, digits=0))
    Rec:    $(round(Rec, digits=0))
""")

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

# Aerodynamic solver
VehicleType     = uns.UVLMVehicle           # Unsteady solver
# VehicleType     = uns.QVLMVehicle         # Quasi-steady solver
const_solution  = VehicleType==uns.QVLMVehicle  # Whether to assume that the
                                                # solution is constant or not
# Time parameters
nrevs           = 10                        # Number of revolutions in simulation
nsteps_per_rev  = 36                        # Time steps per revolution, at 180 simulation dose not advance
                                            # Total time step is ?
nsteps          = const_solution ? 2 : nrevs*nsteps_per_rev # Number of time steps
#nsteps          = nrevs*nsteps_per_rev # Number of time steps, solution is assumed to b
ttot            = nsteps/nsteps_per_rev / (RPM/60)       # (s) total simulation time

# VPM particle shedding
p_per_step      = 5                         # Sheds per time step
shed_starting   = false                      # Whether to shed starting vortex, (NOTE: starting vortex might make simulation unstable with AOA>8)
shed_unsteady   = false                      # Whether to shed vorticity from unsteady loading
unsteady_shedcrit = 0.001                   # Shed unsteady loading whenever circulation
                                            #  fluctuates by more than this ratio
max_particles   = 10*nrotors*((2*n_rotor+1)*B)*nsteps*p_per_step + 1 # Maximum number of particles, can I turn it off or define manually?
max_particles   += (nsteps+1)*(2*n_wing*(p_per_step+1) + p_per_step)
#max_particles   = 10000000

# ---------- Vortex sheet parameters (ASM activation)---------------
vlm_vortexsheet = true                          # Whether to spread the wing circulation as a vortex sheet
vlm_vortexsheet_overlap = 2.125/10                 # Overlap of the particles that make the vortex sheet
vlm_vortexsheet_distribution = uns.g_pressure   # Distribution of the vortex sheet

# Regularization of embedded vorticity
# Note strong wake interactions (e.g., a rotor wake impinging on a wing surface), the size of the particles used to immerse the surface vorticity
# can change the effects of the interaction.
# the size of the embedded particles needs to be about two orders of magnitude smaller than the relevant length scale.
# the recommended sigma is 0.01 of length scale 
sigma_vlm_surf  = b/200                     # VLM-on-VPM smoothing radius of the lifting bound vorticity of wings 
sigma_rotor_surf= R/80                      # Rotor-on-VPM smoothing radius of the rotor bound vorticity
vlm_vortexsheet_sigma_tbv = thickness_w*(b/ar) / 100 # Smoothing radius of trailing bound vorticity of wing (ASM activation)
vlm_vortexsheet_maxstaticparticle = 1000000 # particles to preallocate for the vortex sheet (ASM activation) 

lambda_vpm      = 2.125                     # VPM core overlap

# VPM smoothing radius
sigma_vpm_overwrite = lambda_vpm * 2*pi*R/(nsteps_per_rev*p_per_step) # Attention, sigma should be 0.01 of R, that is 2.125*2*pi/(72*20)=0.00925

sigmafactor_vpmonvlm= 1                     # Shrink particles by this factor when
                                            #  calculating VPM-on-VLM/Rotor induced velocities
sigma_vlm_solver= -1                        # VLM-on-VLM smoothing radius (deactivated with <0) (ASM activation)
                                        
# Rotor solver
vlm_rlx         = 0.3                       # VLM relaxation <-- this also applied to rotors
#hubtiploss_correction = vlm.hubtiploss_nocorrection # Hub and tip correction
hubtiploss_correction = ( (0.75, 10, 0.5, 0.05), (1, 1, 1, 1.0) )
# VPM solver
#vpm_integration = vpm.rungekutta3           # VPM temporal integration scheme
vpm_integration = vpm.euler

vpm_viscous     = vpm.Inviscid()            # VPM viscous diffusion scheme
# vpm_viscous   = vpm.CoreSpreading(-1, -1, vpm.zeta_fmm; beta=100.0, itmax=20, tol=1e-1)

vpm_SFS         = vpm.SFS_none              # VPM LES subfilter-scale model
#vpm_SFS       = vpm.SFS_Cd_threelevel_nobackscatter #Should be turned on since high turbulence is expected in the wake 
# vpm_SFS       = vpm.DynamicSFS(vpm.Estr_fmm, vpm.pseudo3level_positive;
#                                   alpha=0.999, maxC=1.0,
#                                   clippings=[vpm.clipping_backscatter])

if VehicleType == uns.QVLMVehicle
    # Mute warnings regarding potential colinear vortex filaments. This is
    # needed since the quasi-steady solver will probe induced velocities at the
    # lifting line of the blade
    uns.vlm.VLMSolver._mute_warning(true)
end

# ----------------- WAKE TREATMENT ---------------------------------------------
# NOTE: It is known in the CFD community that rotor simulations with an
#       impulsive RPM start (*i.e.*, 0 to RPM in the first time step, as opposed
#       to gradually ramping up the RPM) leads to the hub "fountain effect",
#       with the root wake reversing the flow near the hub.
#       The fountain eventually goes away as the wake develops, but this happens
#       very slowly, which delays the convergence of the simulation to a steady
#       state. To accelerate convergence, here we define a wake treatment
#       procedure that suppresses the hub wake for the first three revolutions,
#       avoiding the fountain effect altogether.
#       This is especially helpful in low and mid-fidelity simulations.

suppress_fountain   = true                 # Toggle

# Supress wake shedding on blade elements inboard of this r/R radial station
no_shedding_Rthreshold = suppress_fountain ? 0.35 : 0.0

# Supress wake shedding for this many time steps
no_shedding_nstepsthreshold = 3*nsteps_per_rev

omit_shedding = []          # Index of blade elements to supress wake shedding

# Function to suppress or activate wake shedding
function wake_treatment_supress(sim, args...; optargs...)

    # Case: start of simulation -> suppress shedding
    if sim.nt == 1

        # Identify blade elements on which to suppress shedding
        for i in 1:vlm.get_m(rotor)
            HS = vlm.getHorseshoe(rotor, i)
            CP = HS[5]

            if uns.vlm.norm(CP - vlm._get_O(rotor)) <= no_shedding_Rthreshold*R
                push!(omit_shedding, i)
            end
        end
    end

    # Case: sufficient time steps -> enable shedding
    if sim.nt == no_shedding_nstepsthreshold

        # Flag to stop suppressing
        omit_shedding .= -1

    end

    return false
end

println("""
    Resolving wake for $(round(ttot*magVref/b, digits=1)) span distances
    Revolution: $(nrevs)
    Total simulation time :   $(ttot) s for $(nsteps) steps
""")

# ----------------- 1) VEHICLE DEFINITION --------------------------------------

# -------- Generate components
println("Generating geometry...")

# Generate wing
wing = vlm.simpleWing(b, ar, tr, twist_root, lambda, gamma;
                                    twist_tip=twist_tip, n=n_wing, r=r_wing);

# Pitch wing to its angle of attack
O = [0.0, 0.0, 0.0]                                 # New position
Oaxis = uns.gt.rotation_matrix2(0, -AOAwing, 0)     # New orientation
vlm.setcoordsystem(wing, O, Oaxis)

# Generate rotors
rotors = vlm.Rotor[]
for ri in 1:nrotors #nrotors is the number of rotors

    # Generate rotor
    rotor = uns.generate_rotor(rotor_file;
                                        pitch=pitch,
                                        n=n_rotor, CW=CWs[ri], blade_r=r_rotor,
                                        altReD=[RPM, J, mu/rho],
                                        xfoil=xfoil,
                                        #ncrit=ncrit,#it should be turned off when xfoil is false
                                        read_polar=read_polar, #it should be here when xfoil is false
                                        data_path=data_path,
                                        verbose=true,
                                        verbose_xfoil=false,
                                        plot_disc=true
                                        );

    # Determine position along wing LE
    y = spanpos[ri]*b/2
    x = abs(y)*tand(lambda) + xpos[ri]*b/ar
    z = abs(y)*tand(gamma) + zpos[ri]*b/ar

    # Account for angle of attack of wing
    nrm = sqrt(x^2 + z^2)
    snx = sign(x)
    #x = nrm*cosd(AOAwing)
    x = snx * nrm * cosd(AOAwing)
    z = -nrm*sind(AOAwing)

    # Translate rotor to its position along wing
    O = [x, y, z]                                       # New position
    Oaxis = uns.gt.rotation_matrix2(0, 0, 180)          # New orientation
    vlm.setcoordsystem(rotor, O, Oaxis)

    push!(rotors, rotor)
end


# -------- Generate vehicle
println("Generating vehicle...")

# System of all FLOWVLM objects
system = vlm.WingSystem()

vlm.addwing(system, "Wing", wing)

for (ri, rotor) in enumerate(rotors)
    vlm.addwing(system, "Rotor$(ri)", rotor)
end

# System solved through VLM solver
vlm_system = vlm.WingSystem()
vlm.addwing(vlm_system, "Wing", wing)

# Systems of rotors
rotor_systems = (rotors, );

# System that will shed a VPM wake
wake_system = vlm.WingSystem()              # System that will shed a VPM wake
vlm.addwing(wake_system, "Wing", wing)
                                            # NOTE: Do NOT include rotor when using the quasi-steady solver
if VehicleType != uns.QVLMVehicle
    for (ri, rotor) in enumerate(rotors)
        vlm.addwing(wake_system, "Rotor$(ri)", rotor)
    end
end

# Pitch vehicle to its angle of attack
O = [0.0, 0.0, 0.0]                                 # New position
Oaxis = uns.gt.rotation_matrix2(0, -AOA, 0)         # New orientation
vlm.setcoordsystem(system, O, Oaxis)

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


# ------------- 2) MANEUVER DEFINITION -----------------------------------------
# Non-dimensional translational velocity of vehicle over time
Vvehicle(t) = [-1, 0, 0]        # <---- Vehicle is traveling in the -x direction

# Angle of the vehicle over time
anglevehicle(t) = zeros(3)

# RPM control input over time (RPM over `RPMref`)
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 = magVvehicle                          # 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)  # Initial angular velocity

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

# ---------- Force calculation parameters (ASM activation)----------
KJforce_type                = "regular"     # KJ force evaluated at middle of bound vortices
# KJforce_type              = "averaged"  # KJ force evaluated at average vortex sheet
# KJforce_type              = "weighted"  # KJ force evaluated at strength-weighted vortex sheet

include_trailingboundvortex = false        # Include trailing bound vortices in force calculations

include_freevortices        = false         # Include free vortices in force calculation
include_freevortices_TBVs   = false         # Include trailing bound vortex in free-vortex force

include_unsteadyforce       = true          # Include unsteady force
add_unsteadyforce           = false         # Whether to add the unsteady force to Ftot or to simply output it

include_parasiticdrag       = true          # Include parasitic-drag force
add_skinfriction            = true          # If false, the parasitic drag is purely parasitic, meaning no skin friction
calc_cd_from_cl             = false         # Whether to calculate cd from cl or effective AOA
wing_polar_file             = "xf-n64015-il-1000000.csv"    # Airfoil polar for parasitic drag


# ---------- Aerodynamic forces (ASM activation) --------------

forces = []

# Calculate Kutta-Joukowski force
kuttajoukowski = uns.generate_aerodynamicforce_kuttajoukowski(KJforce_type,
                                sigma_vlm_surf, sigma_rotor_surf,
                                vlm_vortexsheet, vlm_vortexsheet_overlap,
                                vlm_vortexsheet_distribution,
                                vlm_vortexsheet_sigma_tbv;
                                vehicle=vehicle)
push!(forces, kuttajoukowski)

# Free-vortex force
if include_freevortices
    freevortices = uns.generate_calc_aerodynamicforce_freevortices(
                                            vlm_vortexsheet_maxstaticparticle,
                                            sigma_vlm_surf,
                                            vlm_vortexsheet,
                                            vlm_vortexsheet_overlap,
                                            vlm_vortexsheet_distribution,
                                            vlm_vortexsheet_sigma_tbv;
                                            Ffv=uns.Ffv_direct,
                                            include_TBVs=include_freevortices_TBVs
                                            )
    push!(forces, freevortices)
end

# Force due to unsteady circulation
if include_unsteadyforce
    unsteady(args...; optargs...) = uns.calc_aerodynamicforce_unsteady(args...; add_to_Ftot=add_unsteadyforce, optargs...)

    push!(forces, unsteady)
end

# Parasatic-drag force (form drag and skin friction)
if include_parasiticdrag
    parasiticdrag = uns.generate_aerodynamicforce_parasiticdrag(wing_polar_file;
                                                            read_path=joinpath(data_path, "airfoils"),
                                                            calc_cd_from_cl=calc_cd_from_cl,
                                                            add_skinfriction=add_skinfriction,
                                                            Mach=speedofsound!=nothing ? magVref/speedofsound : nothing)
    push!(forces, parasiticdrag)
end


# Stitch all the forces into one function
function calc_aerodynamicforce_fun(vlm_system, args...; per_unit_span=false, optargs...)

    # Delete any previous force field
    fieldname = per_unit_span ? "ftot" : "Ftot"
    if fieldname in keys(vlm_system.sol)
        pop!(vlm_system.sol, fieldname)
    end

    Ftot = nothing

    for (fi, force) in enumerate(forces)
        Ftot = force(vlm_system, args...; per_unit_span=per_unit_span, optargs...)
    end

    return Ftot
end


# ------------- 4) MONITORS DEFINITIONS ----------------------------------------

# Generate function that computes wing aerodynamic forces
#calc_aerodynamicforce_fun = uns.generate_calc_aerodynamicforce(;
#                                    add_parasiticdrag=true,
#                                    add_skinfriction=true,
#                                    airfoilpolar="xf-n64015-il-1000000.csv"
#                                    )

#L_dir(t) = [0, 0, 1]                        # Direction of lift
#D_dir(t) = [1, 0, 0]                        # Direction of drag

D_dir = [cosd(AOA), 0.0, sind(AOA)]         # Direction of drag
#L_dir = uns.cross(D_dir, [0,1,0])           # Direction of lift
L_dir = [-sind(AOA), 0, cosd(AOA)]          # Direction of lift

figs, figaxs = [], []                       # Figures generated by monitor

# Generate wing monitor
monitor_wing = uns.generate_monitor_wing(wing, Vinf, b, ar,
                                            rho, qinf, nsteps;
                                            include_trailingboundvortex=include_trailingboundvortex, #ASM activation
                                            calc_aerodynamicforce_fun=calc_aerodynamicforce_fun,
                                            L_dir=L_dir,
                                            D_dir=D_dir,
                                            save_path=save_path,
                                            run_name=run_name*"-wing",
                                            figname="wing monitor",
                                            )

# Generate rotors monitor
monitor_rotors = uns.generate_monitor_rotors(rotors, J, rho, RPM, nsteps;
                                            t_scale=RPM/60,        # Scaling factor for time in plots
                                            t_lbl="Revolutions",   # Label for time axis
                                            save_path=save_path,
                                            run_name=run_name*"-rotors",
                                            figname="rotors monitor",
                                            )
# Concatenate monitors
monitors = uns.concatenate(monitor_wing, monitor_rotors)


# ------------- 5) RUN SIMULATION ----------------------------------------------
println("Running simulation...")

#If Matip is different than zero while running XFOIL, you must deactive compressibility corrections in run_simulation by using sound_spd=nothing (sound_spd=speedofsound). Otherwise, compressibility effects will be double accounted for.

# Run simulation
uns.run_simulation(simulation, nsteps;
                    # ----- SIMULATION OPTIONS -------------
                    Vinf=Vinf,
                    rho=rho, mu=mu, sound_spd=nothing,
                    # ----- SOLVERS OPTIONS ----------------
                    p_per_step=p_per_step,
                    max_particles=max_particles,
                    vpm_integration=vpm_integration,
                    vpm_viscous=vpm_viscous,
                    vpm_SFS=vpm_SFS,
                    vlm_vortexsheet=vlm_vortexsheet,#ASM activation
                    vlm_vortexsheet_overlap=vlm_vortexsheet_overlap,#ASM activation
                    vlm_vortexsheet_distribution=vlm_vortexsheet_distribution,#ASM activation
                    vlm_vortexsheet_sigma_tbv=vlm_vortexsheet_sigma_tbv,#ASM activation
                    max_static_particles=vlm_vortexsheet_maxstaticparticle,#ASM activation
                    sigma_vlm_solver=sigma_vlm_solver,
                    sigma_vlm_surf=sigma_vlm_surf,
                    sigma_rotor_surf=sigma_rotor_surf,
                    sigma_vpm_overwrite=sigma_vpm_overwrite,
                    sigmafactor_vpmonvlm=sigmafactor_vpmonvlm,
                    vlm_rlx=vlm_rlx,
                    hubtiploss_correction=hubtiploss_correction,
                    shed_starting=shed_starting,
                    shed_unsteady=shed_unsteady,
                    unsteady_shedcrit=unsteady_shedcrit,
                    extra_runtime_function=monitors,
                    # ----- OUTPUT OPTIONS ------------------
                    save_path=save_path,
                    run_name=run_name,
                    save_wopwopin=false,  # <--- Generates input files for PSU-WOPWOP noise analysis
                    save_code=@__FILE__
                    );




# ----------------- 6) Extract sectional Cl and Cd -------------------------------------------
import FLOWUnsteady: cross, dot, norm, plt, @L_str

println("Saving sectional data...")

Xac         = [0.25*b/ar, 0, 0]             # (m) aerodynamic center for moment calculation
ls, ds      = [], []                        # Load and drag distributions
spanposs    = []                            # Spanwise positions for load distributions
# Integrate total lift and drag
L = sum(wing.sol["L"])
D = sum(wing.sol["D"])

# Lift and drag coefficients
CL = norm(L) / (qinf*b^2/ar)
CD = norm(D) / (qinf*b^2/ar)

# Control point of each element
Xs = [vlm.getControlPoint(wing, i) for i in 1:vlm.get_m(wing)]
#Xs is an array with the location of the control points in space in the global coordinate system

# # Force of each element
Fs = wing.sol["Ftot"]

# Integrate the total moment with respect to aerodynamic center
M = sum( cross(X - Xac, F) for (X, F) in zip(Xs, Fs) )

# Integrated moment decomposed into rolling, pitching, and yawing moments
Shat = [0, 1, 0]                            # Spanwise direction
Dhat = [cosd(AOA), 0.0, sind(AOA)]          # Direction of drag
Lhat = cross(Dhat, Shat)                    # Direction of lift

lhat = Dhat                   # Rolling direction
mhat = Shat                   # Pitching direction
nhat = Lhat                   # Yawing direction

roll = dot(M, lhat)
pitch = dot(M, mhat)
yaw = dot(M, nhat)

# Sectional loading (in vector form) at each control point
fs = wing.sol["ftot"]

# Decompose vectors into lift and drag distribution
l = [ dot(f, Lhat) for f in fs ]
d = [ dot(f, Dhat) for f in fs ]

# Span position of each control point
spanpos = [ dot(X, Shat) / (b/2) for X in Xs ]

nondim = 0.5*rho*magVinf^2*b/ar   # Normalization factor

cls = l / nondim
cds = d / nondim

using DelimitedFiles

# saving the section aerodynamic data
writedlm(save_path*"_cl_cd.dat", zip(cls, cds, spanpos))

[Turbulence1116]

No worries I'll figure it out!

[EdoAlvarezR]

On a second glance, I think I found what the issue is. In these lines you are pitching the vehicle by the angle of attack

# Pitch vehicle to its angle of attack
O = [0.0, 0.0, 0.0]                                 # New position
Oaxis = uns.gt.rotation_matrix2(0, -AOA, 0)         # New orientation
vlm.setcoordsystem(system, O, Oaxis)

while in this line you also have defined the freestream to come at an angle of attack

Vinf(X, t)      = magVinf*[cosd(AOA), 0.0, sind(AOA)]  # Freestream function

so you are effectively running the simulation at twice the AOA that you were intending to. I recommend not pitching the vehicle and just let the freestream come at an AOA.

[Turbulence1116]

Thank you for spotting that, it is exactly where my mistake was!

[Turbulence1116]

And the normal force coef = lift coef / cos (AoA) ?

[inse0918]

Hi jon.

Muhammad Nazmul Haque et al.: Experimental Investigation on the Performance of NACA 4412 Aerofoil with Curved Leading Edge Planform

[inse0918]

Multi Propeller-Wing Interaction Analysis

Analysis Condition

• Vinf : 30.13 m/s
• AoA : 4 deg.
• RPM : 1,310
• Diameter : 1.78 m
• J : 0.7753
• Propeller : Beaver with 2 blades
• Wing Airfoil : NACA0012 ("naca0012polar.csv" was used)
• b = 5.7*2 m, c = 1.295 m

# Vehicle assembly
AOAwing         = 0.0                       # (deg) wing angle of attack
spanpos         = [-1.000, -0.6146, -0.2293, 0.2293, 0.6146, 1.000]       # Semi-span position of each rotor, 2*y/b
xpos            = [-0.2317, -0.2317, -0.2317, -0.2317,-0.2317, -0.2317]        # x-position of rotors relative to LE, x/c
zpos            = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]                # z-position of rotors relative to LE, z/c
CWs             = [true, true, true, false, false, false]           # Rotation direction of each rotor: inboard up

Isolated Wing Cases

• isowing1 : ran on prowim.jl based script / p_per_step=5
• isowing2 : ran on wing.jl based script
• vspaero : ran on vspaero of the OpenVSP

The sectional lift coefficient along the half-span is not fit with CFD-RANS results when the propeller is not equipped(isowing)

Also, because of the propeller diameter size is bigger then a PROWIM, the airfoil DB cannot used in this cases(xfoil=false).

Besides, when the xfoil option turned on, the calculation is not going well likewise xfoil=false case with this error.

ERROR: LoadError: The maximal number of iterations maxit (set to 20 by the program) allowed for finding a smoothing spline with fp=s has been reached: s too small. There is an approximation returned but the corresponding weighted sum of squared residuals does not satisfy the condition abs(fp-s)/s < tol.

Should I make the beaver propeller airfoil DB for using this cases ?

Thank you.
Best regards,
Inseo.

[EdoAlvarezR]

The simulation should be spot on for an isolated with a symmetric airfoil (it's simply a VLM with an unsteady wake). There has to be something off in the simulation setup. How far are you developing the wake for? Typically you need to resolve the wake 3 or 4 wing-span distances downstream of the trailing edge for the loading to converge.

[EdoAlvarezR]

[inse0918]

Hi Edo, sorry for late reply

This is a results about the reolved spans was set as 2.8.

[inse0918]

Hi Edo, these are the results about the P3COEQ.

VLM-ADT(VSPAERO)
VLM-LLM(FLOWUnsteady)

In my cases, the DEP case was predicted well using FLOWUnsteady.

Thank you.
Best regards,
Inseo.

[EdoAlvarezR]

Thanks for sharing your results. I look forward to your paper!

[jayjiejie]

I tired to generate an INSEAN E779a propeller and ended up with the same error when generating rotor geometry.

ERROR: LoadError: The maximal number of iterations maxit (set to 20 by the program) allowed for finding a smoothing spline with fp=s has been reached: stoo small. There is an approximation returned but the corresponding weighted sum of squared residuals does not satisfy the condition abs(fp-s)/s < tol.

It has something to do with Dierckx.Spline1D, any suggestions?

Thank you very much.

[EdoAlvarezR]

I would have to see the full error stracktrace to know what spline is failing, but my guess is that is the spline used in the airfoil polar. I'd recommend using a different spline_s value: LINK

@inse0918

[jayjiejie]

Hi Edo

Thank you very much for your advice. I am actually using the propeller.jl from the examples/propeller case as the template. After adding in the extra Discretization Setting parameters, I am able to generate the propeller blade geometry. However, the geometry is not quite the same as the CAD file. I will try to go over again with my chorddist, twistdist, sweepdist, heightdist files.

Kind regards,
Jay

[jayjiejie]

Hi Edo

I tried to play around with the geometrical files, e.g. chorddist.csv, weistdist.csv, sweepdist.csv, heightdist.csv and airfoils.csv and my blade profile in FLOWUnsteady is still not the same as the CAD file. I am getting a bit confused with the followings:

1. twistdist vs. pitchdist
   According to the web page for Generating components - Rotor https://flow.byu.edu/FLOWUnsteady/api/flowunsteady-vehicle-components/#Rotor, it uses twistdist.csv. However, the /examples/propeller tutorial uses apc10x7_pitchdist.csv (located in /database/rotors directory) and its content still use the column name of twist inside. Does it matter if twistdist or pitchdist is used in the blade.csv file? My understanding (which could be wrong) are the twist of the propeller blade is the angle between the profile chord of the blade tip and the profile chord of the blade root and the pitch is the angle between the chord profile and the plane of rotation. Hence, the two definitions are different.

2. sweepdist
   For sweepdist, the definition (found in the FLOWUnsteady) is y/R and y distance of LE from the middle point of hub. Would you be able to clarify the definition for the middle point of hub?

3. heightdist
   This is a similar question as the item 2 above. The heightdist is defined (in the FLOWUnsteady) as z/R height of leading edge from top face of hub. How do we quantify this parameter from the CAD file? Which part of the leading edge we should take as the measurement point? Where is the top surface of hub?

4. airfoils.csv
   For all the aero files in the airfoil.csv, presumably they are the airfoil profiles normalised with the chrod length (c becomes 1), LE at c = 0 and TE at c = 1, AoA = 0 deg. Does it matter if the points start from LE (TE), or from the top to bottom surface (bottom to top surface)?

I apologise for many questions but I feel it is worth spending time to get those definitions clarified so that we are not making silly mistakes when generating the geometrical files from the CAD.

Thank you very much.

Jay