christianhauschel
commented
5 months ago n0012.csv Here the full code:
# ==============================================================================
# Setup
# ==============================================================================
import FLOWUnsteady as uns
import FLOWUnsteady: vlm, vpm, gt, Im
fidelity = :low
delete_previous_sim = true
nsteps_restart = -1
n_revs = 2 # number of revolutions to simulate
restart_file = nothing
run_name = "sim_12"
nsteps_save = 5
if fidelity == :low
nsteps_save *= 1
elseif fidelity == :mid
nsteps_save *= 5
else
nsteps_save *= 10
end
stride_fluiddomain = max(25, nsteps_save)
maneuver_plot = true
paraview = true # Whether to visualize with Paraview
dir_base = @__DIR__
save_path = joinpath("out", run_name)
if delete_previous_sim
if isdir(save_path)
rm(save_path; force=true, recursive=true)
end
end
mkpath(save_path)
# ==============================================================================
# Operating Point
# ==============================================================================
Vinf(X, t) = [-18.5, 0.0001, 0.0001]
ρ = 1.2
μ = 0.000018
c = 340
Vcruise = 0.0 # (m/s) cruise speed (reference)
RPMh_w = 2000 # (RPM) hover RPM
ttot = n_revs / (RPMh_w / 60) # (s) total time to perform maneuver
# ==============================================================================
# Vehicle
# ==============================================================================
# Drone parameters
pos_LE = [0.0, 0.0, 0.0]
pos_G_LE = [-0.093, 0.0, -0.055] # center of gravity w.r.t. LE
tiphub_expansion = 1 / 10
m = 22.636 # (kg) mass
# Wing
AR = 5.0
b = 4.95
c_tip = 0.19200
c_hub = 0.32
AR = b / c_tip
n_span = 40
pos = [0, 0.8251649999999999, b / 2] / (b / 2)
clen = [c_hub / c_tip, c_hub / c_tip, c_tip / c_tip]
twist = -[-2.540, -2.540, -2.650]
sweep = [0.0, 0.0]
dihed = [0.0, 0.0]
wing = uns.vlm.complexWing(b, AR, n_span, pos, clen, twist, sweep, dihed; symmetric=true)
wing.Oaxis = gt.rotation_matrix(180, 0, 0)
wing.O = -pos_G_LE
system = uns.vlm.WingSystem()
vlm_system = uns.vlm.WingSystem()
wake_system = uns.vlm.WingSystem()
uns.vlm.addwing(vlm_system, "Wing", wing)
uns.vlm.addwing(wake_system, "SolveVLM", vlm_system)
grids = []
vehicle = uns.VLMVehicle(
system;
vlm_system=vlm_system,
wake_system=wake_system,
grids=grids
)
# ==============================================================================
# Maneuver
# ==============================================================================
angle_vehicle_x(t) = 20
angle_vehicle_y(t) = 10
angle_vehicle_z(t) = 30
angle_vehicle(t) = [angle_vehicle_x(t), angle_vehicle_y(t), angle_vehicle_z(t)]
Vvehicle(t) = [0, 0, 0]
RPM = ()
angle_tilts = ()
maneuver = uns.KinematicManeuver(angle_tilts, RPM, Vvehicle, angle_vehicle)
# ==============================================================================
# Simulation
# ==============================================================================
use_variable_pitch = true # Whether to use variable pitch in cruise
tstart = 0.00 * ttot # (s) start simulation at this point in time
tquit = 1.00 * ttot # (s) end simulation at this point in time
start_kinmaneuver = true # If true, it starts the maneuver with the
# velocity and angles of tstart.
# If false, starts with velocity=0
# and angles as initiated by the geometry
# -------------------------------------
# Solver
# -------------------------------------
# Aerodynamic solver
VehicleType = uns.UVLMVehicle # Unsteady solver
# VehicleType = uns.QVLMVehicle # Quasi-steady solver
# Time parameters
if fidelity == :low
nsteps_per_rev = 36 # 10 °/step
elseif fidelity == :mid
nsteps_per_rev = 72 # 5 °/step
else
nsteps_per_rev = 360 # 1 °/step
end
resolution_deg = 360 / nsteps_per_rev
nsteps = Int(round(RPMh_w / 60 * 360 * ttot / resolution_deg))
# nsteps = 250 # Time steps for entire maneuver
dt = ttot / nsteps # (s) time step
# VPM particle shedding
p_per_step = 4 # Sheds per time step
shed_starting = fidelity == :low || fidelity == :mid ? false : 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
# Regularization of embedded vorticity
r_tip = 0.2
sigma_rotor_surf = fidelity == :low || fidelity == :mid ? r_tip / 10 : r_tip / 80
sigma_vlm_surf = b / 400 # sigma wing
lambda_vpm = 2.125 # VPM core overlap
# VPM smoothing radius
sigma_vpm_overwrite = lambda_vpm * 2 * pi * r_tip / (nsteps_per_rev * p_per_step)
sigmafactor_vpmonvlm = fidelity == :low || fidelity == :mid ? 1 : 5.5 # Shrink particles by this factor when
# calculating VPM-on-VLM/Rotor induced velocities
# Rotor solver
vlm_rlx = 0.2 # VLM relaxation <-- this also applied to rotors
hubtiploss_correction = vlm.hubtiploss_correction_prandtl # Hub and tip correction
# ==============================================================================
# Wing Solver Settings
# ==============================================================================
# Wing solver: actuator surface model (ASM)
vlm_vortexsheet = fidelity != :low ? true : false # Whether to spread the wing surface vorticity as a vortex sheet (activates ASM)
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
chordwing_mean = 0.3
chord = chordwing_mean
thickness_norm = 0.13698
vlm_vortexsheet_sigma_tbv = thickness_norm * chord / 100 # Size of particles in trailing bound vortices
# vlm_vortexsheet_sigma_tbv = sigma_vpm_overwrite
# How many particles to preallocate for the vortex sheet
vlm_vortexsheet_maxstaticparticle = vlm_vortexsheet == false ? nothing : 6000000
# Wing solver: force calculation
KJforce_type = "regular" # KJ force evaluated at middle of bound vortices_vortexsheet also true
include_trailingboundvortex = false # Include trailing bound vortices in force calculations
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
# VPM solver
vpm_integration = fidelity == :mid || fidelity == :high ? vpm.rungekutta3 : vpm.euler
# Viscosous particles
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 LES subfilter-scale model
if fidelity == :high
# vpm_SFS = vpm.DynamicSFS(vpm.Estr_fmm, vpm.pseudo3level_positive;
# alpha=0.999, maxC=1.0,
# clippings=[vpm.clipping_backscatter],
# controls=[vpm.control_directional, vpm.control_magnitude]
# )
vpm_SFS = vpm.SFS_Cd_twolevel_nobackscatter()
else
vpm_SFS = vpm.SFS_none
end
if VehicleType == uns.QVLMVehicle
uns.vlm.VLMSolver._mute_warning(true)
end
# -------------------------------------
# Acoustics
# -------------------------------------
save_wopwop_input = false
save_init_plots = false
if maneuver_plot
uns.plot_maneuver(maneuver; ti=0, tf=tquit, save_path=save_path, size_factor=2 / 3)
end
# ==============================================================================
# Monitors
# ==============================================================================
# -------------------- WING MONITORS ---------------------------------------
# Reference parameters for calculating coefficients
# NOTE: make b, ar, and qinf equals to 1.0 to obtain dimensional force
b_ref, ar_ref = 1.0, 1.0
qinf = 1.0
Jref = 1.0
# Force axis labels
CL_lbl = "Lift (N)"
CD_lbl = "Drag (N)"
# Directions of force components
L_dir = [0, 0, 1]
D_dir = [-1, 0, 0]
calc_aerodynamicforce_fun_wing2 = uns.generate_calc_aerodynamicforce(;
add_parasiticdrag=include_parasiticdrag,
add_skinfriction=add_skinfriction,
airfoilpolar="n0012.csv",
)
monitor_name = "Wing"
wingL = vlm.get_wing(vehicle.vlm_system, "Wing")
monitor_wing = uns.generate_monitor_wing(
wingL, Vinf, b_ref, ar_ref,
ρ, qinf, nsteps;
calc_aerodynamicforce_fun=calc_aerodynamicforce_fun_wing2,
save_path=save_path,
run_name=monitor_name,
figname=monitor_name,
CL_lbl=CL_lbl,
CD_lbl=CD_lbl,
L_dir=L_dir,
D_dir=D_dir,
# monitor_optargs_wing...
)
runtime_function = monitor_wing
# -------------------- OTHER MONITORS --------------------------------------
# State-variable monitor
statevariable_monitor = uns.generate_monitor_statevariables(; save_path=save_path)
push!(monitors, statevariable_monitor)
# Global enstrophy monitor (numerical stability)
monitor_enstrophy = uns.generate_monitor_enstrophy(; save_path=save_path)
push!(monitors, monitor_enstrophy)
# Monitor of SFS model coefficient Cd
if vpm_SFS != vpm.SFS_none
monitor_Cd = uns.generate_monitor_Cd(; save_path=save_path)
push!(monitors, monitor_Cd)
end
# -------------------- CONCATENATE MONITORS --------------------------------
monitors = uns.concatenate(monitors)
# Reference parameters
Vref = Vcruise # Reference velocity to scale maneuver by
RPMref = RPMh_w # Reference RPM to scale maneuver by
# Initial conditions
t0 = tstart / ttot * start_kinmaneuver # Time at start of simulation
Vinit = Vref * maneuver.Vvehicle(t0) # Initial vehicle velocity
Winit = pi / 180 * (maneuver.anglevehicle(t0 + 1e-12) - maneuver.anglevehicle(t0)) / (ttot * 1e-12) # Initial angular velocity
# Define simulation
simulation = uns.Simulation(
vehicle,
maneuver,
Vref,
RPMref,
ttot;
Vinit=Vinit,
Winit=Winit,
t=tstart
);
# Maximum number of particles (for pre-allocating memory)
max_particles = ceil(Int, (nsteps + 2) * (2 * vlm.get_m(vehicle.wake_system) * (p_per_step + 1) + p_per_step))
max_particles = tquit != Inf ? ceil(Int, max_particles * (tquit - tstart) / ttot) : max_particles
max_particles = min(10000000, max_particles)
max_particles = VehicleType == uns.QVLMVehicle ? 10000 : max_particles
omit_shedding = []
# ==============================================================================
# RUN SIMULATION
# ==============================================================================
@time uns.run_simulation(
simulation, nsteps;
# ----- SIMULATION OPTIONS -------------
Vinf=Vinf,
rho=ρ, mu=μ, tquit=tquit, sound_spd=c,
# ----- SOLVERS OPTIONS ----------------
p_per_step=p_per_step,
max_particles=max_particles,
max_static_particles=vlm_vortexsheet_maxstaticparticle,
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_vortexsheet=vlm_vortexsheet,
vlm_vortexsheet_overlap=vlm_vortexsheet_overlap,
vlm_vortexsheet_distribution=vlm_vortexsheet_distribution,
vlm_vortexsheet_sigma_tbv=vlm_vortexsheet_sigma_tbv,
vlm_rlx=vlm_rlx,
hubtiploss_correction=hubtiploss_correction,
shed_starting=shed_starting,
shed_unsteady=shed_unsteady,
unsteady_shedcrit=unsteady_shedcrit,
omit_shedding=omit_shedding,
extra_runtime_function=runtime_function,
# ----- RESTART OPTIONS -----------------
restart_vpmfile=restart_file,
# ----- OUTPUT OPTIONS ------------------
save_path=save_path,
run_name=run_name,
save_wopwopin=save_wopwop_input, # <--- Generates input files for PSU-WOPWOP noise analysis
# save_code=splitdir(@__FILE__)[1]
nsteps_save=nsteps_save,
nsteps_restart=nsteps_restart,
prompt=false,
create_savepath=false,
);
I want to simulate a wing in cruise flight. I have set a vehicle rotation maneuver as follows:
For some reason, the rotation is not visible in Paraview and not applied in the FU simulation.
Why could this be the case?