| Line | Exclusive | Inclusive | Code |
|---|---|---|---|
| 1 | #= MIT License | ||
| 2 | |||
| 3 | Copyright (c) 2020, 2021, 2022 Uwe Fechner | ||
| 4 | |||
| 5 | Permission is hereby granted, free of charge, to any person obtaining a copy | ||
| 6 | of this software and associated documentation files (the "Software"), to deal | ||
| 7 | in the Software without restriction, including without limitation the rights | ||
| 8 | to use, copy, modify, merge, publish, distribute, sublicense, and/or sell | ||
| 9 | copies of the Software, and to permit persons to whom the Software is | ||
| 10 | furnished to do so, subject to the following conditions: | ||
| 11 | |||
| 12 | The above copyright notice and this permission notice shall be included in all | ||
| 13 | copies or substantial portions of the Software. | ||
| 14 | |||
| 15 | THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR | ||
| 16 | IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, | ||
| 17 | FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE | ||
| 18 | AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER | ||
| 19 | LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, | ||
| 20 | OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE | ||
| 21 | SOFTWARE. =# | ||
| 22 | |||
| 23 | #= Model of a kite-power system in implicit form: residual = f(y, yd) | ||
| 24 | |||
| 25 | This model implements a 3D mass-spring system with reel-out. It uses six tether segments (the number can be | ||
| 26 | configured in the file data/settings.yaml). The kite is modelled using 4 point masses and 3 aerodynamic | ||
| 27 | surfaces. The spring constant and the damping decrease with the segment length. The aerodynamic kite forces | ||
| 28 | are acting on three of the four kite point masses. | ||
| 29 | |||
| 30 | Four point kite model, included from KiteModels.jl. | ||
| 31 | |||
| 32 | Scientific background: http://arxiv.org/abs/1406.6218 =# | ||
| 33 | |||
| 34 | # Array of connections of bridlepoints. | ||
| 35 | # First point, second point, unstressed length. | ||
| 36 | const SPRINGS_INPUT = [0. 1. 150. | ||
| 37 | 1. 2. -1. # s1, p7, p8 | ||
| 38 | 4. 2. -1. # s2, p10, p8 | ||
| 39 | 4. 5. -1. # s3, p10, p11 | ||
| 40 | 3. 4. -1. # s4, p9, p10 | ||
| 41 | 5. 1. -1. # s5, p11, p7 | ||
| 42 | 4. 1. -1. # s6, p10, p7 | ||
| 43 | 3. 5. -1. # s7, p9, p11 | ||
| 44 | 5. 2. -1. # s8, p11, p8 | ||
| 45 | 2. 3. -1.] # s9, p8, p9 | ||
| 46 | |||
| 47 | # struct, defining the phyical parameters of one spring | ||
| 48 | @with_kw struct Spring{I, S} | ||
| 49 | p1::I = 1 # number of the first point | ||
| 50 | p2::I = 2 # number of the second point | ||
| 51 | length::S = 1.0 # current unstressed spring length | ||
| 52 | c_spring::S = 1.0 # spring constant [N/m] | ||
| 53 | damping::S = 0.1 # damping coefficent [Ns/m] | ||
| 54 | end | ||
| 55 | |||
| 56 | const SP = Spring{Int16, Float64} | ||
| 57 | const KITE_PARTICLES = 4 | ||
| 58 | const KITE_SPRINGS = 9 | ||
| 59 | const KITE_ANGLE = 3.83 # angle between the kite and the last tether segment due to the mass of the control pod | ||
| 60 | const PRE_STRESS = 0.9998 # Multiplier for the initial spring lengths. | ||
| 61 | const KS = deg2rad(16.565 * 1.064 * 0.875 * 1.033 * 0.9757 * 1.083) # max steering | ||
| 62 | const DRAG_CORR = 0.93 # correction of the drag for the 4-point model | ||
| 63 | function zero(::Type{SP}) | ||
| 64 | SP(0,0,0,0,0) | ||
| 65 | end | ||
| 66 | |||
| 67 | """ | ||
| 68 | mutable struct KPS4{S, T, P, Q, SP} <: AbstractKiteModel | ||
| 69 | |||
| 70 | State of the kite power system, using a 4 point kite model. Parameters: | ||
| 71 | - S: Scalar type, e.g. SimFloat | ||
| 72 | In the documentation mentioned as Any, but when used in this module it is always SimFloat and not Any. | ||
| 73 | - T: Vector type, e.g. MVector{3, SimFloat} | ||
| 74 | - P: number of points of the system, segments+1 | ||
| 75 | - Q: number of springs in the system, P-1 | ||
| 76 | - SP: struct type, describing a spring | ||
| 77 | Normally a user of this package will not have to access any of the members of this type directly, | ||
| 78 | use the input and output functions instead. | ||
| 79 | |||
| 80 | $(TYPEDFIELDS) | ||
| 81 | """ | ||
| 82 | @with_kw mutable struct KPS4{S, T, P, Q, SP} <: AbstractKiteModel | ||
| 83 | "Reference to the settings struct" | ||
| 84 | set::Settings = se() | ||
| 85 | "Reference to the KCU model (Kite Control Unit as implemented in the package KitePodModels" | ||
| 86 | kcu::KCU = KCU() | ||
| 87 | "Reference to the atmospheric model as implemented in the package AtmosphericModels" | ||
| 88 | am::AtmosphericModel = AtmosphericModel() | ||
| 89 | "Reference to winch model as implemented in the package WinchModels" | ||
| 90 | wm::AbstractWinchModel = AsyncMachine() | ||
| 91 | "Iterations, number of calls to the function residual!" | ||
| 92 | iter:: Int64 = 0 | ||
| 93 | "Function for calculation the lift coefficent, using a spline based on the provided value pairs." | ||
| 94 | calc_cl = Spline1D(se().alpha_cl, se().cl_list) | ||
| 95 | "Function for calculation the drag coefficent, using a spline based on the provided value pairs." | ||
| 96 | calc_cd = Spline1D(se().alpha_cd, se().cd_list) | ||
| 97 | "wind vector at the height of the kite" | ||
| 98 | v_wind::T = zeros(S, 3) | ||
| 99 | "wind vector at reference height" | ||
| 100 | v_wind_gnd::T = zeros(S, 3) | ||
| 101 | "wind vector used for the calculation of the tether drag" | ||
| 102 | v_wind_tether::T = zeros(S, 3) | ||
| 103 | "apparent wind vector at the kite" | ||
| 104 | v_apparent::T = zeros(S, 3) | ||
| 105 | "drag force of kite and bridle; output of calc_aero_forces!" | ||
| 106 | drag_force::T = zeros(S, 3) | ||
| 107 | "lift force of the kite; output of calc_aero_forces!" | ||
| 108 | lift_force::T = zeros(S, 3) | ||
| 109 | "spring force of the current tether segment, output of calc_particle_forces!" | ||
| 110 | spring_force::T = zeros(S, 3) | ||
| 111 | "last winch force" | ||
| 112 | last_force::T = zeros(S, 3) | ||
| 113 | "a copy of the residual one (pos,vel) for debugging and unit tests" | ||
| 114 | res1::SVector{P, KVec3} = zeros(SVector{P, KVec3}) | ||
| 115 | "a copy of the residual two (vel,acc) for debugging and unit tests" | ||
| 116 | res2::SVector{P, KVec3} = zeros(SVector{P, KVec3}) | ||
| 117 | "a copy of the actual positions as output for the user" | ||
| 118 | pos::SVector{P, KVec3} = zeros(SVector{P, KVec3}) | ||
| 119 | "velocity vector of the kite" | ||
| 120 | vel_kite::T = zeros(S, 3) | ||
| 121 | "unstressed segment length [m]" | ||
| 122 | segment_length::S = 0.0 | ||
| 123 | "lift coefficient of the kite, depending on the angle of attack" | ||
| 124 | param_cl::S = 0.2 | ||
| 125 | "drag coefficient of the kite, depending on the angle of attack" | ||
| 126 | param_cd::S = 1.0 | ||
| 127 | "azimuth angle in radian; inital value is zero" | ||
| 128 | psi::S = zero(S) | ||
| 129 | "depower angle [deg]" | ||
| 130 | alpha_depower::S = 0.0 | ||
| 131 | "relative start time of the current time interval" | ||
| 132 | t_0::S = 0.0 | ||
| 133 | "reel out speed of the winch" | ||
| 134 | v_reel_out::S = 0.0 | ||
| 135 | "reel out speed at the last time step" | ||
| 136 | last_v_reel_out::S = 0.0 | ||
| 137 | "unstretched tether length" | ||
| 138 | l_tether::S = 0.0 | ||
| 139 | "air density at the height of the kite" | ||
| 140 | rho::S = 0.0 | ||
| 141 | "actual relative depower setting, must be between 0 .. 1.0" | ||
| 142 | depower::S = 0.0 | ||
| 143 | "actual relative steering setting, must be between -1.0 .. 1.0" | ||
| 144 | steering::S = 0.0 | ||
| 145 | "multiplier for the stiffniss of tether and bridle" | ||
| 146 | stiffness_factor::S = 1.0 | ||
| 147 | "initial masses of the point masses" | ||
| 148 | initial_masses::MVector{P, S} = ones(P) | ||
| 149 | "current masses, depending on the total tether length" | ||
| 150 | masses::MVector{P, S} = zeros(P) | ||
| 151 | "vector of the springs, defined as struct" | ||
| 152 | springs::MVector{Q, SP} = zeros(SP, Q) | ||
| 153 | "vector of the forces, acting on the particles" | ||
| 154 | forces::SVector{P, KVec3} = zeros(SVector{P, KVec3}) | ||
| 155 | "synchronous speed of the motor/ generator" | ||
| 156 | sync_speed::S = 0.0 | ||
| 157 | "x vector of kite reference frame" | ||
| 158 | x::T = zeros(S, 3) | ||
| 159 | "y vector of kite reference frame" | ||
| 160 | y::T = zeros(S, 3) | ||
| 161 | "z vector of kite reference frame" | ||
| 162 | z::T = zeros(S, 3) | ||
| 163 | end | ||
| 164 | |||
| 165 | @inline @inbounds function norm(vec::SVector{3, Float64}) | ||
| 166 | sqrt(vec[1]*vec[1]+vec[2]*vec[2]+vec[3]*vec[3]) | ||
| 167 | end | ||
| 168 | @inline @inbounds function norm(vec::MVector{3, Float64}) | ||
| 169 | sqrt(vec[1]*vec[1]+vec[2]*vec[2]+vec[3]*vec[3]) | ||
| 170 | end | ||
| 171 | |||
| 172 | """ | ||
| 173 | clear!(s::KPS4) | ||
| 174 | |||
| 175 | Initialize the kite power model. | ||
| 176 | """ | ||
| 177 | function clear!(s::KPS4) | ||
| 178 | s.t_0 = 0.0 # relative start time of the current time interval | ||
| 179 | s.v_reel_out = 0.0 | ||
| 180 | s.last_v_reel_out = 0.0 | ||
| 181 | s.v_wind_gnd .= [s.set.v_wind, 0, 0] # wind vector at reference height | ||
| 182 | s.v_wind_tether .= [s.set.v_wind, 0, 0] | ||
| 183 | s.v_apparent .= [s.set.v_wind, 0, 0] | ||
| 184 | height = sin(deg2rad(s.set.elevation)) * (s.set.l_tether) | ||
| 185 | s.v_wind .= s.v_wind_gnd * calc_wind_factor(s.am, height) | ||
| 186 | |||
| 187 | s.l_tether = s.set.l_tether | ||
| 188 | s.segment_length = s.l_tether / s.set.segments | ||
| 189 | init_masses!(s) | ||
| 190 | init_springs!(s) | ||
| 191 | for i in 1:s.set.segments + KiteModels.KITE_PARTICLES + 1 | ||
| 192 | s.forces[i] .= zeros(3) | ||
| 193 | end | ||
| 194 | s.drag_force .= [0.0, 0, 0] | ||
| 195 | s.lift_force .= [0.0, 0, 0] | ||
| 196 | s.rho = s.set.rho_0 | ||
| 197 | s.calc_cl = Spline1D(s.set.alpha_cl, s.set.cl_list) | ||
| 198 | s.calc_cd = Spline1D(s.set.alpha_cd, s.set.cd_list) | ||
| 199 | s.kcu.depower = s.set.depower/100.0 | ||
| 200 | s.kcu.set_depower = s.kcu.depower | ||
| 201 | KiteModels.set_depower_steering!(s, get_depower(s.kcu), get_steering(s.kcu)) | ||
| 202 | end | ||
| 203 | |||
| 204 | function KPS4(kcu::KCU) | ||
| 205 | s = KPS4{SimFloat, KVec3, kcu.set.segments+KITE_PARTICLES+1, kcu.set.segments+KITE_SPRINGS, SP}() | ||
| 206 | s.set = kcu.set | ||
| 207 | s.kcu = kcu | ||
| 208 | s.calc_cl = Spline1D(s.set.alpha_cl, s.set.cl_list) | ||
| 209 | s.calc_cd = Spline1D(s.set.alpha_cd, s.set.cd_list) | ||
| 210 | clear!(s) | ||
| 211 | return s | ||
| 212 | end | ||
| 213 | |||
| 214 | """ | ||
| 215 | calc_particle_forces!(s::KPS4, pos1, pos2, vel1, vel2, spring, segments, d_tether, rho, i) | ||
| 216 | |||
| 217 | Calculate the drag force of the tether segment, defined by the parameters pos1, pos2, vel1 and vel2 | ||
| 218 | and distribute it equally on the two particles, that are attached to the segment. | ||
| 219 | The result is stored in the array s.forces. | ||
| 220 | """ | ||
| 221 | @inline function calc_particle_forces!(s, pos1, pos2, vel1, vel2, spring, segments, d_tether, rho, i) | ||
| 222 | l_0 = spring.length # Unstressed length | ||
| 223 | k = spring.c_spring * s.stiffness_factor # Spring constant | ||
| 224 | c = spring.damping # Damping coefficient | ||
| 225 | segment = pos1 - pos2 | ||
| 226 | rel_vel = vel1 - vel2 | ||
| 227 | av_vel = 0.5 * (vel1 + vel2) | ||
| 228 | norm1 = norm(segment) | ||
| 229 | unit_vector = segment / norm1 | ||
| 230 | |||
| 231 | k1 = 0.25 * k # compression stiffness kite segments | ||
| 232 | k2 = 0.1 * k # compression stiffness tether segments | ||
| 233 | c1 = 6.0 * c # damping kite segments | ||
| 234 | spring_vel = unit_vector â‹… rel_vel | ||
| 235 | if (norm1 - l_0) > 0.0 | ||
| 236 | if i > segments # kite springs | ||
| 237 | s.spring_force .= (k * (norm1 - l_0) + (c1 * spring_vel)) * unit_vector | ||
| 238 | else | ||
| 239 | s.spring_force .= (k * (norm1 - l_0) + (c * spring_vel)) * unit_vector | ||
| 240 | end | ||
| 241 | elseif i > segments # kite spring | ||
| 242 | s.spring_force .= (k1 * (norm1 - l_0) + (c * spring_vel)) * unit_vector | ||
| 243 | else | ||
| 244 | s.spring_force .= (k2 * (norm1 - l_0) + (c * spring_vel)) * unit_vector | ||
| 245 | end | ||
| 246 | |||
| 247 | s.v_apparent .= s.v_wind_tether - av_vel | ||
| 248 | if s.set.version == 1 | ||
| 249 | area = norm1 * d_tether | ||
| 250 | else | ||
| 251 | area = norm1 * s.set.d_line * 0.001 | ||
| 252 | end | ||
| 253 | 1 (2 %) |
1 (2 %)
samples spent in calc_particle_forces!
1 (100 %) (incl.) when called from inner_loop! line 336
1 (100 %)
samples spent calling
dot
v_app_perp = s.v_apparent - s.v_apparent â‹… unit_vector * unit_vector
|
|
| 254 | half_drag_force = (-0.25 * rho * s.set.cd_tether * norm(v_app_perp) * area) * v_app_perp | ||
| 255 | |||
| 256 | @inbounds s.forces[spring.p1] .+= half_drag_force + s.spring_force | ||
| 257 | @inbounds s.forces[spring.p2] .+= half_drag_force - s.spring_force | ||
| 258 | if i == 1 s.last_force .= s.forces[spring.p1] end | ||
| 259 | nothing | ||
| 260 | end | ||
| 261 | |||
| 262 | """ | ||
| 263 | calc_aero_forces!(s::KPS4, pos, vel, rho, alpha_depower, rel_steering) | ||
| 264 | |||
| 265 | Calculates the aerodynamic forces acting on the kite particles. | ||
| 266 | |||
| 267 | Parameters: | ||
| 268 | - pos: vector of the particle positions | ||
| 269 | - vel: vector of the particle velocities | ||
| 270 | - rho: air density [kg/m^3] | ||
| 271 | - rel_depower: value between 0.0 and 1.0 | ||
| 272 | - alpha_depower: depower angle [degrees] | ||
| 273 | - rel_steering: value between -1.0 and +1.0 | ||
| 274 | |||
| 275 | Updates the vector s.forces of the first parameter. | ||
| 276 | """ | ||
| 277 | function calc_aero_forces!(s::KPS4, pos, vel, rho, alpha_depower, rel_steering) | ||
| 278 | rel_side_area = s.set.rel_side_area/100.0 # defined in percent | ||
| 279 | K = 1 - rel_side_area # correction factor for the drag | ||
| 280 | # pos_B, pos_C, pos_D: position of the kite particles B, C, and D | ||
| 281 | # v_B, v_C, v_D: velocity of the kite particles B, C, and D | ||
| 282 | pos_B, pos_C, pos_D = pos[s.set.segments+3], pos[s.set.segments+4], pos[s.set.segments+5] | ||
| 283 | v_B, v_C, v_D = vel[s.set.segments+3], vel[s.set.segments+4], vel[s.set.segments+5] | ||
| 284 | va_2, va_3, va_4 = s.v_wind - v_B, s.v_wind - v_C, s.v_wind - v_D | ||
| 285 | |||
| 286 | pos_centre = 0.5 * (pos_C + pos_D) | ||
| 287 | delta = pos_B - pos_centre | ||
| 288 | z = -normalize(delta) | ||
| 289 | y = normalize(pos_C - pos_D) | ||
| 290 | x = y × z | ||
| 291 | s.x .= x; s.y .= y; s.z .= z # save the kite reference frame in the state | ||
| 292 | |||
| 293 | va_xz2 = va_2 - (va_2 â‹… y) * y | ||
| 294 | va_xy3 = va_3 - (va_3 â‹… z) * z | ||
| 295 | va_xy4 = va_4 - (va_4 â‹… z) * z | ||
| 296 | |||
| 297 | alpha_2 = rad2deg(Ï€ - acos2(normalize(va_xz2) â‹… x) - alpha_depower) + s.set.alpha_zero | ||
| 298 | alpha_3 = rad2deg(Ï€ - acos2(normalize(va_xy3) â‹… x) - rel_steering * KS) + s.set.alpha_ztip | ||
| 299 | alpha_4 = rad2deg(Ï€ - acos2(normalize(va_xy4) â‹… x) + rel_steering * KS) + s.set.alpha_ztip | ||
| 300 | |||
| 301 | 3 (5 %) |
3 (100 %)
samples spent calling
Spline1D
CL2, CD2 = calc_cl(alpha_2), DRAG_CORR * calc_cd(alpha_2)
|
|
| 302 | 1 (2 %) |
1 (100 %)
samples spent calling
Spline1D
CL3, CD3 = calc_cl(alpha_3), DRAG_CORR * calc_cd(alpha_3)
|
|
| 303 | 3 (5 %) |
3 (100 %)
samples spent calling
Spline1D
CL4, CD4 = calc_cl(alpha_4), DRAG_CORR * calc_cd(alpha_4)
|
|
| 304 | L2 = (-0.5 * rho * (norm(va_xz2))^2 * s.set.area * CL2) * normalize(va_2 × y) | ||
| 305 | L3 = (-0.5 * rho * (norm(va_xy3))^2 * s.set.area * rel_side_area * CL3) * normalize(va_3 × z) | ||
| 306 | L4 = (-0.5 * rho * (norm(va_xy4))^2 * s.set.area * rel_side_area * CL4) * normalize(z × va_4) | ||
| 307 | D2 = (-0.5 * K * rho * norm(va_2) * s.set.area * CD2) * va_2 | ||
| 308 | D3 = (-0.5 * K * rho * norm(va_3) * s.set.area * rel_side_area * CD3) * va_3 | ||
| 309 | D4 = (-0.5 * K * rho * norm(va_4) * s.set.area * rel_side_area * CD4) * va_4 | ||
| 310 | s.lift_force .= L2 | ||
| 311 | s.drag_force .= D2 + D3 + D4 | ||
| 312 | |||
| 313 | s.forces[s.set.segments + 3] .+= (L2 + D2) | ||
| 314 | s.forces[s.set.segments + 4] .+= (L3 + D3) | ||
| 315 | s.forces[s.set.segments + 5] .+= (L4 + D4) | ||
| 316 | end | ||
| 317 | |||
| 318 | """ | ||
| 319 | inner_loop!(s::KPS4, pos, vel, v_wind_gnd, segments, d_tether) | ||
| 320 | |||
| 321 | Calculate the forces, acting on all particles. | ||
| 322 | |||
| 323 | Output: | ||
| 324 | - s.forces | ||
| 325 | - `s.v_wind_tether` | ||
| 326 | """ | ||
| 327 | @inline function inner_loop!(s::KPS4, pos, vel, v_wind_gnd, segments, d_tether) | ||
| 328 | for i in eachindex(s.springs) | ||
| 329 | p1 = s.springs[i].p1 # First point nr. | ||
| 330 | p2 = s.springs[i].p2 # Second point nr. | ||
| 331 | height = 0.5 * (pos[p1][3] + pos[p2][3]) | ||
| 332 | rho = calc_rho(s.am, height) | ||
| 333 | @assert height > 0 | ||
| 334 | |||
| 335 | 2 (3 %) |
2 (100 %)
samples spent calling
calc_wind_factor
s.v_wind_tether .= calc_wind_factor(s.am, height) * v_wind_gnd
|
|
| 336 | 4 (7 %) |
3 (75 %)
samples spent calling
calc_particle_forces!
calc_particle_forces!(s, pos[p1], pos[p2], vel[p1], vel[p2], s.springs[i], segments, d_tether, rho, i)
1 (25 %) samples spent calling calc_particle_forces! |
|
| 337 | end | ||
| 338 | nothing | ||
| 339 | end | ||
| 340 | |||
| 341 | """ | ||
| 342 | loop!(s::KPS4, pos, vel, posd, veld) | ||
| 343 | |||
| 344 | Calculate the vectors s.res1 and calculate s.res2 using loops | ||
| 345 | that iterate over all tether segments. | ||
| 346 | """ | ||
| 347 | function loop!(s::KPS4, pos, vel, posd, veld) | ||
| 348 | L_0 = s.l_tether / s.set.segments | ||
| 349 | if s.set.version == 1 | ||
| 350 | # for compatibility with the python code and paper | ||
| 351 | mass_per_meter = 0.011 | ||
| 352 | else | ||
| 353 | mass_per_meter = s.set.rho_tether * π * (s.set.d_tether/2000.0)^2 | ||
| 354 | end | ||
| 355 | s.res1[1] .= pos[1] | ||
| 356 | s.res2[1] .= vel[1] | ||
| 357 | particles = s.set.segments + KITE_PARTICLES + 1 | ||
| 358 | for i in 2:particles | ||
| 359 | s.res1[i] .= vel[i] - posd[i] | ||
| 360 | end | ||
| 361 | # Compute the masses and forces | ||
| 362 | m_tether_particle = mass_per_meter * s.segment_length | ||
| 363 | s.masses[s.set.segments+1] = s.set.kcu_mass + 0.5 * m_tether_particle | ||
| 364 | # TODO: check if the next two lines are correct | ||
| 365 | damping = s.set.damping / L_0 | ||
| 366 | c_spring = s.set.c_spring/L_0 | ||
| 367 | for i in 1:s.set.segments | ||
| 368 | @inbounds s.masses[i] = m_tether_particle | ||
| 369 | @inbounds s.springs[i] = SP(s.springs[i].p1, s.springs[i].p2, s.segment_length, c_spring, damping) | ||
| 370 | end | ||
| 371 | 6 (10 %) | inner_loop!(s, pos, vel, s.v_wind_gnd, s.set.segments, s.set.d_tether/1000.0) | |
| 372 | for i in 2:particles | ||
| 373 | s.res2[i] .= veld[i] - (SVector(0, 0, -G_EARTH) - s.forces[i] / s.masses[i]) | ||
| 374 | end | ||
| 375 | end | ||
| 376 | |||
| 377 | """ | ||
| 378 | residual!(res, yd, y::MVector{S, SimFloat}, s::KPS4, time) where S | ||
| 379 | |||
| 380 | N-point tether model, four points for the kite on top: | ||
| 381 | Inputs: | ||
| 382 | State vector y = pos1, pos2, ... , posn, vel1, vel2, . .., veln, length, v_reel_out | ||
| 383 | Derivative yd = posd1, posd2, ..., posdn, veld1, veld2, ..., veldn, lengthd, v_reel_outd | ||
| 384 | Output: | ||
| 385 | Residual res = res1, res2 = vel1-posd1, ..., veld1-acc1, ..., | ||
| 386 | |||
| 387 | Additional parameters: | ||
| 388 | s: Struct with work variables, type KPS4 | ||
| 389 | S: The dimension of the state vector | ||
| 390 | The number of the point masses of the model N = S/6, the state of each point | ||
| 391 | is represented by two 3 element vectors. | ||
| 392 | """ | ||
| 393 |
48 (83 %)
samples spent in residual!
function residual!(res, yd, y::Vector{SimFloat}, s::KPS4, time)
28 (58 %) (ex.), 48 (100 %) (incl.) when called from DAEFunction line 2207 |
||
| 394 | S = length(y) | ||
| 395 | 13 (22 %) | 14 (24 %) |
1 (100 %)
samples spent calling
StaticArray
y_ = MVector{S, SimFloat}(y)
|
| 396 | 14 (24 %) | 15 (26 %) |
1 (100 %)
samples spent calling
StaticArray
yd_ = MVector{S, SimFloat}(yd)
|
| 397 | 1 (2 %) | 19 (33 %) |
18 (100 %)
samples spent calling
residual!
residual!(res, yd_, y_, s, time)
|
| 398 | end | ||
| 399 | function residual!(res, yd, y::MVector{S, SimFloat}, s::KPS4, time) where S | ||
| 400 | T = S-2 # T: three times the number of particles excluding the origin | ||
| 401 | segments = div(T,6) - KITE_PARTICLES | ||
| 402 | |||
| 403 | # Reset the force vector to zero. | ||
| 404 | for i in 1:segments+KITE_PARTICLES+1 | ||
| 405 | 1 (2 %) |
1 (100 %)
samples spent calling
materialize!
s.forces[i] .= SVector(0.0, 0, 0)
|
|
| 406 | end | ||
| 407 | # extract the data for the winch simulation | ||
| 408 | length, v_reel_out = y[end-1], y[end] | ||
| 409 | lengthd, v_reel_outd = yd[end-1], yd[end] | ||
| 410 | # extract the data of the particles | ||
| 411 | y_ = @view y[1:end-2] | ||
| 412 | yd_ = @view yd[1:end-2] | ||
| 413 | # unpack the vectors y and yd | ||
| 414 | part = reshape(SVector{T}(y_), Size(3, div(T,6), 2)) | ||
| 415 | partd = reshape(SVector{T}(yd_), Size(3, div(T,6), 2)) | ||
| 416 | pos1, vel1 = part[:,:,1], part[:,:,2] | ||
| 417 | pos = SVector{div(T,6)+1}(if i==1 SVector(0.0,0,0) else SVector(pos1[:,i-1]) end for i in 1:div(T,6)+1) | ||
| 418 | 1 (2 %) |
1 (100 %)
samples spent calling
StaticArray
vel = SVector{div(T,6)+1}(if i==1 SVector(0.0,0,0) else SVector(vel1[:,i-1]) end for i in 1:div(T,6)+1)
|
|
| 419 | posd1, veld1 = partd[:,:,1], partd[:,:,2] | ||
| 420 | posd = SVector{div(T,6)+1}(if i==1 SVector(0.0,0,0) else SVector(posd1[:,i-1]) end for i in 1:div(T,6)+1) | ||
| 421 | veld = SVector{div(T,6)+1}(if i==1 SVector(0.0,0,0) else SVector(veld1[:,i-1]) end for i in 1:div(T,6)+1) | ||
| 422 | @assert ! isnan(pos[2][3]) | ||
| 423 | |||
| 424 | # core calculations | ||
| 425 | s.l_tether = length | ||
| 426 | s.segment_length = length / segments | ||
| 427 | 7 (12 %) |
7 (100 %)
samples spent calling
calc_aero_forces!
calc_aero_forces!(s, pos, vel, s.rho, s.alpha_depower, s.steering)
|
|
| 428 | 6 (10 %) |
6 (100 %)
samples spent calling
loop!
loop!(s, pos, vel, posd, veld)
|
|
| 429 | |||
| 430 | # winch calculations | ||
| 431 | res[end-1] = lengthd - v_reel_out | ||
| 432 | res[end] = v_reel_outd - calc_acceleration(s.wm, s.sync_speed, v_reel_out, norm(s.forces[1]), true) | ||
| 433 | |||
| 434 | # copy and flatten result | ||
| 435 | for i in 2:div(T,6)+1 | ||
| 436 | for j in 1:3 | ||
| 437 | 2 (3 %) | @inbounds res[3*(i-2)+j] = s.res1[i][j] | |
| 438 | @inbounds res[3*(div(T,6))+3*(i-2)+j] = s.res2[i][j] | ||
| 439 | end | ||
| 440 | end | ||
| 441 | # copy the position vector for easy debugging | ||
| 442 | for i in 1:div(T,6)+1 | ||
| 443 | @inbounds s.pos[i] .= pos[i] | ||
| 444 | 1 (2 %) |
1 (100 %)
samples spent calling
iterate
end
|
|
| 445 | s.vel_kite .= vel[end-2] | ||
| 446 | s.v_reel_out = v_reel_out | ||
| 447 | |||
| 448 | @assert ! isnan(norm(res)) | ||
| 449 | s.iter += 1 | ||
| 450 | |||
| 451 | nothing | ||
| 452 | end | ||
| 453 | |||
| 454 | # =================== getter functions ==================================================== | ||
| 455 | |||
| 456 | """ | ||
| 457 | calc_height(s::KPS4) | ||
| 458 | |||
| 459 | Determine the height of the topmost kite particle above ground. | ||
| 460 | """ | ||
| 461 | function calc_height(s::KPS4) | ||
| 462 | pos_kite(s)[3] | ||
| 463 | end | ||
| 464 | |||
| 465 | """ | ||
| 466 | pos_kite(s::KPS4) | ||
| 467 | |||
| 468 | Return the position of the kite (top particle). | ||
| 469 | """ | ||
| 470 | function pos_kite(s::KPS4) | ||
| 471 | s.pos[end-2] | ||
| 472 | end | ||
| 473 | |||
| 474 | """ | ||
| 475 | kite_ref_frame(s::KPS4) | ||
| 476 | |||
| 477 | Returns a tuple of the x, y, and z vectors of the kite reference frame. | ||
| 478 | """ | ||
| 479 | function kite_ref_frame(s::KPS4) | ||
| 480 | s.x, s.y, s.z | ||
| 481 | end | ||
| 482 | |||
| 483 | """ | ||
| 484 | winch_force(s::KPS4) | ||
| 485 | |||
| 486 | Return the absolute value of the force at the winch as calculated during the last timestep. | ||
| 487 | """ | ||
| 488 | function winch_force(s::KPS4) norm(s.last_force) end | ||
| 489 | |||
| 490 | # ==================== end of getter functions ================================================ | ||
| 491 | |||
| 492 | function spring_forces(s::KPS4) | ||
| 493 | forces = zeros(SimFloat, s.set.segments+KITE_SPRINGS) | ||
| 494 | for i in 1:s.set.segments | ||
| 495 | forces[i] = s.springs[i].c_spring * (norm(s.pos[i+1] - s.pos[i]) - s.segment_length) * s.stiffness_factor | ||
| 496 | if forces[i] > 4000.0 | ||
| 497 | println("Tether raptures for segment $i !") | ||
| 498 | end | ||
| 499 | end | ||
| 500 | for i in 1:KITE_SPRINGS | ||
| 501 | p1 = s.springs[i+s.set.segments].p1 # First point nr. | ||
| 502 | p2 = s.springs[i+s.set.segments].p2 # Second point nr. | ||
| 503 | pos1, pos2 = s.pos[p1], s.pos[p2] | ||
| 504 | spring = s.springs[i+s.set.segments] | ||
| 505 | l_0 = spring.length # Unstressed length | ||
| 506 | k = spring.c_spring * s.stiffness_factor # Spring constant | ||
| 507 | segment = pos1 - pos2 | ||
| 508 | norm1 = norm(segment) | ||
| 509 | k1 = 0.25 * k # compression stiffness kite segments | ||
| 510 | if (norm1 - l_0) > 0.0 | ||
| 511 | spring_force = k * (norm1 - l_0) | ||
| 512 | else | ||
| 513 | spring_force = k1 * (norm1 - l_0) | ||
| 514 | end | ||
| 515 | forces[i+s.set.segments] = spring_force | ||
| 516 | if norm(s.spring_force) > 4000.0 | ||
| 517 | println("Bridle brakes for spring $i !") | ||
| 518 | end | ||
| 519 | end | ||
| 520 | forces | ||
| 521 | end | ||
| 522 | |||
| 523 | """ | ||
| 524 | find_steady_state!(s::KPS4; prn=false, delta = 0.0, stiffness_factor=0.035) | ||
| 525 | |||
| 526 | Find an initial equilibrium, based on the inital parameters | ||
| 527 | `l_tether`, elevation and `v_reel_out`. | ||
| 528 | """ | ||
| 529 | function find_steady_state!(s::KPS4; prn=false, delta = 0.0, stiffness_factor=0.035) | ||
| 530 | s.stiffness_factor = stiffness_factor | ||
| 531 | res = zeros(MVector{6*(s.set.segments+KITE_PARTICLES)+2, SimFloat}) | ||
| 532 | iter = 0 | ||
| 533 | |||
| 534 | # helper function for the steady state finder | ||
| 535 | function test_initial_condition!(F, x::Vector) | ||
| 536 | x1 = copy(x) | ||
| 537 | y0, yd0 = init(s, x1) | ||
| 538 | residual!(res, yd0, y0, s, 0.0) | ||
| 539 | for i in 1:s.set.segments+KITE_PARTICLES-1 | ||
| 540 | if i != s.set.segments+KITE_PARTICLES-1 | ||
| 541 | j = i | ||
| 542 | else | ||
| 543 | j = i + 1 | ||
| 544 | end | ||
| 545 | # copy the x-component of the residual res2 (acceleration) | ||
| 546 | F[i] = res[1 + 3*(j-1) + 3*(s.set.segments+KITE_PARTICLES)] | ||
| 547 | # copy the z-component of the residual res2 | ||
| 548 | F[i+s.set.segments+KITE_PARTICLES] = res[3 + 3*(j-1) + 3*(s.set.segments+KITE_PARTICLES)] | ||
| 549 | end | ||
| 550 | # copy the acceleration of point KCU in x direction | ||
| 551 | i = s.set.segments+1 | ||
| 552 | F[end-1] = res[1 + 3*(i-1) + 3*(s.set.segments+KITE_PARTICLES)] | ||
| 553 | # copy the acceleration of point C in y direction | ||
| 554 | i = s.set.segments+3 | ||
| 555 | F[end] = res[2 + 3*(i-1) + 3*(s.set.segments+KITE_PARTICLES)] | ||
| 556 | iter += 1 | ||
| 557 | return nothing | ||
| 558 | end | ||
| 559 | if prn println("\nStarted function test_nlsolve...") end | ||
| 560 | X00 = zeros(SimFloat, 2*(s.set.segments+KITE_PARTICLES-1)+2) | ||
| 561 | results = nlsolve(test_initial_condition!, X00, autoscale=true, xtol=2e-7, ftol=2e-7, iterations=s.set.max_iter) | ||
| 562 | if prn println("\nresult: $results") end | ||
| 563 | init(s, results.zero) | ||
| 564 | end |