-
Notifications
You must be signed in to change notification settings - Fork 0
/
generate_piecewise_inputs.jl
130 lines (104 loc) · 4.69 KB
/
generate_piecewise_inputs.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
import PyPlot
PyPlot.pygui(true)
import FLOWMath
# Inputs
chord_length = 1.0 #m
windspeed = 15.0 #m/s
NchordFlow = 1.5
steady_state_time = 1/windspeed * NchordFlow/chord_length# sec, at 15m/s for a 1m chord, 1 chords is .67 seconds
# Potentially repeating/array of numbers
sin_duration = 2.0 # sec, for 15m/s and 1m chord, then 1/15 = 0.067 seconds per chord, and if (below) we are doing 20 chords in a full cycle, if we did 5 cycles, then that's 100 chords of flow, or 6.7 seconds, so just say 10 seconds to match up with the sim steps?
N_sin_points = 100
frequency_rot = 0.75 #hz -> 1 hz means we go from 0 to amplitude to negative amplitude to zero in one second. for 15m/s and 1m chord, if we want 5 chords of flow when it hits the high side (1/4 stroke), then 1/15 = 0.067 seconds per chord, multiply by 5 = 0.33 seconds, so want full cycle in 4x that time, so 1 cycle in 1.33 seconds, 1/1.33 = 0.75 cycles/second
frequency_edge = 0.75 #hz
frequency_flap = 0.75 #hz
rot_amp_deg = 10.0 # deg
edge_amp_m = 0.25 # m
flap_amp_m = 0.125 # m
N_sample_out = 50
# Get to steady state
N_steady_state = 10
time_array = collect(LinRange(0.0, steady_state_time, N_steady_state))
theta_deg_array = collect(LinRange(0.0, 0.0, N_steady_state))
edge_offset_array = collect(LinRange(0.0, 0.0, N_steady_state))
flap_offset_array = collect(LinRange(0.0, 0.0, N_steady_state))
# Define motion, could do multiple input in a loop to get multiple different frequencies, amplitudes in a run, etc
# Intermediate time array
time_sin = collect(LinRange(1e-6, sin_duration, N_sin_points))
# Define sinusoidal rotation
sin_rotation_deg = rot_amp_deg*sin.(2*pi*frequency_rot*time_sin)
# Define sinusoidal edge motion
sin_edge_m = edge_amp_m*sin.(2*pi*frequency_edge*time_sin)
# define sinusoidal flap motion
sin_flap_m = flap_amp_m*sin.(2*pi*frequency_flap*time_sin)
time_array = [time_array; time_array[end].+time_sin]
theta_deg_array = [theta_deg_array; sin_rotation_deg]
edge_offset_array = [edge_offset_array; sin_edge_m]
flap_offset_array = [flap_offset_array; sin_flap_m]
span_offset_array = zeros(length(time_array))
# Spline
time_out = LinRange(time_array[1],time_array[end],N_sample_out)
theta_deg_out = FLOWMath.akima(time_array,theta_deg_array,time_out)
edge_offset_out = FLOWMath.akima(time_array,edge_offset_array,time_out)
flap_offset_out = FLOWMath.akima(time_array,flap_offset_array,time_out)
# Plot to verify
PyPlot.figure()
PyPlot.plot(time_array,theta_deg_array,"r",label="Theta (Deg)")
PyPlot.plot(time_array,edge_offset_array,"g",label="Edge offset (m)")
PyPlot.plot(time_array,flap_offset_array,"b",label="Flap offset (m)")
PyPlot.plot(time_out,theta_deg_out,"r--",label="Theta (Deg)")
PyPlot.plot(time_out,edge_offset_out,"g--",label="Edge offset (m)")
PyPlot.plot(time_out,flap_offset_out,"b--",label="Flap offset (m)")
PyPlot.legend()
# theta_deg_out = cumsum(theta_deg_out)
# edge_offset_out = cumsum(edge_offset_out)
# flap_offset_out = cumsum(flap_offset_out)
# PyPlot.figure()
# for itime = 1:length(time_out)-1
# myx0 = [-chord_length/2,chord_length/2]
# myy0 = [0,0]
# myx = myx0*cosd(theta_deg_out[itime])-myy0*sind(theta_deg_out[itime])
# myy = myx0*sind(theta_deg_out[itime])+myy0*cosd(theta_deg_out[itime])
# PyPlot.clf()
# PyPlot.title("Time: $(time_out[itime])")
# PyPlot.plot(myx.+edge_offset_out[itime],myy.+flap_offset_out[itime])
# PyPlot.xlim((-1.5,5.5))
# PyPlot.ylim((-1.5,5.5))
# PyPlot.savefig("./figs/motionframe$(lpad(itime,3,"0")).jpg",transparent = true)
# sleep(0.001)
# if time_out[itime]>0.693
# break
# end
# end
# cd("./figs/")
# run(`ffmpeg -i motionframe%03d.jpg -vcodec libx264 -vf "pad=ceil(iw/2)*2:ceil(ih/2)*2" -r 24 -y -an -pix_fmt yuv420p airfoilmotion.mp4`)
# run(`ffmpeg -ss 0 -t 26 -i airfoilmotion.mp4 \
# -vf "fps=10,scale=320:-1:flags=lanczos,split[s0][s1];[s0]palettegen[p];[s1][p]paletteuse" \
# -loop 0 airfoilmotion.gif`)
# cd("../")
# Now Create the output definition
theta_deg_out = diff(theta_deg_out)
edge_offset_out = diff(edge_offset_out)
flap_offset_out = diff(flap_offset_out)
println(" mesh_motion:
- name: arbitrary_motion_airfoil
mesh_parts:
- airfoil-HEX
frame: non_inertial
motion:")
# Rotations
for itime = 1:length(time_out)-1
println(" - type: rotation
angle: $(theta_deg_out[itime])
start_time: $(time_out[itime]+1e-6)
end_time: $(time_out[itime+1])
axis: [0.0, 0.0, 1.0]
origin: [0.0, 0.0, 0.0]")
end
# Displacements
for itime = 1:length(time_out)-1
println(" - type: translation
start_time: $(time_out[itime]+1e-6)
end_time: $(time_out[itime+1])
displacement: [$(flap_offset_out[itime]), $(edge_offset_out[itime]), 0.0]")
end