Propeller wing interaction #98
Answered
by
EdoAlvarezR
Turbulence1116
asked this question in
Q&A
-
|
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.
Answered by
EdoAlvarezR
Jan 8, 2024
Replies: 10 comments 46 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.
2 replies
-
|
Also, read this note in the theory section of the docs: LINK |
Beta Was this translation helpful? Give feedback.
-
|
Hi Ed, thank you for the feedback, the note is also helpful. I have switched to ASM, with reference to the link. The simulation did not take super longer that I expected. However, the results from ASM do not seem to improve, I attached both ALM and ASM here:
The cl is even higher for ASM. By observing the particle field and the monitor window, I found that the predicted cl is normal before time step 60, when the propeller wake does not touch the wing wake. Once the propeller wake begins to interact with wing wake, the cl increases immediately. As the wakes travel further downstream, the propeller wakes begin to converge and finally becomes one large wake. In the process, the cl increases again. This observation also applies for ALM. See following shots:
This should not be too related to temporal resolution as I use Euler solver for a faster simulation. The flow is assumed to be inviscid, this should be fine. I also tried to use SFS model, to include taillingboundvortex, the results are still the same. Maybe I have missed something again? I'll go back to your PhD thesis to read more. I attach my Prowim with ASM in the following: #=#############################################################################
# DESCRIPTION
Simulation of prowim, Veldhuis 2005
Full wing span 1.28m
Euler, no SFS, ASM, 40 elements on wing, 20 elements on blade: 20 revos took ? mins and ? particles, killed due to max particle reached.
dt: 0.000113
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 = "prowim_4_asm_" # Name of this simulation
save_path = "prowim_4_asm_" # 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)
# 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 = 4.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, can I turn it off or define manually?
#max_particles += (nsteps+1)*(2*n_wing*(p_per_step+1) + p_per_step)
max_particles = 1000000
# ---------- Vortex sheet parameters (ASM activation)---------------
vlm_vortexsheet = true # Whether to spread the wing circulation as a vortex sheet
vlm_vortexsheet_overlap = 2.125 # 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/100 # VLM-on-VPM smoothing radius of the lifting bound vorticity of wings
sigma_rotor_surf= R/100 # 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)
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.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.SFS_Cd_threelevel_nobackscatter
# 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
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,
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);
# ---------- 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 = true # 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
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=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.
-
|
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.
1 reply
-
|
Thank for the feedback. Yes indeed I missed the sigma_vpm_overwrite, then I refined the overall discretisation accordingly with the reference to the convergence criteria that you described. Now the rule of sigma=0.01*length scale should be satisfied. Additionally, I used RK temporal solver with SFS enabled. The number of particles at the final step is 1.8M. However, the cl is still too high like the following: #=#############################################################################
# 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
Full wing span 1.28m is simulated to remove the reflection plate in the experiment.
Euler, no SFS, ASM, 40 elements on wing, 20 elements on blade, 36 time steps per revo
20 revos took ? mins and ? particles
Total 720 time steps and dt = 0.000113
Too coarse to get reliable results
Rk, SFS(cst Cd), ASM, 100 elements on wing, 50 elements on blade, 72 time steps per revo,
15 revos, RAM is just sufficient (<2M particles), first 100 steps took 1h, 1080 steps = 34h15mins
Total 0.0814s, 1440 steps and dt = 0.00005652777
Still unable to get correct forces
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 = "prowim_4_asm_hr" # Name of this simulation
save_path = "prowim_4_asm_hr" # 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)
# Wing Discretization
n_wing = 100 # 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 = 2.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 = 50 # 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 = 4.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 = 15 # Number of revolutions in simulation
nsteps_per_rev = 72 # 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 = 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 # 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/100 # VLM-on-VPM smoothing radius of the lifting bound vorticity of wings
sigma_rotor_surf= R/100 # 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.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.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
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,
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);
# ---------- 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 = true # 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
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=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.
-
|
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.
24 replies
-
|
Hi, your results of the propeller testing with variation of the incidence of angle of attack is quite fit with my results.
But, the when I ran your ASM code, the sectional lift coefficient distribution was too high as you mentioned. How was it going about analyzing the propeller-wing interaction about PROWIM cases ? Thank you. |
Beta Was this translation helpful? Give feedback.
-
|
Hi Jon! It is odd that CT changes so much as you go to a smaller time step, especially for the CT at 0deg incidence angle. My experience is that CT converges quite nicely for a prop in the J regimes that you are simulating (it might not be as nicely behaved as you approach hover conditions though). In particular, CT in this prop simulation should change very little between coarse and fine time-stepping, as shown in the CT vs Nsteps plot below. The key to have smooth convergence is to change one parameter at a time, while being careful to keep the other parameters constant, as I describe in Sec. 7.2.3 of my dissertation. For instance, as you doubled the number of time steps per rev, it's likely that you also double the sheddings per rev unintentionally.
|
Beta Was this translation helpful? Give feedback.
-
|
Also, check the convergence section in the docs: flow.byu.edu/FLOWUnsteady/theory/convergence |
Beta Was this translation helpful? Give feedback.
-
|
Hello author, what method did you use to accurately determine the numerical values of t1, t2, t3, h1, h2, h3? I hope you can give me a reply. Thank you! |
Beta Was this translation helpful? Give feedback.
-
|
Hi everyone, I am also looking into the t1, t2, t3, h1, h2, h3 tunable parameters. Has anyone found a way to get those? |
Beta Was this translation helpful? Give feedback.
-
|
Hi Inseo, here are results for the propeller at AoAs: |
Beta Was this translation helpful? Give feedback.
3 replies
-
|
I notice here that you are still plotting cl, while the experimental data is for cn (hence, it includes some of the drag force). See my previous comment a month ago. It might not make much of a difference though. It's good to see that at least now the case at 0deg AOA is working well. That's good progress! I wish I had more time to go through this and help you guys out---I go over this forum in my spare time between work and family life. |
Beta Was this translation helpful? Give feedback.
-
|
thanks for the remind, to plot the force components, is it like the way for cls in spanpos? Is it possible to post-process after the simulation is finished? I have turned the xfoil off and used the polar file for the propellers, then the results at AoA =0 are good even with coarse simulation. Thank you for your time, Ed, I really appreciate that! |
Beta Was this translation helpful? Give feedback.
-
|
Hello author, what method did you use to accurately determine the numerical values of t1, t2, t3, h1, h2, h3? I hope you can give me a reply. Thank you! |
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.
2 replies
-
|
Hi Inseo, I think your script is ASM since you have defined the thickness_w = 0.15 and you have vlm_vortexsheet in your uns.run_simulation(). I have turned the xfoil off (which makes a big difference) and have the hubtiploss_correction = ( (0.75, 10, 0.5, 0.05), (1, 1, 1, 1.0) ). A collective pitch at 2 is used. Here is my script which is working fine at AoA=0: #=#############################################################################
# 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
Full wing span 1.28m is simulated to remove the reflection plate in the experiment.
Euler, no SFS, ASM, 40 elements on wing, 20 elements on blade, 36 time steps per revo
20 revos took ? mins and ? particles
Xfoil on
Total 720 time steps and dt = 0.000113, "asm"
Too coarse to get reliable results
Rk, SFS(cst Cd), ASM, 100 elements on wing, 50 elements on blade, 72 time steps per revo (2.5 degress each step)
15 revos, RAM is just sufficient 25G/40G (<2M particles), first 100 steps took 1h, 1080 steps = 34h15mins
Total ?s, 1080 steps and dt = 0.00005652777, "asm_hr"
Still unable to get correct forces
Rk, SFS(cst Cd), ASM, 100 elements on wing, 50 elements on blade, 180 time steps per revo (1 degress each step)
15 revos, tip&hub correction included, RAM at 26G/40G should be fine but simulation seems stuck at first step
for more than 2h, no further progress
Xfoil on, sigma /100
Total 0.061045s, 2700 steps and dt = 0.00002261, "asm_h"
Still unable to get correct forces?
Rk, no SFS, ASM, 50 elements on wing, 20 elements on blade, 36 time steps per revo (10 degress each step),
4 sheds per step, total 144 sheds per revo
15 revos, RAM is sufficient 25G/40G (< 700 000 particles), 3.5h, Cl of wing not converged
Xfoil off, sigma /40
Total 0.061s, 540 steps and dt = 0.000113, hub correction on, denotes "asm_c"
At AoA=0, the sectional lift is well predicted, at AoA=4, still much overpredicted.
Rk, no SFS, ASM, 50 elements on wing, 20 elements on blade, 72 time steps per revo (5 degress each step),
2 sheds per step, total 144 sheds per revo
15 revos, RAM is less than 25G (< 820 000 particles), 8h, Cl of wing not converged
Xfoil off, sigma /40
Total 0.061s, 1080 steps and dt = 0.0000565, hub correction on, denotes "asm_cc"
At AoA=4, still much overpredicted.
Rk, no SFS, ASM, 50 elements on wing, 40 elements on blade, 72 time steps per revo (5 degress each step),
2 sheds per step, total 144 sheds per revo
15 revos, RAM is not sufficient >40G (< 820 000 particles), ?h, Cl of wing not converged?
Xfoil off, sigma /80
Total 0.061s, 1080 steps and dt = 0.0000565, hub correction on, denotes "asm_m"
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 = "prowim_4_asm_cc" # Name of this simulation
save_path = "prowim_4_asm_cc" # 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)
# Wing Discretization
n_wing = 50 # 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 = 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
read_polar = vlm.ap.read_polar2 # What polar reader to use
# 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 = 0.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 = 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 = 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 = 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 # 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/40 # VLM-on-VPM smoothing radius of the lifting bound vorticity of wings
sigma_rotor_surf= R/40 # Rotor-on-VPM smoothing radius of the rotor bound vorticity
vlm_vortexsheet_sigma_tbv = thickness_w*(b/ar) / 40 # 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.5 # 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
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=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)]
#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))Thanks for sharing your script, I'll give it a try :) |
Beta Was this translation helpful? Give feedback.
-
|
Hi, I define the thickness as 0.15, but I did not define the vortex sheet is activating... By the way, how is it going about validating the PROWIM case ? Thank you. |
Beta Was this translation helpful? Give feedback.
-
|
@Turbulence1116, any updates? |
Beta Was this translation helpful? Give feedback.
7 replies
-
|
Hi, Jack. As you know, the airfoil NACA 642-015A and NACA 642-015 are symmetric airfoil. And you can find the difference between NACA 642-015A and NACA 642-015 at the below figure. Maybe, the NACA 642-015A looks there is no sharpen trailing edge.
And the script does not consider the camber effect because the method(Vortex Lattice : 2D). Maybe, the crucial point is the why the upwash and downwash effect are not fit with experimental data and the previous Alvarez's Doctoral thesis ... Thank you. |
Beta Was this translation helpful? Give feedback.
-
|
It's been a while I am busy with something else. @inse0918 if you turn off the Xfoil and use the polar files for the propeller, the results are fit well with experiment (even with very coarse resolution), however, this works with an AoA=0, not for AoA=4 or 8. Maybe the polar file needs to be updated accordingly? |
Beta Was this translation helpful? Give feedback.
-
|
Do you have any issue using ASM with number of element > 50 (say 100)? I tried on 3 high spec work stations, the simulation got stuck at 0 step time, CPU loads are fine, RAM is more than enough. |
Beta Was this translation helpful? Give feedback.
-
|
Hi Jon ! As you mentioned before, the up/downwash effect was fit well with exp. at the 0 deg. like below figure. Then, did you fit the AoA=4 and AoA=8 cases ? or still higher than the exp. data ? I think that the problem is the upwash and downwash effect induced swirl velocity from blade is something wrong or we were mistake in using the scripts. And, as you mentioned, the airfoil polar data should be validated. But, when I study the effect of the tilt angle on propeller performance, the thrust and torque was fit well with exp. data. So, there is no important issue in the airfoil polar data. BTW, the ASM is not appropriate when we check the propeller-wing interaction validation cases, because of a lot of the computational costs. Did you ran the script that I wrote before ? (on Oct 9) and did you compared it ? I'm still running the PROWIM cases using ALM not a ASM ... Thank you. |
Beta Was this translation helpful? Give feedback.
-
|
Hi Inseo, I could not fit for the AOA=4 and AOA=8, they are still overpredicted. Yeah maybe not from the polar file, same observations with me for isolated beaver propeller: |
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.
0 replies
-
|
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.
6 replies
-
|
Hi ! Nice work !, I checked your script that can run the prowim code just now. But, the calculation does not going well.
The first problem is the 'No sound speed has been provided.' So, I delete the 'sound_spd=nothing' line. Next, the error does not pop uped. But, the calculation also does not going.
Did you check this ? Thank you. |
Beta Was this translation helpful? Give feedback.
-
|
Same with me. Simulation just got stuck at time step 0, no issue with RAM and CPUs are fully charged. |
Beta Was this translation helpful? Give feedback.
-
|
Hi, @inse0918 , @Turbulence1116 , this code works fine on my workstation (although it requires a lot of computational effort). |
Beta Was this translation helpful? Give feedback.
-
|
Hey, can you pls let me know the specs of your workstation. Mine has CPU i7-13700, 32G RAM. |
Beta Was this translation helpful? Give feedback.
-
|
Laptop : 12th Gen. Intel(R) Core(TM) i5-12500H 2.5 GHz / 40 GB RAM Inseo |
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.
0 replies
-
|
Closing this thread and continuing the discussion in #129. |
Beta Was this translation helpful? Give feedback.
1 reply
-
|
Thank you a lot for the help! |
Beta Was this translation helpful? Give feedback.
Answer selected by
EdoAlvarezR
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Closing this thread and continuing the discussion in #129.