Propeller wing interaction #98
-
Dear all, I am setting up a simulation to replicate Veldhuis's Prowim experiment, now I got some questions regarding to my results. First, the results of total lift coef vs AoA on the clean wing are reasonable, despite slightly over-predicted values when AoA is greater than 4 degrees. When the propeller is added to the wing, the predicted lift is way too big as the AoA increases. Second, the sectional lift coef distribution is plotted, the trend looks good, however, it is obvious that the lift coef cannot be that big. Regarding to the results, I think there is something I must have missed that is related to the propeller, since the results are good on the clean wing. I checked that the compressibility is turned off by setting speed of sound = nothing and there is no difference even the compressibility is taken into account. I then checked the revelant parameters for propeller, they are inline with Veldhuis's Prowim experiment. So far, I haven't found anything that I misunderstood related to the calculation of the Cl. Next, I think I am going to check the propeller's performance. I also attach the code here, maybe I did make some mistakes and I did not find them. Your help is much appreciated! #=#############################################################################
# DESCRIPTION
Simulation of prowim, Veldhuis 2005
Full wing span 1.28m
Euler time, no SFS, ALM, 40 elements on wing, 20 elements on blade: 20 revolutions (3.1-span wake) took 70 mins and 500 000 particles
dt: 0.000113
Freesteam: 49.5m/s
PropRadi: 0.118m
J: 0.85
RPM: 14743 is is too fast?
TipMach: 0.534
ReD: 2.86e6
ReC: 784985 (<8e5)
=###############################################################################
import FLOWUnsteady as uns
import FLOWVLM as vlm
import FLOWVPM as vpm
run_name = "prowim_2" # Name of this simulation
save_path = "prowim_2" # 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
# Wing Discretization
n_wing = 40 # Number of spanwise elements per side
r_wing = 1.0 # Geometric expansion of elements
# Rotor geometry
rotor_file = "beaver.csv" # Rotor geometry
data_path = uns.def_data_path # Path to rotor database
pitch = 0.0 # (deg) collective pitch of blades?
CW = false # Clock-wise rotation
xfoil = true # Whether to run XFOIL
ncrit = 9 # Turbulence criterion for XFOIL
# Rotor Discretization
n_rotor = 20 # Number of blade elements per blade
r_rotor = 1/5 # 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 = 2.0 # (deg) vehicle angle of attack
# Freestream
#magVinf = 1e-8 # (m/s) freestream velocity
magVinf = 49.5
rho = 1.225 # (kg/m^3) air density
mu = 1.85508e-5 # (kg/ms) air dynamic viscosity
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 = 20 # Number of revolutions in simulation
nsteps_per_rev = 36 # Time steps per revolution
# Total time step is 720
nsteps = const_solution ? 2 : nrevs*nsteps_per_rev # Number of time steps
ttot = nsteps/nsteps_per_rev / (RPM/60) # (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
unsteady_shedcrit = 0.001 # Shed unsteady loading whenever circulation
# fluctuates by more than this ratio
max_particles = nrotors*((2*n_rotor+1)*B)*nsteps*p_per_step + 1 # Maximum number of particles
max_particles += (nsteps+1)*(2*n_wing*(p_per_step+1) + p_per_step)
# Regularization
sigma_vlm_surf = b/100 # VLM-on-VPM smoothing radius
sigma_rotor_surf= R/50 # Rotor-on-VPM smoothing radius
lambda_vpm = 2.125 # VPM core overlap
# VPM smoothing radius
sigma_vpm_overwrite = lambda_vpm * 2*pi*R/(nsteps_per_rev*p_per_step)
sigmafactor_vpmonvlm= 1 # Shrink particles by this factor when
# calculating VPM-on-VLM/Rotor induced velocities
# Rotor solver
vlm_rlx = 0.5 # VLM relaxation <-- this also applied to rotors
hubtiploss_correction = vlm.hubtiploss_nocorrection # Hub and tip correction
# 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.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
println("""
Resolving wake for $(round(ttot*magVref/b, digits=1)) span distances
""")
# ----------------- 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,
data_path=data_path,
verbose=true,
verbose_xfoil=false,
plot_disc=false
);
# 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);
# ------------- 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
# Generate wing monitor
monitor_wing = uns.generate_monitor_wing(wing, Vinf, b, ar,
rho, qinf, nsteps;
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,
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=true, # <--- 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)]
# # 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)) |
Beta Was this translation helpful? Give feedback.
Replies: 10 comments 45 replies
-
Related to this other discussion: LINK
Let me know if you have questions as you set up the actuator surface model. |
Beta Was this translation helpful? Give feedback.
-
I have seen the unphysical wake behavior that you are encountering when using too coarse spatial discretization in cases where there should a lot of turbulence (like in this case that the wing should be breaking down the wake structure), which in turn messes up the aero loadings. As a reference, this is the level of spatial discretization I used in the rotor-wing interaction results presented in my dissertation: I see that you are using
In the docs/theory I mention the rule of |
Beta Was this translation helpful? Give feedback.
-
mmmmh, any chance you could upload the particle field (.h5 and .xmf files) from the last time step and flow field that you generated for me to take a look at the simulation? |
Beta Was this translation helpful? Give feedback.
-
Beta Was this translation helpful? Give feedback.
-
Hi. I write the script about propeller-wing interaction using ALM(not a ASM). In my cases, the increased lift coefficient caused by upwash was estimated well at 4 deg. But, the decreased lift coefficient by downwash cannot caputured well. Most important thing is, what is differ with your ALM code with my script. Could you run this script in your pc or laptop ? for comparing the results ? Thank you. Best regards,
|
Beta Was this translation helpful? Give feedback.
-
@Turbulence1116, any updates? |
Beta Was this translation helpful? Give feedback.
-
Hey Inseo@inse0918, I have tried your ALM script, the curve seems to be similar to my ALM. It was ran for 8 revolution, on my workstation it took 15mins: |
Beta Was this translation helpful? Give feedback.
-
Hi guys! Could one of you provide me with the latest version of your script so I can run the simulation and debug from there? (hopefully I can find time during the holidays to take a look at it) |
Beta Was this translation helpful? Give feedback.
-
Just letting you guys know that I was able to re-run the full PROWIM validation, obtaining exactly the same results shown in my dissertation (relatively good agreement with the experiment at 0deg, 4deg, and 10deg). Each simulation takes one or two hours on an EPYC node. I'll convert that script into a tutorial in the documentation in the next few days. |
Beta Was this translation helpful? Give feedback.
-
Closing this thread and continuing the discussion in #129. |
Beta Was this translation helpful? Give feedback.
Closing this thread and continuing the discussion in #129.