import FLOWUnsteady as uns import FLOWVLM as vlm import FLOWVPM as vpm run_name = "heaveingPWI-example" # Name of this simulation BasePath = "/home/inse0918/.julia/packages/FLOWUnsteady/Jl5A6/Results/PWI" save_path = run_name # Where to save this simulation paraview = false # Whether to visualize with Paraview # ----------------- GEOMETRY PARAMETERS ---------------------------------------- # Wing geometry b = 1.28 # (m) span length ar = 1.28/0.24 # 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 frequency = 100.0 # (Hz) oscillation frequency in heaving motion amplitude = 5.0 # (deg) AOA amplitude in heaving motion # Rotor geometry rotor_file = "apc10x7.csv" # Rotor geometry data_path = uns.def_data_path # Path to rotor database pitch = 0.0 # (deg) collective pitch of blades xfoil = false # Whether to run XFOIL ncrit = 9 # Turbulence criterion for XFOIL # Read radius of this rotor and number of blades R, Rhub, B = uns.read_rotor(rotor_file; data_path=data_path)[[1, 2, 3]] # Vehicle assembly AOAwing = 0 # (deg) wing angle of attack spanpos = [-0.3, 0.3] # Semi-span position of each rotor, 2*y/b xpos = [0, 0] # x-position of rotors relative to LE, x/c zpos = [0, 0] # z-position of rotors relative to LE, z/c CWs = [true, false] # Clockwise rotation for each rotor PhaseLagAng = [0, 0] # angle of the phase lag angle nrotors = length(spanpos) # Number of rotors # Discretization n_rotor = 30 # Number of blade elements per blade r_rotor = 1/10 # Geometric expansion of elements n_wing = 50 r_wing = 1.0 # 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 # Rotor operation J = 0.85 # Advance ratio Vref/(nD) RPM = 14800 # RPM # RPMSeries = [] # Freestream magVinf = J * (RPM / 60) * (2 * R) # (m/s) freestream velocity rho = 1.225 # (kg/m^3) air density mu = 1.81e-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 # Rotor operation 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 J : $(round(J, digits=1)) RPM: $(RPM) Matip: $(round(Matip, digits=3)) ReD: $(round(ReD, 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 = 10 # Time steps per revolution 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 = 4 # 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_viscous = vpm.Inviscid() # VPM viscous diffusion scheme vpm_SFS = vpm.SFS_none # VPM LES subfilter-scale model 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_wing = [0.2018, 0.0, 0.0] # New position Oaxis_wing = uns.gt.rotation_matrix2(0, -AOAwing, 0) # New orientation vlm.setcoordsystem(wing, O_wing, Oaxis_wing) # Generate rotors rotors = vlm.Rotor[] for ri in 1:nrotors # 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] x = xpos[ri] z = zpos[ri] # # Account for angle of attack of wing # nrm = sqrt(x^2 + z^2) # x = nrm*cosd(AOAwing) # z = -nrm*sind(AOAwing) # Translate rotor to its position along wing O_rotor = [x, y, z] # New position Oaxis_rotor = uns.gt.rotation_matrix2(PhaseLagAng[ri], 0, 180) # New orientation vlm.setcoordsystem(rotor, O_rotor, Oaxis_rotor) push!(rotors, rotor) end # -------- Generate vehicle println("Generating vehicle...") heavingsystem = vlm.WingSystem() # Tilting system vlm.addwing(heavingsystem, "FrontWing", wing) # System of all FLOWVLM objects system = vlm.WingSystem() vlm.addwing(system, "Heaving", heavingsystem) for (ri, rotor) in enumerate(rotors) vlm.addwing(system, "Rotor$(ri)", rotor) end tilting_systems = (heavingsystem, ); # Tilting systems # Systems of rotors rotor_systems = (rotors, ); # System solved through VLM solver vlm_system = vlm.WingSystem() vlm.addwing(vlm_system, "Heaving", heavingsystem) # System that will shed a VPM wake wake_system = vlm.WingSystem() # System that will shed a VPM wake vlm.addwing(wake_system, "Heaving", heavingsystem) # 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; tilting_systems=tilting_systems, 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) = [0, 0, 0] # <---- Vehicle is traveling in the -x direction # Angle of the vehicle over time anglevehicle(t) = zeros(3) # Control inputs angle_frontwing(t) = [0, amplitude*sin(2*pi*frequency * t*ttot), 0] # Tilt angle of front wing # RPM control input over time (RPM over `RPMref`) RPMcontrol(t) = 1.0 angles = (angle_frontwing,) # 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-n64015a-il-1000000.csv" ) L_dir(t) = [0, 0, 1] # Direction of lift D_dir(t) = [1, 0, 0] # Direction of drag # 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...") # Run simulation uns.run_simulation(simulation, nsteps; # ----- SIMULATION OPTIONS ------------- Vinf=Vinf, rho=rho, mu=mu, sound_spd=speedofsound, # ----- 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 ); import FLOWUnsteady: cross, dot, norm, plt, @L_str 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) ) Shat = [0, 1, 0] # Spanwise direction Dhat = [cosd(AOA), 0.0, sind(AOA)] # Direction of drag Lhat = cross(Dhat, Shat) # Direction of lift # Integrated moment decomposed into rolling, pitching, and yawing moments 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(BasePath*"pos.dat", spanpos) writedlm(BasePath*"cl.dat", cls) writedlm(BasePath*"cd.dat", cds)