diff --git a/.gitignore b/.gitignore index 55a910501..8e989afc8 100644 --- a/.gitignore +++ b/.gitignore @@ -48,3 +48,5 @@ data/settings_v9f.yaml Manifest-v1.10.toml.bak Manifest-v1.11.toml.bak Manifest-v1.11.toml.backup +data/prob_dynamic_1.10_3_seg.bin.default +data/prob_dynamic_1.11_3_seg.bin.default diff --git a/Project.toml b/Project.toml index 44aa56e4b..dc8a8eb42 100644 --- a/Project.toml +++ b/Project.toml @@ -60,7 +60,7 @@ DocStringExtensions = "0.9.4" Documenter = "1.10.1" Interpolations = "0.15.1" KitePodModels = "0.3.8" -KiteUtils = "0.10" +KiteUtils = "0.10.3" LaTeXStrings = "1.4.0" LinearSolve = "~2.39.0" ModelingToolkit = "~9.71.0" diff --git a/data/settings.yaml b/data/settings.yaml index a4b00b7f7..5aee4f242 100644 --- a/data/settings.yaml +++ b/data/settings.yaml @@ -60,18 +60,6 @@ kps4: smc: 0.0 # steering moment coefficient [-] cmq: 0.0 # pitch rate dependant moment coefficient [-] cord_length: 2.0 # average aerodynamic cord length of the kite [m] - -kps4_3l: - radius: 2.0 # the radius of the circle shape on which the kite lines, viewed - # from the front [m] - bridle_center_distance: 4.0 # the distance from point the center bridle connection point of - # the middle line to the kite [m] - middle_length: 1.5 # the cord length of the kite in the middle [m] - tip_length: 0.62 # the cord length of the kite at the tips [m] - min_steering_line_distance: 1.0 # the distance between the left and right steering bridle [m] - # line connections on the kite that are closest to each other [m] - width_3l: 4.1 # width of the kite [m] - aero_surfaces: 3 # the number of aerodynamic surfaces to use per mass point [-] bridle: d_line: 2.5 # bridle line diameter [mm] diff --git a/data/settings_kps5.yaml b/data/settings_kps5.yaml new file mode 100644 index 000000000..5d2e6c300 --- /dev/null +++ b/data/settings_kps5.yaml @@ -0,0 +1,142 @@ +system: + log_file: "data/log_8700W_8ms" # filename without extension [replay only] + # use / as path delimiter, even on Windows + log_level: 2 # 0: no logging + time_lapse: 1.0 # relative replay speed + sim_time: 5 # simulation time [sim only] + segments: 6 # number of tether segments + sample_freq: 20 # sample frequency in Hz + zoom: 0.03 # zoom factor for the system view + kite_scale: 3.0 # relative zoom factor for the 4 point kite + fixed_font: "" # name or filepath+filename of alternative fixed pitch font, e.g. Liberation Mono + +initial: + l_tether: 10 # initial tether length [m] + elevation: 45 # initial elevation angle [deg] + v_reel_out: 0.5 # initial reel out speed [m/s] + depower: 25.0 # initial depower settings [%] + +solver: + abs_tol: 0.0006 # absolute tolerance of the DAE solver [m, m/s] + rel_tol: 0.001 # relative tolerance of the DAE solver [-] + solver: "DFBDF" # DAE solver, IDA or DFBDF or DImplicitEuler + linear_solver: "GMRES" # can be GMRES or LapackDense or Dense (only for IDA) + max_order: 4 # maximal order, usually between 3 and 5 (IDA and DFBDF) + max_iter: 200 # max number of iterations of the steady-state-solver + +steering: + c0: 0.0 # steering offset -0.0032 [-] + c_s: 2.59 # steering coefficient one point model; 2.59 was 0.6; TODO: check if it must be divided by kite_area + c2_cor: 0.93 # correction factor one point model + k_ds: 1.5 # influence of the depower angle on the steering sensitivity + delta_st: 0.02 # steering increment (when pressing RIGHT) + max_steering: 16.834 # max. steering angle of the side planes for four point model [degrees] + cs_4p: 1.0 # correction factor for the steering coefficient of the four point model + +depower: + alpha_d_max: 31.0 # max depower angle [deg] + depower_offset: 23.6 # at rel_depower=0.236 the kite is fully powered [%] + +kite: + model: "data/kite.obj" # 3D model of the kite + physical_model: "KPS4" # name of the kite model to use (KPS3 or KPS4) + version: 1 # version of the model to use + mass: 6.2 # kite mass incl. sensor unit [kg] + area: 10.18 # projected kite area [m²] + rel_side_area: 30.6 # relative side area [%] + height_k: 2.23 # height of the kite [m] + alpha_cl: [-180.0, -160.0, -90.0, -20.0, -10.0, -5.0, 0.0, 20.0, 40.0, 90.0, 160.0, 180.0] + cl_list: [ 0.0, 0.5, 0.0, 0.08, 0.125, 0.15, 0.2, 1.0, 1.0, 0.0, -0.5, 0.0] + alpha_cd: [-180.0, -170.0, -140.0, -90.0, -20.0, 0.0, 20.0, 90.0, 140.0, 170.0, 180.0] + cd_list: [ 0.5, 0.5, 0.5, 1.0, 0.2, 0.1, 0.2, 1.0, 0.5, 0.5, 0.5] + +kps4: + width: 5.77 # width of the kite [m] + alpha_zero: 4.0 # should be 4 .. 10 [degrees] + alpha_ztip: 10.0 # [degrees] + m_k: 0.2 # relative nose distance; increasing m_k increases C2 of the turn-rate law + rel_nose_mass: 0.47 # relative nose mass + rel_top_mass: 0.4 # mass of the top particle relative to the sum of top and side particles + smc: 0.0 # steering moment coefficient [-] + cmq: 0.0 # pitch rate dependant moment coefficient [-] + cord_length: 2.0 # average aerodynamic cord length of the kite [m] + +kps5: + c_spring_kite: 614600.0 # unit spring constant coefficient of the kite springs [N] + damping_kite_springs: 473.0 # unit damping coefficient [Ns] + rel_mass_p2: 0.25 # rel mass of p2 + rel_mass_p3: 0.25 # rel mass of p3 + rel_mass_p4: 0.25 # rel mass of p4 and p5 + +bridle: + d_line: 2.5 # bridle line diameter [mm] + l_bridle: 33.4 # sum of the lengths of the bridle lines [m] + h_bridle: 4.9 # height of bridle [m] + rel_compr_stiffness: 0.25 # relative compression stiffness of the kite springs [-] + rel_damping: 6.0 # relative damping of the kite spring (relative to main tether) [-] + +kcu: + kcu_model: "KCU1" # name of the kite control unit model, KCU1 or KCU2 + kcu_mass: 8.4 # mass of the kite control unit [kg] + kcu_diameter: 0.4 # diameter of the KCU for drag calculation [m] + cd_kcu: 0.47 #nonzero, changed form 0.0 tp 0.47 # drag coefficient of the KCU [-] + power2steer_dist: 1.3 # [m] + depower_drum_diameter: 0.069 # [m] + tape_thickness: 0.0006 # [m] + v_depower: 0.075 # max velocity of depowering in units per second (full range: 1 unit) + v_steering: 0.2 # max velocity of steering in units per second (full range: 2 units) + depower_gain: 3.0 # 3.0 means: more than 33% error -> full speed + steering_gain: 3.0 + +tether: + d_tether: 4 # tether diameter [mm] + cd_tether: 0.958 # drag coefficient of the tether + damping: 473.0 # unit damping coefficient [Ns] + c_spring: 614600.0 # unit spring constant coefficient [N] + rho_tether: 724.0 # density of Dyneema [kg/m³] + e_tether: 55000000000.0 # axial tensile modulus of Dyneema (M.B. Ruppert) [Pa] + # SK75: 109 to 132 GPa according to datasheet + +winch: + winch_model: "AsyncMachine" # or TorqueControlledMachine + max_force: 4000 # maximal (nominal) tether force; short overload allowed [N] + v_ro_max: 8.0 # maximal reel-out speed [m/s] + v_ro_min: -8.0 # minimal reel-out speed (=max reel-in speed) [m/s] + drum_radius: 0.1615 # radius of the drum [m] + max_acc: 4.0 # maximal acceleration of the winch [m/s²] + gear_ratio: 6.2 # gear ratio of the winch [-] + inertia_total: 0.204 # total inertia, as seen from the motor/generator [kgm²] + f_coulomb: 122.0 # coulomb friction [N] + c_vf: 30.6 # coefficient for the viscous friction [Ns/m] + p_speed: 10000.0 # proportional gain of the winch speed controller [-] + i_speed: 2500.0 # integral gain of the winch speed controller [-] + +environment: + v_wind: 9.5 # wind speed at reference height [m/s] + upwind_dir: -180 # upwind direction (0 is form North, clockwise positive) [deg] + temp_ref: 15.0 # temperature at reference height [°C] + height_gnd: 0.0 # height of groundstation above see level [m] + h_ref: 6.0 # reference height for the wind speed [m] + + rho_0: 1.225 # air density at zero height and 15 °C [kg/m³] + alpha: 0.08163 # exponent of the wind profile law + z0: 0.0002 # surface roughness [m] + profile_law: 3 # 1=EXP, 2=LOG, 3=EXPLOG, 4=FAST_EXP, 5=FAST_LOG, 6=FAST_EXPLOG + # the following parameters are for calculating the turbulent wind field using the Mann model + use_turbulence: 0.0 # turbulence intensity relative to Cabauw, NL + v_wind_gnds: [3.483, 5.324, 8.163] # wind speeds at ref height for calculating the turbulent wind field [m/s] + avg_height: 200.0 # average height during reel out [m] + rel_turbs: [0.342, 0.465, 0.583] # relative turbulence at the v_wind_gnds + i_ref: 0.14 # is the expected value of the turbulence intensity at 15 m/s. + v_ref: 42.9 # five times the average wind speed in m/s at hub height over the full year [m/s] + # Cabauw: 8.5863 m/s * 5.0 = 42.9 m/s + height_step: 2.0 # use a grid with 2m resolution in z direction [m] + grid_step: 2.0 # grid resolution in x and y direction [m] + + + + + # needed: + # springcomstamt_bridle + # springcomstamt_tether + # rel_damping_kite: \ No newline at end of file diff --git a/data/system_kps5.yaml b/data/system_kps5.yaml new file mode 100644 index 000000000..171ec6e75 --- /dev/null +++ b/data/system_kps5.yaml @@ -0,0 +1,3 @@ +system: + sim_settings: "settings_kps5.yaml" # simulator settings + \ No newline at end of file diff --git a/examples/KPS5onecode.jl b/examples/KPS5onecode.jl new file mode 100644 index 000000000..bf778add5 --- /dev/null +++ b/examples/KPS5onecode.jl @@ -0,0 +1,453 @@ +using ModelingToolkit: t_nounits as t, D_nounits as D +using KiteUtils +using Timers +using Pkg +#using DAEProblemLibrary +# if ! ("ControlPlots" ∈ keys(Pkg.project().dependencies)) +# using TestEnv; TestEnv.activate() +# end +using ControlPlots +import OrdinaryDiffEqCore.init +import OrdinaryDiffEqCore.step! +import OrdinaryDiffEqCore.solve +import OrdinaryDiffEqCore +using Dierckx, ModelingToolkit, LinearAlgebra, Statistics, Parameters, + OrdinaryDiffEqCore, OrdinaryDiffEqBDF#, OrdinaryDiffEqSDIRK, NonlinearSolve, FiniteDiff, DifferentiationInterface +toc() +@with_kw mutable struct Settings2 @deftype Float64 + g_earth::Vector{Float64} = [0.0, 0.0, -9.81] # gravitational acceleration [m/s²] + v_wind_tether::Vector{Float64} = [13.5, 0.0, 0.0] # wind velocity [m/s] + rho::Float64 = 1.225 # air density [kg/m³] + sim_time::Float64 = 5 # simulation duration [s] + dt = 0.05 # time step [s] + tol = 0.0006 # tolerance for the solver + save::Bool = false # save animation frames + # kite + mass::Float64 = 6.2 # mass of kite [kg] + Area::Float64 = 10.18 # surface area [m²] + width::Float64 = 5.77 # width of kite [m] + height::Float64 = 2.23 # height of kite [m] + chord_length::Float64 = 2.0 # chord length [m] + # KCU + bridle + h_bridle = 4.9 # height of bridle [m] + d_line::Float64 = 2.5 # bridle line diameter [mm] + l_bridle::Float64 = 33.4 # sum of the lengths of the bridle lines [m] + kcu_cd::Float64 = 0.47 # KCU drag coefficient + kcu_diameter::Float64 = 0.4 # KCU diameter [m] + kcu_mass::Float64 = 8.4 # mass of KCU (at bridle point)[kg] + # spring constants + springconstant_tether::Float64 = 614600.0 # TETHER unit spring constant [N] + springconstant_bridle::Float64 = 614600.0 # BRIDLE unit spring constant [N] + springconstant_kite::Float64 = 614600.0 # KITE unit spring constant [N] + # damping + damping_tether::Float64 = 473 # TETHER unit damping coefficient [Ns] + rel_damping_kite::Float64 = 6.0 # KITE relative unit damping coefficient [-] + rel_damping_bridle:: Float64 = 6.0 # BRIDLE relative unit damping coefficient [-] + # tether + l_tether::Float64 = 10.0 # tether length [m] + v_reel_out ::Float64 = 0.5 # reel-out speed [m/s] + rho_tether::Float64 = 724.0 # density of tether [kg/m³] + cd_tether::Float64 = 0.958 # drag coefficient of tether + d_tether::Float64 = 4 # tether diameter [mm] + elevation::Float64 = deg2rad(70.8) # Elevation angle, angle XZ plane, between origin and bridle point [rad] + segments::Int64 = 6 # number of tether segments [-] +end +@with_kw mutable struct KPS5 + "Reference to the settings2 struct" + set::Settings2 = Settings2() + sys::Union{ModelingToolkit.ODESystem, Nothing} = nothing + t_0::Float64 = 0.0 + iter::Int64 = 0 + prob::Union{OrdinaryDiffEqCore.ODEProblem, Nothing} = nothing + integrator::Union{OrdinaryDiffEqCore.ODEIntegrator, Nothing} = nothing + get_state::Function = () -> nothing +end +# ------------------------------------------------------------------------------------------------------------------------------------------------------------ +# deriving a constants for points, total segments, and connections +# --------- +function points(s) + return 5 + s.set.segments +end +function total_segments(s) + return 9 + s.set.segments +end +function getconnections(s) + conn = [(1,2), (2,3), (3,1), (1,4), (2,4), (3,4), (1,5), (2,5), (3,5)] # connections between KITE points + conn = vcat(conn, [(6+i, 6+i+1) for i in 0:(s.set.segments-2)]...) # connection between tether points + conn = vcat(conn, [(6+s.set.segments-1, 1)]) # connection final tether point to bridle point + return conn +end +# ------------------------------------------------------------------------------------------------------------------------------------------------------------------------ +# Interpolating polars using Dierckx +# ------------------------------------ +alpha_cl = [-180.0, -160.0, -90.0, -20.0, -10.0, -5.0, 0.0, 20.0, 40.0, 90.0, 160.0, 180.0] +cl_list = [ 0.0, 0.5, 0.0, 0.08, 0.125, 0.15, 0.2, 1.0, 1.0, 0.0, -0.5, 0.0] +alpha_cd = [-180.0, -170.0, -140.0, -90.0, -20.0, 0.0, 20.0, 90.0, 140.0, 170.0, 180.0] +cd_list = [ 0.5, 0.5, 0.5, 1.0, 0.2, 0.1, 0.2, 1.0, 0.5, 0.5, 0.5] +function cl_interp(alpha) + cl_spline = Spline1D(alpha_cl, cl_list) + return cl_spline(alpha) +end +function cd_interp(alpha) + cd_spline = Spline1D(alpha_cd, cd_list) + return cd_spline(alpha) +end +@register_symbolic cl_interp(alpha) +@register_symbolic cd_interp(alpha) +# ----------------------------------------------------------------------------------------------------------------------------------------------------------------- +# Initialize the simulation +# ----------------------------------------- +function init_sim!(s::KPS5) + pos, vel = calc_initial_state(s) + simple_sys, pos, vel, e_x, e_y, e_z, v_app_point, alpha1p = model(s, pos, vel) + s.sys = simple_sys + tspan = (0.0, s.set.sim_time) + s.prob = ODEProblem(simple_sys, nothing, tspan) + s.integrator = OrdinaryDiffEqCore.init(s.prob, FBDF(autodiff=false); s.set.dt, abstol=s.set.tol, save_on=false) +end +# ------------------------------ +# Calculate Initial State +# ------------------------------ +function calc_initial_state(s) + p1location = [s.set.l_tether*cos(s.set.elevation) 0 s.set.l_tether*sin(s.set.elevation)] + kitepos0rot = get_kite_points(s) + POS0 = kitepos0rot .+ p1location' + POS0 = hcat(POS0, zeros(3, 1)) + if s.set.segments > 1 + extra_nodes = [POS0[:,6] + (POS0[:,1] - POS0[:,6]) * i / s.set.segments for i in 1:(s.set.segments-1)] + POS0 = hcat(POS0, extra_nodes...) + end + VEL0 = zeros(3, points(s)) + return POS0, VEL0 +end +# ----------------------------------------------------- +# initializing kite points +# ----------------------------------------------------- +function get_kite_points(s) + # Original kite points in local reference frame + kitepos0 = # KITE points + # P1 Bridle P2 P3 P4 P5 + [0.000 s.set.chord_length/2 -s.set.chord_length/2 0 0; + 0.000 0 0 -s.set.width/2 s.set.width/2; + 0.000 s.set.height+s.set.h_bridle s.set.height+s.set.h_bridle s.set.h_bridle s.set.h_bridle] + + beta = s.set.elevation + Y_r = [sin(beta) 0 cos(beta); + 0 1 0; + -cos(beta) 0 sin(beta)] + # Apply rotation to all points + kitepos0rot = Y_r * kitepos0 + + return kitepos0rot +end +# ------------------------------------------------------------------------------------------------------------------------------------------------------------------------- +# Define a function to create reference frame update equations +# --------------------------------------------------------------- +function create_reference_frame_equations(pos, e_x, e_y, e_z) + # Calculate vectors for the reference frame + X = pos[:, 2] - pos[:, 3] # Vector from P3 to P2 + Y = pos[:, 5] - pos[:, 4] # Vector from P4 to P5 + Z = cross(X, Y) # Cross product for Z axis + # Normalize these vectors to get unit vectors + norm_X = sqrt(sum(X .^ 2)) + norm_Y = sqrt(sum(Y .^ 2)) + norm_Z = sqrt(sum(Z .^ 2)) + # Create equations to update the reference frame + ref_frame_eqs = [ + e_x[1] ~ X[1] / norm_X, + e_x[2] ~ X[2] / norm_X, + e_x[3] ~ X[3] / norm_X, + e_y[1] ~ Y[1] / norm_Y, + e_y[2] ~ Y[2] / norm_Y, + e_y[3] ~ Y[3] / norm_Y, + e_z[1] ~ Z[1] / norm_Z, + e_z[2] ~ Z[2] / norm_Z, + e_z[3] ~ Z[3] / norm_Z + ] + return ref_frame_eqs +end +# -------------------------------------------------------------------------- +# computing angle of attack +# -------------------------------------------------------------------------- +function compute_alpha1p(v_a, e_z, e_x) + # Calculate the angle Alpha1p + v_a_z = v_a ⋅ e_z # Use the symbolic dot product + v_a_x = v_a ⋅ e_x # Use the symbolic dot product + + # Use atan for the symbolic computation + alpha1p =rad2deg.(atan(v_a_z, v_a_x)) + return alpha1p +end +# ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- +# Calculate Rest Lengths +# --------------------------- +function calc_rest_lengths(s) + conn = getconnections(s) + POS0, VEL0 = calc_initial_state(s) + lengths = [norm(POS0[:,conn[i][2]] - POS0[:,conn[i][1]]) for i in 1:9] + l10 = norm(POS0[:,1] - POS0[:,6]) + lengths = vcat(lengths, [(l10 + s.set.v_reel_out*t)/s.set.segments for _ in 1:s.set.segments]...) + return lengths, l10 +end +# ------------------------------------------------------------------------------------------------------------------------------------------------------------------ +#Define the Model +# ----------------------------------------------- +function model(s, pos, vel) + POS0, VEL0 = pos, vel + rest_lengths, l_tether = calc_rest_lengths(s) + @parameters K1=s.set.springconstant_tether K2=s.set.springconstant_bridle K3=s.set.springconstant_kite C1=s.set.damping_tether C2=s.set.rel_damping_bridle*s.set.damping_tether C3=s.set.rel_damping_kite*s.set.damping_tether + @parameters m_kite=s.set.mass kcu_mass=s.set.kcu_mass rho_tether=s.set.rho_tether + @parameters rho=s.set.rho cd_tether=s.set.cd_tether d_tether=s.set.d_tether S=s.set.Area + @parameters kcu_cd=s.set.kcu_cd kcu_diameter=s.set.kcu_diameter + @variables pos(t)[1:3, 1:points(s)] = POS0 + @variables vel(t)[1:3, 1:points(s)] = VEL0 + @variables acc(t)[1:3, 1:points(s)] + @variables v_app_point(t)[1:3, 1:points(s)] + @variables segment(t)[1:3, 1:total_segments(s)] + @variables unit_vector(t)[1:3, 1:total_segments(s)] + @variables norm1(t)[1:total_segments(s)] + @variables rel_vel(t)[1:3, 1:total_segments(s)] + @variables spring_vel(t)[1:total_segments(s)] + @variables k_spring(t)[1:total_segments(s)] + @variables spring_force(t)[1:3, 1:total_segments(s)] + @variables v_apparent(t)[1:3, 1:total_segments(s)] + @variables v_app_perp(t)[1:3, 1:total_segments(s)] + @variables norm_v_app(t)[1:total_segments(s)] + @variables half_drag_force(t)[1:3, 1:total_segments(s)] + @variables drag_force(t)[1:3, 1:total_segments(s)] + @variables total_force(t)[1:3, 1:points(s)] + # local kite reference frame + @variables e_x(t)[1:3] + @variables e_y(t)[1:3] + @variables e_z(t)[1:3] + @variables alpha1p(t)[1:1] + + eqs1 = vcat(D.(pos) .~ vel, + D.(vel) .~ acc) + eqs2 = vcat(eqs1...) + eqs2 = vcat(eqs2, acc[:,6] .~ [0.0, 0.0, 0.0]) # origin is six, make 6 not being hardcoded + # ----------------------------- + # defining the connections and their respective rest lengths, unit spring constants, damping and masses + # ----------------------------- + # connections adding segment connections, from origin to bridle + conn = getconnections(s) + # final connection last tether point to bridle point + # unit spring constants (K1 tether, K2 bridle, K3 kite) + k_segments = [K2, K3, K2, K2, K3, K3, K2, K3, K3] + k_segments = vcat(k_segments, [K1 for _ in 1:s.set.segments]...) + # unit damping constants (C1 tether, C2 bridle, C3 kite) + c_segments = [C2, C3, C2, C2, C3, C3, C2, C3, C3] + c_segments = vcat(c_segments, [C1 for _ in 1:s.set.segments]...) + # masses + mass_bridlelines = ((s.set.d_line/2000)^2)*pi*rho_tether*s.set.l_bridle #total mass entire bridle + mass_halfbridleline = mass_bridlelines/8 # half the connection of bridle line to kite (to assign to each kitepoint) so the other 4 halves get assigned to bridlepoint + mass_tether = ((d_tether/2000)^2)*pi*rho_tether*l_tether + mass_tetherpoints = mass_tether/(s.set.segments+1) + mass_bridlepoint = 4*mass_halfbridleline + kcu_mass + mass_tetherpoints # 4 bridle connections, kcu and tether + m_kitepoints = (m_kite/4) + mass_halfbridleline + PointMasses = [mass_bridlepoint, m_kitepoints, m_kitepoints, m_kitepoints, m_kitepoints] + PointMasses = vcat(PointMasses, [mass_tetherpoints for _ in 1:s.set.segments]...) + # ----------------------------- + # Equations for Each Segment (Spring Forces, Drag, etc.) + # ----------------------------- + for i in 1:total_segments(s) + eqs = [ + segment[:, i] ~ pos[:, conn[i][2]] - pos[:, conn[i][1]], + norm1[i] ~ norm(segment[:, i]), + unit_vector[:, i] ~ -segment[:, i] / norm1[i], + rel_vel[:, i] ~ vel[:, conn[i][2]] - vel[:, conn[i][1]], + spring_vel[i] ~ -unit_vector[:, i] ⋅ rel_vel[:, i], + k_spring[i] ~ (k_segments[i]/rest_lengths[i]) * (0.1 + 0.9*(norm1[i] > rest_lengths[i])), #define relative compresions stiffnes + spring_force[:, i] ~ (k_spring[i]*(norm1[i] - rest_lengths[i]) + c_segments[i] * spring_vel[i]) * unit_vector[:, i], + v_apparent[:, i] ~ s.set.v_wind_tether .- (vel[:, conn[i][1]] + vel[:, conn[i][2]]) / 2, + v_app_perp[:, i] ~ v_apparent[:, i] - (v_apparent[:, i] ⋅ unit_vector[:, i]) .* unit_vector[:, i], + norm_v_app[i] ~ norm(v_app_perp[:, i]) + ] + if i > 9 # tether segments + push!(eqs, half_drag_force[:, i] ~ 0.25 * rho * cd_tether * norm_v_app[i] * (rest_lengths[i]*d_tether/1000) * v_app_perp[:, i]) + elseif i in [1, 3, 4, 7] # bridle lines, try to find Cd_bridlelines later + push!(eqs, half_drag_force[:, i] ~ 0.25 * rho * cd_tether * norm_v_app[i] * (rest_lengths[i]*(s.set.d_line/1000)) * v_app_perp[:, i]) + else # kite + push!(eqs, half_drag_force[:, i] ~ zeros(3)) + end + eqs2 = vcat(eqs2, reduce(vcat, eqs)) + end + # ----------------------------- + # Reference Frame and Aerodynamic Coefficients + # ----------------------------- + ref_frame_eqs = create_reference_frame_equations(pos, e_x, e_y, e_z) + eqs2 = vcat(eqs2, ref_frame_eqs) + # only 1 AOA + v_a_kite = s.set.v_wind_tether - (vel[:, 2] + vel[:, 3])/2 # computing AOA only for center chord # Appas.set.t wind velocity + alpha1p = compute_alpha1p(v_a_kite, e_z, e_x) # Calculate Alpha1p at this time step + eqs2 = vcat(eqs2, alpha1p[1] ~ alpha1p) # Add the equation for Alpha1p for each of 4 kite points (first bering bridle so i-1) + # getting Cl and Cd + Cl = cl_interp(alpha1p) + Cd = cd_interp(alpha1p) + + # ----------------------------- + # Force Balance at Each Point + # ----------------------------- + for i in 1:points(s) + eqs = [] + force = sum([spring_force[:, j] for j in 1:total_segments(s) if conn[j][2] == i]; init=zeros(3)) - + sum([spring_force[:, j] for j in 1:total_segments(s) if conn[j][1] == i]; init=zeros(3)) + + sum([half_drag_force[:, j] for j in 1:total_segments(s) if conn[j][1] == i]; init=zeros(3)) + + sum([half_drag_force[:, j] for j in 1:total_segments(s) if conn[j][2] == i]; init=zeros(3)) + v_app_point[:, i] ~ s.set.v_wind_tether - vel[:, i] + if i == 1 # KCU drag at bridle point + area_kcu = pi * ((kcu_diameter / 2) ^ 2) + Dx_kcu = 0.5*rho*kcu_cd *area_kcu*(v_app_point[1, i]*v_app_point[1, i]) + Dy_kcu = 0.5*rho*kcu_cd *area_kcu*(v_app_point[2, i]*v_app_point[2, i]) + Dz_kcu = 0.5*rho*kcu_cd *area_kcu*(v_app_point[3, i]*v_app_point[3, i]) + D = [Dx_kcu, Dy_kcu, Dz_kcu] + push!(eqs, total_force[:, i] ~ force + D) + elseif i in 2:5 # the kite points that get Aero Forces + + v_app_mag_squared = v_app_point[1, i]^2 + v_app_point[2, i]^2 + v_app_point[3, i]^2 + # Lift calculation + L_perpoint = (1/4) * 0.5 * rho * Cl * S * (v_app_mag_squared) + # Cross product and normalization + cross_vapp_X_e_y = cross(v_app_point[:, i], e_y) + normcross_vapp_X_e_y = norm(cross_vapp_X_e_y) + L_direction = cross_vapp_X_e_y / normcross_vapp_X_e_y + # Final lift force vector + L = L_perpoint * L_direction + + # Drag calculation + D_perpoint = (1/4) * 0.5 * rho * Cd * S * v_app_mag_squared + # Create drag direction components + D_direction = [v_app_point[1, i] / norm(v_app_point[:, i]), v_app_point[2, i] / norm(v_app_point[:, i]), v_app_point[3, i] / norm(v_app_point[:, i])] + # Final drag force vector components + D = D_perpoint * D_direction + + # Total aerodynamic force + Fa = [L[1]+ D[1], L[2]+ D[2], L[3]+ D[3]] + push!(eqs, total_force[1, i] ~ force[1] + Fa[1]) + push!(eqs, total_force[2, i] ~ force[2] + Fa[2]) + push!(eqs, total_force[3, i] ~ force[3] + Fa[3]) + elseif i != 6 + push!(eqs, total_force[:, i] ~ force) + end + push!(eqs, acc[:, i] ~ s.set.g_earth + total_force[:, i] / PointMasses[i]) + eqs2 = vcat(eqs2, reduce(vcat, eqs)) + eqs2 = vcat(eqs2, v_app_point[:, i] ~ s.set.v_wind_tether - vel[:, i]) + end + @named sys = ODESystem(reduce(vcat, Symbolics.scalarize.(eqs2)), t) + simple_sys = structural_simplify(sys) + simple_sys, pos, vel, e_x, e_y, e_z, v_app_point, alpha1p +end +# ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- +# next step function +# ----------------------------- +function next_step!(s::KPS5; dt=s.set.dt) + s.t_0 = s.integrator.t + steptime = @elapsed OrdinaryDiffEqCore.step!(s.integrator, dt, true) + s.iter += 1 + s.integrator.t, steptime +end +# ----------------------------- +# Simulation Function +# ----------------------------- +function simulate(s, logger) + dt = s.set.dt + tol = s.set.tol + tspan = (0.0, dt) + time_range = 0:dt:s.set.sim_time-dt + steps = length(time_range) + iter = 0 + for i in 1:steps + next_step!(s; dt=s.set.dt) + u = s.get_state(s.integrator) + x = u[1][1, :] + y = u[1][2, :] + z = u[1][3, :] + iter += s.iter + sys_state = SysState{points(s)}() + sys_state.X .= x + sys_state.Y .= y + sys_state.Z .= z + println("iter: $iter", " steps: $steps") + log!(logger, sys_state) + println(x[end], ", ", sys_state.X[end]) + end + println("iter: $iter", " steps: $steps") + return nothing +end +# ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- +# Generate Getters +# ----------------------------- +function generate_getters!(s) + sys = s.sys + c = collect + get_state = ModelingToolkit.getu(sys, + [c(sys.pos)] + ) + s.get_state = (integ) -> get_state(integ) + return nothing +end +function play(s, lg) + conn = getconnections(s) + sl = lg.syslog + total_segmentsvector = Vector{Int64}[] + for conn_pair in conn + push!(total_segmentsvector, Int64[conn_pair[1], conn_pair[2]]) + end + # Add tether segments + for i in 0:(s.set.segments-2) + push!(total_segmentsvector, [6+i, 6+i+1]) + end + # Add final connection from last tether point to bridle point + push!(total_segmentsvector, [6+s.set.segments-1, 1]) + for step in 1:length(0:s.set.dt:s.set.sim_time)-1 #-s.set.dt + # Get positions at this time step + x = sl.X[step] + y = sl.Y[step] + z = sl.Z[step] + # Create points array for all points in the system + pointsvector = Vector{Float64}[] + for i in 1:points(s) #FIX THIS! + push!(pointsvector, Float64[x[i], y[i], z[i]]) + end + # Calculate appropriate limits for the plot + x_min, x_max = 0, 20 + z_min, z_max = 0, 20 + t = s.set.dt * (step-1) + # Plot the kite system at this time step + plot2d(pointsvector, total_segmentsvector, t; + zoom = false, + xlim = (x_min, x_max), + ylim = (z_min, z_max) + ) + # Add a small delay to control animation speed + sleep(0.05) + end + nothing +end +function plot_front_view3(lg) + display(plotxy(lg.y, lg.z; + xlabel="pos_y [m]", + ylabel="height [m]", + fig="front_view")) + nothing +end +# ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- +# Running the Main Function +# ----------------------------- +function main() + global lg, s + set = Settings2() + s = KPS5(set=set) + time_range = 0:set.dt:set.sim_time-set.dt + steps = length(time_range) + logger = Logger(points(s), steps) + init_sim!(s) + generate_getters!(s) + simulate(s, logger) + save_log(logger, "tmp") + lg = load_log("tmp") + play(s, lg) +end + +main() +# ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- \ No newline at end of file diff --git a/examples/kps5_example.jl b/examples/kps5_example.jl new file mode 100644 index 000000000..feba1d1e0 --- /dev/null +++ b/examples/kps5_example.jl @@ -0,0 +1,167 @@ +using KiteModels +using Timers +using Pkg +if ! ("ControlPlots" ∈ keys(Pkg.project().dependencies)) + using TestEnv; TestEnv.activate() +end +using ControlPlots +import OrdinaryDiffEqCore.init +import OrdinaryDiffEqCore.step! +import OrdinaryDiffEqCore.solve +import OrdinaryDiffEqCore +using Dierckx, Interpolations, Serialization, StaticArrays, LinearAlgebra, Statistics, Parameters, NLsolve, + DocStringExtensions, OrdinaryDiffEqCore, OrdinaryDiffEqBDF, OrdinaryDiffEqSDIRK, NonlinearSolve, FiniteDiff, DifferentiationInterface +# using ModelingToolkit: Symbolics, @register_symbolic +# using ModelingToolkit, ControlPlots, Parameters +set = deepcopy(load_settings("system_kps5.yaml")) +kcu::KCU = KCU(set) +s = KPS5(kcu) + +dt = 1/s.set.sample_freq +time_range = 0:dt:set.sim_time-dt +steps = length(time_range) +logger = Logger(KiteModels.points(s), steps) +#solver = DImplicitEuler(autodiff=false)#@AutoFiniteDiff()) +KiteModels.get_kite_points(s) +KiteModels.calc_initial_state(s) +KiteModels.init_sim!(s) #(integrate gen getter in initsim) +KiteModels.generate_getters!(s) +KiteModels.simulate(s, logger) +save_log(logger, "tmp") +lg = load_log("tmp") +function play(s, lg, front_view = false, side_view = true) + dt = 1/s.set.sample_freq + conn = KiteModels.getconnections(s) + sl = lg.syslog + total_segmentsvector = Vector{Int64}[] + for conn_pair in conn + push!(total_segmentsvector, Int64[conn_pair[1], conn_pair[2]]) + end + # Add tether segments + for i in 0:(s.set.segments-2) + push!(total_segmentsvector, [6+i, 6+i+1]) + end + # Add final connection from last tether point to bridle point + push!(total_segmentsvector, [6+s.set.segments-1, 1]) + + for step in 1:length(0:dt:s.set.sim_time)-1 #-s.set.dt + # Get positions at this time step + x = sl.X[step] + y = sl.Y[step] + z = sl.Z[step] + + # Create points array for all points in the system + pointsvector = Vector{Float64}[] + + # Use comparison operator == instead of assignment = + if front_view == true + for i in 1:KiteModels.points(s) + # For front view: Y on horizontal axis, Z on vertical axis + push!(pointsvector, Float64[y[i], z[i], x[i]]) + end + + # Calculate appropriate limits for front view + y_min, y_max = -10, 10 + z_min, z_max = 0, 50 + t = (dt) * (step-1) + + # Plot the kite system at this time step (front view) + plot2d(pointsvector, total_segmentsvector, t; + zoom = false, + xlim = (y_min, y_max), + ylim = (z_min, z_max) + ) + elseif side_view == true + for i in 1:KiteModels.points(s) + # For side view: X on horizontal axis, Z on vertical axis + push!(pointsvector, Float64[x[i], y[i], z[i]]) + end + + # Calculate appropriate limits for side view + x_min, x_max = 0, 20 + z_min, z_max = 0, 20 + t = (dt) * (step-1) + + # Plot the kite system at this time step (side view) + plot2d(pointsvector, total_segmentsvector, t; + zoom = false, + xlim = (x_min, x_max), + ylim = (z_min, z_max) + ) + end + + # Add a small delay to control animation speed + sleep(0.05) + end + nothing +end +function plot_front_view3(lg) + display(plotxy(lg.y, lg.z; + xlabel="pos_y [m]", + ylabel="height [m]", + fig="front_view")) + nothing +end +play(s, lg) + + +# @with_kw mutable struct Settings_not_there_yet @deftype Float64 +# g_earth::Vector{Float64} = [0.0, 0.0, -9.81] # gravitational acceleration [m/s²] +# save::Bool = false # save animation frames +# # spring constants +# springconstant_tether::Float64 = 614600.0 # TETHER unit spring constant [N] +# springconstant_bridle::Float64 = 10000.0 # BRIDLE unit spring constant [N] +# springconstant_kite::Float64 = 60000.0 # KITE unit spring constant [N] +# # damping +# damping_tether::Float64 = 473 # TETHER unit damping coefficient [Ns] +# rel_damping_kite::Float64 = 6.0 # KITE relative unit damping coefficient [-] +# rel_damping_bridle:: Float64 = 6.0 # BRIDLE relative unit damping coefficient [-] +# end +# @with_kw mutable struct KPS5 +# "Reference to the settings2 struct" +# #set::Settings2 = Settings2() +# sys::Union{ModelingToolkit.ODESystem, Nothing} = nothing +# t_0::Float64 = 0.0 +# iter::Int64 = 0 +# prob::Union{OrdinaryDiffEqCore.ODEProblem, Nothing} = nothing +# integrator::Union{OrdinaryDiffEqCore.ODEIntegrator, Nothing} = nothing +# get_state::Function = () -> nothing +# end + +# # ----------------------------- +# # Simulation Function +# # ----------------------------- +# function simulate(s, logger) +# dt = 1/s.set.sample_freq +# tol = s.set.tol +# tspan = (0.0, dt) +# time_range = 0:dt:s.set.sim_time-dt +# steps = length(time_range) +# iter = 0 +# for i in 1:steps +# next_step!(s; dt=s.set.dt) +# u = s.get_state(s.integrator) +# x = u[1][1, :] +# y = u[1][2, :] +# z = u[1][3, :] +# iter += s.iter +# sys_state = SysState{points(s)}() +# sys_state.X .= x +# sys_state.Y .= y +# sys_state.Z .= z +# println("iter: $iter", " steps: $steps") +# log!(logger, sys_state) +# println(x[end], ", ", sys_state.X[end]) +# end +# println("iter: $iter", " steps: $steps") +# return nothing +# end + + +function plot_front_view3(lg) + display(plotxy(lg.y, lg.z; + xlabel="pos_y [m]", + ylabel="height [m]", + fig="front_view")) + nothing +end \ No newline at end of file diff --git a/examples/kps5_example_uwe.jl b/examples/kps5_example_uwe.jl new file mode 100644 index 000000000..e0e76c6c9 --- /dev/null +++ b/examples/kps5_example_uwe.jl @@ -0,0 +1,96 @@ +using Timers +tic() +using KiteModels +using Pkg +if ! ("ControlPlots" ∈ keys(Pkg.project().dependencies)) + using TestEnv; TestEnv.activate() +end +using ControlPlots +toc() + +include("plotting_kps5.jl") + +set = deepcopy(load_settings("system_kps5.yaml")) +kcu::KCU = KCU(set) +s = KPS5(kcu) + +dt = 1/s.set.sample_freq +time_range = 0:dt:set.sim_time-dt +steps = length(time_range) +logger = Logger(KiteModels.points(s), steps) +KiteModels.get_kite_points(s) +KiteModels.calc_initial_state(s) +init_sim!(s) +KiteModels.generate_getters!(s) +KiteModels.simulate(s, logger) +save_log(logger, "tmp") + +function play(s, lg, front_view = false, side_view = true) + dt = 1/s.set.sample_freq + conn = KiteModels.getconnections(s) + sl = lg.syslog + total_segmentsvector = Vector{Int64}[] + for conn_pair in conn + push!(total_segmentsvector, Int64[conn_pair[1], conn_pair[2]]) + end + # Add tether segments + for i in 0:(s.set.segments-2) + push!(total_segmentsvector, [6+i, 6+i+1]) + end + # Add final connection from last tether point to bridle point + push!(total_segmentsvector, [6+s.set.segments-1, 1]) + + for step in 1:length(0:dt:s.set.sim_time)-1 #-s.set.dt + # Get positions at this time step + x = sl.X[step] + y = sl.Y[step] + z = sl.Z[step] + + # Create points array for all points in the system + pointsvector = Vector{Float64}[] + + # Use comparison operator == instead of assignment = + if front_view == true + for i in 1:KiteModels.points(s) + # For front view: Y on horizontal axis, Z on vertical axis + push!(pointsvector, Float64[y[i], z[i], x[i]]) + end + + # Calculate appropriate limits for front view + y_min, y_max = -10, 10 + z_min, z_max = 0, 50 + t = (dt) * (step-1) + + # Plot the kite system at this time step (front view) + plot2d(pointsvector, total_segmentsvector, t; + zoom = false, + xlim = (y_min, y_max), + ylim = (z_min, z_max) + ) + elseif side_view == true + for i in 1:KiteModels.points(s) + # For side view: X on horizontal axis, Z on vertical axis + push!(pointsvector, Float64[x[i], y[i], z[i]]) + end + + # Calculate appropriate limits for side view + x_min, x_max = 0, 20 + z_min, z_max = 0, 20 + t = (dt) * (step-1) + + # Plot the kite system at this time step (side view) + plot2d(pointsvector, total_segmentsvector, t; + zoom = false, + xlim = (x_min, x_max), + ylim = (z_min, z_max) + ) + end + + # Add a small delay to control animation speed + sleep(0.05) + end + nothing +end + +lg = load_log("tmp") +play(s, lg) \ No newline at end of file diff --git a/examples/plotting_kps5.jl b/examples/plotting_kps5.jl new file mode 100644 index 000000000..facbc3642 --- /dev/null +++ b/examples/plotting_kps5.jl @@ -0,0 +1,7 @@ +function plot_front_view3(lg) + display(plotxy(lg.y, lg.z; + xlabel="pos_y [m]", + ylabel="height [m]", + fig="front_view")) + nothing +end diff --git a/src/KPS5.jl b/src/KPS5.jl new file mode 100644 index 000000000..271e67ccd --- /dev/null +++ b/src/KPS5.jl @@ -0,0 +1,546 @@ +#= MIT License + +Copyright (c) 2020, 2021, 2022, 2024 Uwe Fechner + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. =# + +#= Model of a kite-power system in implicit form: residual = f(y, yd) + +This model implements a 3D mass-spring system with reel-out. It uses six tether segments (the number can be +configured in the file data/settings.yaml). The kite is modelled using 4 point masses and 3 aerodynamic +surfaces. The spring constant and the damping decrease with the segment length. The aerodynamic kite forces +are acting on three of the four kite point masses. + +Four point kite model, included from KiteModels.jl. + +Scientific background: http://arxiv.org/abs/1406.6218 =# +# KiteModels.stiffness=1000 +""" + mutable struct KPS5{S, T, P, Q, SP} <: AbstractKiteModel + +State of the kite power system, using a 4 point kite model. Parameters: +- S: Scalar type, e.g. SimFloat + In the documentation mentioned as Any, but when used in this module it is always SimFloat and not Any. +- T: Vector type, e.g. MVector{3, SimFloat} +- P: number of points of the system, segments+1 +- Q: number of springs in the system, P-1 +- SP: struct type, describing a spring +Normally a user of this package will not have to access any of the members of this type directly, +use the input and output functions instead. + +$(TYPEDFIELDS) +""" +@with_kw mutable struct KPS5{S, T} <: AbstractKiteModel + "Reference to the settings struct" + set::Settings + "Reference to the KCU model (Kite Control Unit as implemented in the package KitePodModels" + kcu::KCU + "Reference to the atmospheric model as implemented in the package AtmosphericModels" + am::AtmosphericModel = AtmosphericModel() + "Reference to winch model as implemented in the package WinchModels" + wm::AbstractWinchModel + "Iterations, number of calls to the function residual!" + iter:: Int64 = 0 + "x vector of kite reference frame" + x::T = zeros(S, 3) + "y vector of kite reference frame" + y::T = zeros(S, 3) + "z vector of kite reference frame" + z::T = zeros(S, 3) + sys::Union{ModelingToolkit.ODESystem, Nothing} = nothing + t_0::Float64 = 0.0 + #iter::Int64 = 0 + prob::Union{OrdinaryDiffEqCore.ODEProblem, Nothing} = nothing + integrator::Union{OrdinaryDiffEqCore.ODEIntegrator, Nothing} = nothing + get_state::Function = () -> nothing +end +function KPS5(kcu::KCU) + if kcu.set.winch_model == "AsyncMachine" + wm = AsyncMachine(kcu.set) + elseif kcu.set.winch_model == "TorqueControlledMachine" + wm = TorqueControlledMachine(kcu.set) + end + # wm.last_set_speed = kcu.set.v_reel_out + s = KPS5{SimFloat, KVec3}(set=kcu.set, + kcu=kcu, wm=wm) + return s +end +# ------------------------------------------------------------------------------------------------------------------------------------------------------------ +# deriving a constants for points, total segments, and connections +# --------- +function points(s) + return 5 + s.set.segments +end +function total_segments(s) + return 9 + s.set.segments +end +function getconnections(s) + conn = [(1,2), (2,3), (3,1), (1,4), (2,4), (3,4), (1,5), (2,5), (3,5)] # connections between KITE points + conn = vcat(conn, [(6+i, 6+i+1) for i in 0:(s.set.segments-2)]...) # connection between tether points + conn = vcat(conn, [(6+s.set.segments-1, 1)]) # connection final tether point to bridle point + return conn +end +# ------------------------------------------------------------------------------------------------------------------------------------------------------------------------ +# Interpolating polars using Dierckx +# ------------------------------------ +alpha_cl = [-180.0, -160.0, -90.0, -20.0, -10.0, -5.0, 0.0, 20.0, 40.0, 90.0, 160.0, 180.0] +cl_list = [ 0.0, 0.5, 0.0, 0.08, 0.125, 0.15, 0.2, 1.0, 1.0, 0.0, -0.5, 0.0] +alpha_cd = [-180.0, -170.0, -140.0, -90.0, -20.0, 0.0, 20.0, 90.0, 140.0, 170.0, 180.0] +cd_list = [ 0.5, 0.5, 0.5, 1.0, 0.2, 0.1, 0.2, 1.0, 0.5, 0.5, 0.5] +function cl_interp(alpha) + cl_spline = Spline1D(alpha_cl, cl_list) + return cl_spline(alpha) +end +function cd_interp(alpha) + cd_spline = Spline1D(alpha_cd, cd_list) + return cd_spline(alpha) +end +@register_symbolic cl_interp(alpha) +@register_symbolic cd_interp(alpha) +# ----------------------------------------------------------------------------------------------------------------------------------------------------------------- +# Initialize the simulation +# ----------------------------------------- +function init_sim!(s::KPS5) + pos, vel = calc_initial_state(s) + dt = 1/s.set.sample_freq + simple_sys, pos, vel, e_x, e_y, e_z, v_app_point, alpha1p, height, wind_vector = model(s, pos, vel) + s.sys = simple_sys + tspan = (0.0, s.set.sim_time) + s.prob = ODEProblem(simple_sys, nothing, tspan) + #differential_vars = ones(Bool, length(y0)) + #prob = DAEProblem{true}(residual!, yd0, y0, tspan, s; differential_vars) + #solver = DFBDF(autodiff=AutoFiniteDiff(), max_order=Val{s.set.max_order}()) + #s.integrator = OrdinaryDiffEqCore.init(s.prob, Rodas5(autodiff=false); s.set.dt, abstol=s.set.tol, save_on=false) + #s.integrator = OrdinaryDiffEqCore.init(prob, solver; abstol=abstol, reltol=s.set.rel_tol, save_everystep=false, initializealg=OrdinaryDiffEqCore.NoInit()) + s.integrator = OrdinaryDiffEqCore.init(s.prob, FBDF(autodiff=false); dt, abstol=s.set.abs_tol, reltol = s.set.rel_tol, save_on=false) + generate_getters!(s) +end +# ------------------------------ +# Calculate Initial State +# ------------------------------ +function calc_initial_state(s::KPS5) + p1location = [s.set.l_tether*cos(deg2rad(s.set.elevation)) 0 s.set.l_tether*sin(deg2rad(s.set.elevation))] + kitepos0rot = get_kite_points(s) + POS0 = kitepos0rot .+ p1location' + POS0 = hcat(POS0, zeros(3, 1)) + if s.set.segments > 1 + extra_nodes = [POS0[:,6] + (POS0[:,1] - POS0[:,6]) * i / s.set.segments for i in 1:(s.set.segments-1)] + POS0 = hcat(POS0, extra_nodes...) + end + VEL0 = zeros(3, points(s)) + return POS0, VEL0 +end +# ----------------------------------------------------- +# initializing kite points +# ----------------------------------------------------- +function get_kite_points(s::KPS5) + # Original kite points in local reference frame + kitepos0 = # KITE points + # P1 Bridle P2 P3 P4 P5 + [0.000 s.set.cord_length/2 -s.set.cord_length/2 0 0; + 0.000 0 0 -s.set.width/2 s.set.width/2; + 0.000 s.set.height_k+s.set.h_bridle s.set.height_k+s.set.h_bridle s.set.h_bridle s.set.h_bridle] + + beta = deg2rad(s.set.elevation) + Y_r = [sin(beta) 0 cos(beta); + 0 1 0; + -cos(beta) 0 sin(beta)] + # Apply rotation to all points + kitepos0rot = Y_r * kitepos0 + + return kitepos0rot +end +# ------------------------------------------------------------------------------------------------------------------------------------------------------------------------- +# Define a function to create reference frame update equations +# --------------------------------------------------------------- +function create_reference_frame_equations(pos, e_x, e_y, e_z) + # Calculate vectors for the reference frame + X = pos[:, 2] - pos[:, 3] # Vector from P3 to P2 + Y = pos[:, 5] - pos[:, 4] # Vector from P4 to P5 + Z = cross(X, Y) # Cross product for Z axis + # Normalize these vectors to get unit vectors + norm_X = sqrt(sum(X .^ 2)) + norm_Y = sqrt(sum(Y .^ 2)) + norm_Z = sqrt(sum(Z .^ 2)) + # Create equations to update the reference frame + ref_frame_eqs = [ + e_x[1] ~ X[1] / norm_X, + e_x[2] ~ X[2] / norm_X, + e_x[3] ~ X[3] / norm_X, + e_y[1] ~ Y[1] / norm_Y, + e_y[2] ~ Y[2] / norm_Y, + e_y[3] ~ Y[3] / norm_Y, + e_z[1] ~ Z[1] / norm_Z, + e_z[2] ~ Z[2] / norm_Z, + e_z[3] ~ Z[3] / norm_Z + ] + return ref_frame_eqs +end +# -------------------------------------------------------------------------- +# computing angle of attack +# -------------------------------------------------------------------------- +function compute_alpha1p(v_a, e_z, e_x) + # Calculate the angle Alpha1p + v_a_z = v_a ⋅ e_z # Use the symbolic dot product + v_a_x = v_a ⋅ e_x # Use the symbolic dot product + + # Use atan for the symbolic computation + alpha1p =rad2deg.(atan(v_a_z, v_a_x)) + return alpha1p +end +# ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- +# Calculate Rest Lengths +# --------------------------- +function calc_rest_lengths(s::KPS5) + conn = getconnections(s) + POS0, VEL0 = calc_initial_state(s) + lengths = [norm(POS0[:,conn[i][2]] - POS0[:,conn[i][1]]) for i in 1:9] + l10 = norm(POS0[:,1] - POS0[:,6]) + lengths = vcat(lengths, [(l10 + s.set.v_reel_out*t)/s.set.segments for _ in 1:s.set.segments]...) + return lengths, l10 +end +function get_wind_vector(s::KPS5, v_wind_height) + v_wind_magnitude = v_wind_height # add function + wind_angle = deg2rad(s.set.upwind_dir) # angle measuring form North (pos x), CW(+) + wind_vector = v_wind_magnitude*[-cos(wind_angle),-sin(wind_angle), 0.0] + return wind_vector +end +# ------------------------------------------------------------------------------------------------------------------------------------------------------------------ +#Define the Model +# ----------------------------------------------- +function model(s::KPS5, pos, vel) + POS0, VEL0 = pos, vel + rest_lengths, l_tether = calc_rest_lengths(s) + # @parameters K1=s.set.springconstant_tether K2=s.set.springconstant_bridle K3=s.set.springconstant_kite C1=s.set.damping_tether C2=s.set.rel_damping_bridle*s.set.damping_tether C3=s.set.rel_damping_kite*s.set.damping_tether + # same as Uwe here, only rel c and k , wrt tether, can/should be improved, to be for kite/bridle separately + @parameters E_modulus_tether=s.set.e_tether K3=s.set.c_spring_kite C1=s.set.damping C3=s.set.damping_kite_springs + @parameters m_kite=s.set.mass kcu_mass=s.set.kcu_mass rho_tether=s.set.rho_tether rel_compr_k=s.set.rel_compr_stiffness + @parameters rho=s.set.rho_0 g_earth=-9.81 cd_tether=s.set.cd_tether d_tether=s.set.d_tether d_bridle_line=s.set.d_line S=s.set.area l_bridle=s.set.l_bridle + @parameters kcu_cd=s.set.cd_kcu kcu_diameter=s.set.kcu_diameter v_wind_ground=s.set.v_wind + @variables pos(t)[1:3, 1:points(s)] = POS0 + @variables vel(t)[1:3, 1:points(s)] = VEL0 + @variables acc(t)[1:3, 1:points(s)] + @variables v_app_point(t)[1:3, 1:points(s)] + @variables segment(t)[1:3, 1:total_segments(s)] + @variables unit_vector(t)[1:3, 1:total_segments(s)] + @variables norm1(t)[1:total_segments(s)] + @variables height(t)[1:total_segments(s)] + @variables v_wind_height(t)[1:total_segments(s)] + @variables wind_vector(t)[1:3, 1:total_segments(s)] + @variables rel_vel(t)[1:3, 1:total_segments(s)] + @variables spring_vel(t)[1:total_segments(s)] + @variables k_spring(t)[1:total_segments(s)] + @variables spring_force(t)[1:3, 1:total_segments(s)] + @variables v_apparent(t)[1:3, 1:total_segments(s)] + @variables v_app_perp(t)[1:3, 1:total_segments(s)] + @variables norm_v_app(t)[1:total_segments(s)] + @variables half_drag_force(t)[1:3, 1:total_segments(s)] + @variables drag_force(t)[1:3, 1:total_segments(s)] + @variables total_force(t)[1:3, 1:points(s)] + # local kite reference frame + @variables e_x(t)[1:3] + @variables e_y(t)[1:3] + @variables e_z(t)[1:3] + @variables alpha1p(t)[1:1] + + eqs1 = vcat(D.(pos) .~ vel, + D.(vel) .~ acc) + eqs2 = vcat(eqs1...) + eqs2 = vcat(eqs2, acc[:,6] .~ [0.0, 0.0, 0.0]) # origin is six, make 6 not being hardcoded + # ----------------------------- + # defining the connections and their respective rest lengths, unit spring constants, damping and masses + # ----------------------------- + # connections adding segment connections, from origin to bridle + conn = getconnections(s) + # final connection last tether point to bridle point + # unit spring constants (K1 tether, K2 bridle, K3 kite) K3 from settings + K1 = E_modulus_tether*pi*(d_tether/2000)^2 + K2 = E_modulus_tether*pi*(d_bridle_line/2000)^2 + k_segments = [K2, K3, K2, K2, K3, K3, K2, K3, K3] + k_segments = vcat(k_segments, [K1 for _ in 1:s.set.segments]...) + # unit damping constants (C1 tether, C2 bridle, C3 kite) C1 and C3 from settings + totalbridlelengths_model = rest_lengths[1]+rest_lengths[3]+rest_lengths[4]+rest_lengths[7] + lines_per_segment = l_bridle/totalbridlelengths_model + C2 = lines_per_segment*C1*(d_bridle_line/d_tether)^2 + c_segments = [C2, C3, C2, C2, C3, C3, C2, C3, C3] + c_segments = vcat(c_segments, [C1 for _ in 1:s.set.segments]...) + # masses + mass_bridlelines = ((d_bridle_line/2000)^2)*pi*rho_tether*s.set.l_bridle #total mass entire bridle + mass_halfbridleline = mass_bridlelines/8 # half the connection of bridle line to kite (to assign to each kitepoint) so the other 4 halves get assigned to bridlepoint + mass_tether = ((d_tether/2000)^2)*pi*rho_tether*l_tether + mass_tetherpoints = mass_tether/(s.set.segments+1) + mass_bridlepoint = 4*mass_halfbridleline + kcu_mass + mass_tetherpoints # 4 bridle connections, kcu and tether + m_kitepoints = [(m_kite*s.set.rel_mass_p2) + mass_halfbridleline, (m_kite*s.set.rel_mass_p3) + mass_halfbridleline, (m_kite*s.set.rel_mass_p4) + mass_halfbridleline, (m_kite*s.set.rel_mass_p4) + mass_halfbridleline] # mass P4=P5 + PointMasses = [mass_bridlepoint, m_kitepoints[1], m_kitepoints[2], m_kitepoints[3], m_kitepoints[4]] + PointMasses = vcat(PointMasses, [mass_tetherpoints for _ in 1:s.set.segments]...) + # ----------------------------- + # Equations for Each Segment (Spring Forces, Drag, etc.) + # ----------------------------- + for i in 1:total_segments(s) + eqs = [ + segment[:, i] ~ pos[:, conn[i][2]] - pos[:, conn[i][1]], + norm1[i] ~ norm(segment[:, i]), + height[i] ~ max((pos[3, conn[i][2]] + pos[3, conn[i][1]])/2, 6.0), # 6m is the minimum height + v_wind_height[i] ~ calc_wind_factor(s.am, height[i])*v_wind_ground, + wind_vector[:, i] ~ get_wind_vector(s, v_wind_height[i]), + unit_vector[:, i] ~ -segment[:, i] / norm1[i], + rel_vel[:, i] ~ vel[:, conn[i][2]] - vel[:, conn[i][1]], + spring_vel[i] ~ -unit_vector[:, i] ⋅ rel_vel[:, i], + k_spring[i] ~ (k_segments[i]/rest_lengths[i]) * (rel_compr_k + (1-rel_compr_k)*(norm1[i] > rest_lengths[i])), + spring_force[:, i] ~ (k_spring[i]*(norm1[i] - rest_lengths[i]) + c_segments[i] * spring_vel[i]) * unit_vector[:, i], + v_apparent[:, i] ~ wind_vector[:, i] - (vel[:, conn[i][1]] + vel[:, conn[i][2]]) / 2, # change wind for winddirection! integrate + v_app_perp[:, i] ~ v_apparent[:, i] - (v_apparent[:, i] ⋅ unit_vector[:, i]) .* unit_vector[:, i], + norm_v_app[i] ~ norm(v_app_perp[:, i]) + ] + if i > 9 # tether segments + push!(eqs, half_drag_force[:, i] ~ 0.25 * rho * cd_tether * norm_v_app[i] * (rest_lengths[i]*d_tether/1000) * v_app_perp[:, i]) + elseif i in [1, 3, 4, 7] # bridle lines, try to find Cd_bridlelines later + push!(eqs, half_drag_force[:, i] ~ 0.25 * rho * cd_tether * norm_v_app[i] * (rest_lengths[i]*(s.set.d_line/1000)) * v_app_perp[:, i]) + else # kite + push!(eqs, half_drag_force[:, i] ~ zeros(3)) + end + eqs2 = vcat(eqs2, reduce(vcat, eqs)) + end + # ----------------------------- + # Reference Frame and Aerodynamic Coefficients + # ----------------------------- + ref_frame_eqs = create_reference_frame_equations(pos, e_x, e_y, e_z) + eqs2 = vcat(eqs2, ref_frame_eqs) + # only 1 AOA + v_a_kite = wind_vector[:, 2] - (vel[:, 2] + vel[:, 3])/2 # computing AOA only for center chord # Appas.set.t wind velocity + alpha1p = compute_alpha1p(v_a_kite, e_z, e_x) # Calculate Alpha1p at this time step + eqs2 = vcat(eqs2, alpha1p[1] ~ alpha1p) # Add the equation for Alpha1p for each of 4 kite points (first bering bridle so i-1) + # getting Cl and Cd + Cl = cl_interp(alpha1p) + Cd = cd_interp(alpha1p) + + # ----------------------------- + # Force Balance at Each Point + # ----------------------------- + for i in 1:points(s) + eqs = [] + force = sum([spring_force[:, j] for j in 1:total_segments(s) if conn[j][2] == i]; init=zeros(3)) - + sum([spring_force[:, j] for j in 1:total_segments(s) if conn[j][1] == i]; init=zeros(3)) + + sum([half_drag_force[:, j] for j in 1:total_segments(s) if conn[j][1] == i]; init=zeros(3)) + + sum([half_drag_force[:, j] for j in 1:total_segments(s) if conn[j][2] == i]; init=zeros(3)) + v_app_point[:, i] ~ wind_vector[:, i] - vel[:, i] + if i == 1 # KCU drag at bridle point + area_kcu = pi * ((kcu_diameter / 2) ^ 2) + Dx_kcu = 0.5*rho*kcu_cd *area_kcu*(v_app_point[1, i]*v_app_point[1, i]) + Dy_kcu = 0.5*rho*kcu_cd *area_kcu*(v_app_point[2, i]*v_app_point[2, i]) + Dz_kcu = 0.5*rho*kcu_cd *area_kcu*(v_app_point[3, i]*v_app_point[3, i]) + D = [Dx_kcu, Dy_kcu, Dz_kcu] + push!(eqs, total_force[:, i] ~ force + D) + elseif i in 2:5 # the kite points that get Aero Forces + + v_app_mag_squared = v_app_point[1, i]^2 + v_app_point[2, i]^2 + v_app_point[3, i]^2 + # # Lift calculation + L_perpoint = (1/4) * 0.5 * rho * Cl * S * (v_app_mag_squared) + # Cross product and normalization + cross_vapp_X_e_y = cross(v_app_point[:, i], e_y) + normcross_vapp_X_e_y = norm(cross_vapp_X_e_y) + L_direction = cross_vapp_X_e_y / normcross_vapp_X_e_y + # Final lift force vector + L = L_perpoint * L_direction + + # Drag calculation + D_perpoint = (1/4) * 0.5 * rho * Cd * S * v_app_mag_squared + # # Create drag direction components + D_direction = [v_app_point[1, i] / norm(v_app_point[:, i]), v_app_point[2, i] / norm(v_app_point[:, i]), v_app_point[3, i] / norm(v_app_point[:, i])] + # # Final drag force vector components + D = D_perpoint * D_direction + + # Total aerodynamic force + Fa = [L[1]+ D[1], L[2]+ D[2], L[3]+ D[3]] + push!(eqs, total_force[1, i] ~ force[1] + Fa[1]) + push!(eqs, total_force[2, i] ~ force[2] + Fa[2]) + push!(eqs, total_force[3, i] ~ force[3] + Fa[3]) + elseif i != 6 + push!(eqs, total_force[:, i] ~ force) + end + push!(eqs, acc[:, i] ~ [0.0, 0.0, g_earth] + total_force[:, i] / PointMasses[i]) + eqs2 = vcat(eqs2, reduce(vcat, eqs)) + eqs2 = vcat(eqs2, v_app_point[:, i] ~ wind_vector[:, i] - vel[:, i]) + end + @named sys = ODESystem(reduce(vcat, Symbolics.scalarize.(eqs2)), t) + simple_sys = structural_simplify(sys) + simple_sys, pos, vel, e_x, e_y, e_z, v_app_point, alpha1p, height, wind_vector, norm1 +end +# ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- +# next step function +# ----------------------------- +function next_step!(s::KPS5; dt=(1/s.set.sample_freq)) + s.t_0 = s.integrator.t + steptime = @elapsed OrdinaryDiffEqCore.step!(s.integrator, dt, true) + s.iter += 1 + # Get the current state including wind vector + state = s.get_state(s.integrator) + if !successful_retcode(s.integrator.sol) + println("Return code for solution: ", s.integrator.sol.retcode) + end + if length(state) >= 2 # Make sure wind_vector is included in state + height = state[2] # The height from the state + wind_vec = state[3] + norm1 = state[4] # The norm of the segments from the state + # The wind vector from the state + #println("Height at each segment : $(height[:, :])") + #println("Time: $(s.integrator.t), Wind Vector at each segment: $(wind_vec[:, :])") + println("Norm of the segments at each segment: $(norm1[:, :])") + end + s.integrator.t, steptime +end +function simulate(s, logger) + dt = 1/s.set.sample_freq + tol = s.set.abs_tol + tspan = (0.0, dt) + time_range = 0:dt:s.set.sim_time-dt + steps = length(time_range) + iter = 0 + for i in 1:steps + next_step!(s; dt=dt) + u = s.get_state(s.integrator) + x = u[1][1, :] + y = u[1][2, :] + z = u[1][3, :] + iter += s.iter + sys_state = SysState{points(s)}() + sys_state.X .= x + sys_state.Y .= y + sys_state.Z .= z + println("iter: $iter", " steps: $steps") + log!(logger, sys_state) + println(x[end], ", ", sys_state.X[end]) + end + println("iter: $iter", " steps: $steps") + return nothing +end + +# ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- +# Generate Getters +# ----------------------------- +function generate_getters!(s::KPS5) + sys = s.sys + c = collect + get_state = ModelingToolkit.getu(sys, + [c(sys.pos), c(sys.height) , c(sys.wind_vector), c(sys.norm1)] + ) + s.get_state = (integ) -> get_state(integ) + return nothing +end + + +# =================== getter functions ==================================================== + +""" + calc_height(s::KPS5) + +Determine the height of the topmost kite particle above ground. +""" +function calc_height(s::KPS5) + pos_kite(s)[3] +end + +""" + pos_kite(s::KPS5) + +Return the position of the kite (top particle). +""" +function pos_kite(s::KPS5) + s.pos[end-2] +end + +""" + kite_ref_frame(s::KPS5; one_point=false) + +Returns a tuple of the x, y, and z vectors of the kite reference frame. +""" +function kite_ref_frame(s::KPS5; one_point=false) + if one_point + c = s.z + y = normalize(s.v_apparent × c) + x = normalize(y × c) + return x, y, c + else + return s.x, s.y, s.z + end +end + +""" + winch_force(s::KPS5) + +Return the absolute value of the force at the winch as calculated during the last timestep. +""" +function winch_force(s::KPS5) norm(s.last_force) end + +""" + cl_cd(s::KPS5) + +Calculate the lift and drag coefficients of the kite, based on the current angles of attack. +""" +function cl_cd(s::KPS5) + rel_side_area = s.set.rel_side_area/100.0 # defined in percent + K = 1 - rel_side_area # correction factor for the drag + if s.set.version == 3 + drag_corr = 1.0 + else + drag_corr = DRAG_CORR + end + CL2, CD2 = s.calc_cl(s.alpha_2), drag_corr * s.calc_cd(s.alpha_2) + CL3, CD3 = s.calc_cl(s.alpha_3), drag_corr * s.calc_cd(s.alpha_3) + CL4, CD4 = s.calc_cl(s.alpha_4), drag_corr * s.calc_cd(s.alpha_4) + if s.set.version == 3 + return CL2, CD2 + else + return CL2, K*(CD2+CD3+CD4) + end +end + +# ==================== end of getter functions ================================================ + +function spring_forces(s::KPS5; prn=true) + forces = zeros(SimFloat, s.set.segments+KITE_SPRINGS) + for i in 1:s.set.segments + forces[i] = s.springs[i].c_spring * (norm(s.pos[i+1] - s.pos[i]) - s.segment_length) * s.stiffness_factor + if forces[i] > s.set.max_force && prn + println("Tether raptures for segment $i !") + end + end + for i in 1:KITE_SPRINGS + p1 = s.springs[i+s.set.segments].p1 # First point nr. + p2 = s.springs[i+s.set.segments].p2 # Second point nr. + pos1, pos2 = s.pos[p1], s.pos[p2] + spring = s.springs[i+s.set.segments] + l_0 = spring.length # Unstressed length + k = spring.c_spring * s.stiffness_factor # Spring constant + segment = pos1 - pos2 + norm1 = norm(segment) + k1 = 0.25 * k # compression stiffness kite segments + if (norm1 - l_0) > 0.0 + spring_force = k * (norm1 - l_0) + else + spring_force = k1 * (norm1 - l_0) + end + forces[i+s.set.segments] = spring_force + if norm(s.spring_force) > 4000.0 + println("Bridle brakes for spring $i !") + end + end + forces +end + diff --git a/src/KPS5Friso.jl b/src/KPS5Friso.jl new file mode 100644 index 000000000..001b9a566 --- /dev/null +++ b/src/KPS5Friso.jl @@ -0,0 +1,452 @@ +# average AOA is used, so now only 1 value, not 4, making editable +# clean up code +# New Plot +# Using logger +# segments: now all segments +using Timers +tic() +using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra, Timers, Parameters, ControlPlots +using ModelingToolkit: Symbolics, @register_symbolic +using OrdinaryDiffEqCore +using Dierckx +using ModelingToolkit: t_nounits as t, D_nounits as D +using KiteUtils +include("plots.jl") +toc() +@with_kw mutable struct Settings2 @deftype Float64 + g_earth::Vector{Float64} = [0.0, 0.0, -9.81] # gravitational acceleration [m/s²] + v_wind_tether::Vector{Float64} = [13.5, 0.0, 0.0] # wind velocity [m/s] + rho::Float64 = 1.225 # air density [kg/m³] + sim_time::Float64 = 5 # simulation duration [s] + dt = 0.05 # time step [s] + tol = 1e-6 # tolerance for the solver + save::Bool = false # save animation frames + # kite + mass::Float64 = 10.58 # mass of kite [kg] + Area::Float64 = 20.36 # surface area [m²] + width::Float64 = 8.16 # width of kite [m] + height::Float64 = 3.15 # height of kite [m] + chord_length::Float64 = 2.0 # chord length [m] + # KCU + bridle + h_bridle = 4.9 # height of bridle [m] + d_line::Float64 = 2.5 # bridle line diameter [mm] + l_bridle::Float64 = 33.4 # sum of the lengths of the bridle lines [m] + kcu_cd::Float64 = 0.47 # KCU drag coefficient + kcu_diameter::Float64 = 0.4 # KCU diameter [m] + kcu_mass::Float64 = 11 # mass of KCU (at bridle point)[kg] + # spring constants + springconstant_tether::Float64 = 614600.0 # TETHER unit spring constant [N] + springconstant_bridle::Float64 = 10000.0 # BRIDLE unit spring constant [N] + springconstant_kite::Float64 = 60000.0 # KITE unit spring constant [N] + # damping + damping_tether::Float64 = 473 # TETHER unit damping coefficient [Ns] + rel_damping_kite::Float64 = 6.0 # KITE relative unit damping coefficient [-] + rel_damping_bridle:: Float64 = 6.0 # BRIDLE relative unit damping coefficient [-] + # tether + l_tether::Float64 = 50.0 # tether length [m] + v_reel_out ::Float64 = 0.5 # reel-out speed [m/s] + rho_tether::Float64 = 724.0 # density of tether [kg/m³] + cd_tether::Float64 = 0.958 # drag coefficient of tether + d_tether::Float64 = 4 # tether diameter [mm] + elevation::Float64 = 1.22173048 # Elevation angle, angle XZ plane, between origin and bridle point [rad] + segments::Int64 = 14 # number of tether segments [-] +end +@with_kw mutable struct KPS5 + "Reference to the settings2 struct" + set::Settings2 = Settings2() + sys::Union{ModelingToolkit.ODESystem, Nothing} = nothing + t_0::Float64 = 0.0 + iter::Int64 = 0 + prob::Union{OrdinaryDiffEqCore.ODEProblem, Nothing} = nothing + integrator::Union{OrdinaryDiffEqCore.ODEIntegrator, Nothing} = nothing + get_state::Function = () -> nothing +end +# ------------------------------------------------------------------------------------------------------------------------------------------------------------ +# deriving a constants for points, total segments, and connections +# --------- +function points(s) + return 5 + s.set.segments +end +function total_segments(s) + return 9 + s.set.segments +end +function getconnections(s) + conn = [(1,2), (2,3), (3,1), (1,4), (2,4), (3,4), (1,5), (2,5), (3,5)] # connections between KITE points + conn = vcat(conn, [(6+i, 6+i+1) for i in 0:(s.set.segments-2)]...) # connection between tether points + conn = vcat(conn, [(6+s.set.segments-1, 1)]) # connection final tether point to bridle point + return conn +end +# ------------------------------------------------------------------------------------------------------------------------------------------------------------------------ +# Interpolating polars using Dierckx +# ------------------------------------ +alpha_cl = [-180.0, -160.0, -90.0, -20.0, -10.0, -5.0, 0.0, 20.0, 40.0, 90.0, 160.0, 180.0] +cl_list = [ 0.0, 0.5, 0.0, 0.08, 0.125, 0.15, 0.2, 1.0, 1.0, 0.0, -0.5, 0.0] +alpha_cd = [-180.0, -170.0, -140.0, -90.0, -20.0, 0.0, 20.0, 90.0, 140.0, 170.0, 180.0] +cd_list = [ 0.5, 0.5, 0.5, 1.0, 0.2, 0.1, 0.2, 1.0, 0.5, 0.5, 0.5] +function cl_interp(alpha) + cl_spline = Spline1D(alpha_cl, cl_list) + return cl_spline(alpha) +end +function cd_interp(alpha) + cd_spline = Spline1D(alpha_cd, cd_list) + return cd_spline(alpha) +end +@register_symbolic cl_interp(alpha) +@register_symbolic cd_interp(alpha) +# ----------------------------------------------------------------------------------------------------------------------------------------------------------------- +# Initialize the simulation +# ----------------------------------------- +function init_sim!(s::KPS5) + pos, vel = calc_initial_state(s) + simple_sys, pos, vel, e_x, e_y, e_z, v_app_point, alpha1p = model(s, pos, vel) + s.sys = simple_sys + tspan = (0.0, s.set.sim_time) + s.prob = ODEProblem(simple_sys, nothing, tspan) + s.integrator = OrdinaryDiffEqCore.init(s.prob, Rodas5(autodiff=false); s.set.dt, abstol=s.set.tol, save_on=false) +end +# ------------------------------ +# Calculate Initial State +# ------------------------------ +function calc_initial_state(s) + p1location = [s.set.l_tether*cos(s.set.elevation) 0 s.set.l_tether*sin(s.set.elevation)] + kitepos0rot = get_kite_points(s) + POS0 = kitepos0rot .+ p1location' + POS0 = hcat(POS0, zeros(3, 1)) + if s.set.segments > 1 + extra_nodes = [POS0[:,6] + (POS0[:,1] - POS0[:,6]) * i / s.set.segments for i in 1:(s.set.segments-1)] + POS0 = hcat(POS0, extra_nodes...) + end + VEL0 = zeros(3, points(s)) + return POS0, VEL0 +end +# ----------------------------------------------------- +# initializing kite points +# ----------------------------------------------------- +function get_kite_points(s) + # Original kite points in local reference frame + kitepos0 = # KITE points + # P1 Bridle P2 P3 P4 P5 + [0.000 s.set.chord_length/2 -s.set.chord_length/2 0 0; + 0.000 0 0 -s.set.width/2 s.set.width/2; + 0.000 s.set.height+s.set.h_bridle s.set.height+s.set.h_bridle s.set.h_bridle s.set.h_bridle] + + beta = s.set.elevation + Y_r = [sin(beta) 0 cos(beta); + 0 1 0; + -cos(beta) 0 sin(beta)] + # Apply rotation to all points + kitepos0rot = Y_r * kitepos0 + + return kitepos0rot +end +# ------------------------------------------------------------------------------------------------------------------------------------------------------------------------- +# Define a function to create reference frame update equations +# --------------------------------------------------------------- +function create_reference_frame_equations(pos, e_x, e_y, e_z) + # Calculate vectors for the reference frame + X = pos[:, 2] - pos[:, 3] # Vector from P3 to P2 + Y = pos[:, 5] - pos[:, 4] # Vector from P4 to P5 + Z = cross(X, Y) # Cross product for Z axis + # Normalize these vectors to get unit vectors + norm_X = sqrt(sum(X .^ 2)) + norm_Y = sqrt(sum(Y .^ 2)) + norm_Z = sqrt(sum(Z .^ 2)) + # Create equations to update the reference frame + ref_frame_eqs = [ + e_x[1] ~ X[1] / norm_X, + e_x[2] ~ X[2] / norm_X, + e_x[3] ~ X[3] / norm_X, + e_y[1] ~ Y[1] / norm_Y, + e_y[2] ~ Y[2] / norm_Y, + e_y[3] ~ Y[3] / norm_Y, + e_z[1] ~ Z[1] / norm_Z, + e_z[2] ~ Z[2] / norm_Z, + e_z[3] ~ Z[3] / norm_Z + ] + return ref_frame_eqs +end +# -------------------------------------------------------------------------- +# computing angle of attack +# -------------------------------------------------------------------------- +function compute_alpha1p(v_a, e_z, e_x) + # Calculate the angle Alpha1p + v_a_z = v_a ⋅ e_z # Use the symbolic dot product + v_a_x = v_a ⋅ e_x # Use the symbolic dot product + + # Use atan for the symbolic computation + alpha1p =rad2deg.(atan(v_a_z, v_a_x)) + return alpha1p +end +# ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- +# Calculate Rest Lengths +# --------------------------- +function calc_rest_lengths(s) + conn = getconnections(s) + POS0, VEL0 = calc_initial_state(s) + lengths = [norm(POS0[:,conn[i][2]] - POS0[:,conn[i][1]]) for i in 1:9] + l10 = norm(POS0[:,1] - POS0[:,6]) + lengths = vcat(lengths, [(l10 + s.set.v_reel_out*t)/s.set.segments for _ in 1:s.set.segments]...) + return lengths, l10 +end +# ------------------------------------------------------------------------------------------------------------------------------------------------------------------ +#Define the Model +# ----------------------------------------------- +function model(s, pos, vel) + POS0, VEL0 = pos, vel + rest_lengths, l_tether = calc_rest_lengths(s) + @parameters K1=s.set.springconstant_tether K2=s.set.springconstant_bridle K3=s.set.springconstant_kite C1=s.set.damping_tether C2=s.set.rel_damping_bridle*s.set.damping_tether C3=s.set.rel_damping_kite*s.set.damping_tether + @parameters m_kite=s.set.mass kcu_mass=s.set.kcu_mass rho_tether=s.set.rho_tether + @parameters rho=s.set.rho cd_tether=s.set.cd_tether d_tether=s.set.d_tether S=s.set.Area + @parameters kcu_cd=s.set.kcu_cd kcu_diameter=s.set.kcu_diameter + @variables pos(t)[1:3, 1:points(s)] = POS0 + @variables vel(t)[1:3, 1:points(s)] = VEL0 + @variables acc(t)[1:3, 1:points(s)] + @variables v_app_point(t)[1:3, 1:points(s)] + @variables segment(t)[1:3, 1:total_segments(s)] + @variables unit_vector(t)[1:3, 1:total_segments(s)] + @variables norm1(t)[1:total_segments(s)] + @variables rel_vel(t)[1:3, 1:total_segments(s)] + @variables spring_vel(t)[1:total_segments(s)] + @variables k_spring(t)[1:total_segments(s)] + @variables spring_force(t)[1:3, 1:total_segments(s)] + @variables v_apparent(t)[1:3, 1:total_segments(s)] + @variables v_app_perp(t)[1:3, 1:total_segments(s)] + @variables norm_v_app(t)[1:total_segments(s)] + @variables half_drag_force(t)[1:3, 1:total_segments(s)] + @variables drag_force(t)[1:3, 1:total_segments(s)] + @variables total_force(t)[1:3, 1:points(s)] + # local kite reference frame + @variables e_x(t)[1:3] + @variables e_y(t)[1:3] + @variables e_z(t)[1:3] + @variables alpha1p(t)[1:1] + + eqs1 = vcat(D.(pos) .~ vel, + D.(vel) .~ acc) + eqs2 = vcat(eqs1...) + eqs2 = vcat(eqs2, acc[:,6] .~ [0.0, 0.0, 0.0]) # origin is six, make 6 not being hardcoded + # ----------------------------- + # defining the connections and their respective rest lengths, unit spring constants, damping and masses + # ----------------------------- + # connections adding segment connections, from origin to bridle + conn = getconnections(s) + # final connection last tether point to bridle point + # unit spring constants (K1 tether, K2 bridle, K3 kite) + k_segments = [K2, K3, K2, K2, K3, K3, K2, K3, K3] + k_segments = vcat(k_segments, [K1 for _ in 1:s.set.segments]...) + # unit damping constants (C1 tether, C2 bridle, C3 kite) + c_segments = [C2, C3, C2, C2, C3, C3, C2, C3, C3] + c_segments = vcat(c_segments, [C1 for _ in 1:s.set.segments]...) + # masses + mass_bridlelines = ((s.set.d_line/2000)^2)*pi*rho_tether*s.set.l_bridle #total mass entire bridle + mass_halfbridleline = mass_bridlelines/8 # half the connection of bridle line to kite (to assign to each kitepoint) so the other 4 halves get assigned to bridlepoint + mass_tether = ((d_tether/2000)^2)*pi*rho_tether*l_tether + mass_tetherpoints = mass_tether/(s.set.segments+1) + mass_bridlepoint = 4*mass_halfbridleline + kcu_mass + mass_tetherpoints # 4 bridle connections, kcu and tether + m_kitepoints = (m_kite/4) + mass_halfbridleline + PointMasses = [mass_bridlepoint, m_kitepoints, m_kitepoints, m_kitepoints, m_kitepoints] + PointMasses = vcat(PointMasses, [mass_tetherpoints for _ in 1:s.set.segments]...) + # ----------------------------- + # Equations for Each Segment (Spring Forces, Drag, etc.) + # ----------------------------- + for i in 1:total_segments(s) + eqs = [ + segment[:, i] ~ pos[:, conn[i][2]] - pos[:, conn[i][1]], + norm1[i] ~ norm(segment[:, i]), + unit_vector[:, i] ~ -segment[:, i] / norm1[i], + rel_vel[:, i] ~ vel[:, conn[i][2]] - vel[:, conn[i][1]], + spring_vel[i] ~ -unit_vector[:, i] ⋅ rel_vel[:, i], + k_spring[i] ~ (k_segments[i]/rest_lengths[i]) * (0.1 + 0.9*(norm1[i] > rest_lengths[i])), + spring_force[:, i] ~ (k_spring[i]*(norm1[i] - rest_lengths[i]) + c_segments[i] * spring_vel[i]) * unit_vector[:, i], + v_apparent[:, i] ~ s.set.v_wind_tether .- (vel[:, conn[i][1]] + vel[:, conn[i][2]]) / 2, + v_app_perp[:, i] ~ v_apparent[:, i] - (v_apparent[:, i] ⋅ unit_vector[:, i]) .* unit_vector[:, i], + norm_v_app[i] ~ norm(v_app_perp[:, i]) + ] + if i > 9 # tether segments + push!(eqs, half_drag_force[:, i] ~ 0.25 * rho * cd_tether * norm_v_app[i] * (rest_lengths[i]*d_tether/1000) * v_app_perp[:, i]) + elseif i in [1, 3, 4, 7] # bridle lines, try to find Cd_bridlelines later + push!(eqs, half_drag_force[:, i] ~ 0.25 * rho * cd_tether * norm_v_app[i] * (rest_lengths[i]*(s.set.d_line/1000)) * v_app_perp[:, i]) + else # kite + push!(eqs, half_drag_force[:, i] ~ zeros(3)) + end + eqs2 = vcat(eqs2, reduce(vcat, eqs)) + end + # ----------------------------- + # Reference Frame and Aerodynamic Coefficients + # ----------------------------- + ref_frame_eqs = create_reference_frame_equations(pos, e_x, e_y, e_z) + eqs2 = vcat(eqs2, ref_frame_eqs) + # only 1 AOA + v_a_kite = s.set.v_wind_tether - (vel[:, 2] + vel[:, 3])/2 # computing AOA only for center chord # Appas.set.t wind velocity + alpha1p = compute_alpha1p(v_a_kite, e_z, e_x) # Calculate Alpha1p at this time step + eqs2 = vcat(eqs2, alpha1p[1] ~ alpha1p) # Add the equation for Alpha1p for each of 4 kite points (first bering bridle so i-1) + # getting Cl and Cd + Cl = cl_interp(alpha1p) + Cd = cd_interp(alpha1p) + + # ----------------------------- + # Force Balance at Each Point + # ----------------------------- + for i in 1:points(s) + eqs = [] + force = sum([spring_force[:, j] for j in 1:total_segments(s) if conn[j][2] == i]; init=zeros(3)) - + sum([spring_force[:, j] for j in 1:total_segments(s) if conn[j][1] == i]; init=zeros(3)) + + sum([half_drag_force[:, j] for j in 1:total_segments(s) if conn[j][1] == i]; init=zeros(3)) + + sum([half_drag_force[:, j] for j in 1:total_segments(s) if conn[j][2] == i]; init=zeros(3)) + v_app_point[:, i] ~ s.set.v_wind_tether - vel[:, i] + if i == 1 # KCU drag at bridle point + area_kcu = pi * ((kcu_diameter / 2) ^ 2) + Dx_kcu = 0.5*rho*kcu_cd *area_kcu*(v_app_point[1, i]*v_app_point[1, i]) + Dy_kcu = 0.5*rho*kcu_cd *area_kcu*(v_app_point[2, i]*v_app_point[2, i]) + Dz_kcu = 0.5*rho*kcu_cd *area_kcu*(v_app_point[3, i]*v_app_point[3, i]) + D = [Dx_kcu, Dy_kcu, Dz_kcu] + push!(eqs, total_force[:, i] ~ force + D) + elseif i in 2:5 # the kite points that get Aero Forces + + v_app_mag_squared = v_app_point[1, i]^2 + v_app_point[2, i]^2 + v_app_point[3, i]^2 + # Lift calculation + L_perpoint = (1/4) * 0.5 * rho * Cl * S * (v_app_mag_squared) + # Cross product and normalization + cross_vapp_X_e_y = cross(v_app_point[:, i], e_y) + normcross_vapp_X_e_y = norm(cross_vapp_X_e_y) + L_direction = cross_vapp_X_e_y / normcross_vapp_X_e_y + # Final lift force vector + L = L_perpoint * L_direction + + # Drag calculation + D_perpoint = (1/4) * 0.5 * rho * Cd * S * v_app_mag_squared + # Create drag direction components + D_direction = [v_app_point[1, i] / norm(v_app_point[:, i]), v_app_point[2, i] / norm(v_app_point[:, i]), v_app_point[3, i] / norm(v_app_point[:, i])] + # Final drag force vector components + D = D_perpoint * D_direction + + # Total aerodynamic force + Fa = [L[1]+ D[1], L[2]+ D[2], L[3]+ D[3]] + push!(eqs, total_force[1, i] ~ force[1] + Fa[1]) + push!(eqs, total_force[2, i] ~ force[2] + Fa[2]) + push!(eqs, total_force[3, i] ~ force[3] + Fa[3]) + elseif i != 6 + push!(eqs, total_force[:, i] ~ force) + end + push!(eqs, acc[:, i] ~ s.set.g_earth + total_force[:, i] / PointMasses[i]) + eqs2 = vcat(eqs2, reduce(vcat, eqs)) + eqs2 = vcat(eqs2, v_app_point[:, i] ~ s.set.v_wind_tether - vel[:, i]) + end + @named sys = ODESystem(reduce(vcat, Symbolics.scalarize.(eqs2)), t) + simple_sys = structural_simplify(sys) + simple_sys, pos, vel, e_x, e_y, e_z, v_app_point, alpha1p +end +# ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- +# next step function +# ----------------------------- +function next_step!(s::KPS5; dt=s.set.dt) + s.t_0 = s.integrator.t + steptime = @elapsed OrdinaryDiffEqCore.step!(s.integrator, dt, true) + s.iter += 1 + s.integrator.t, steptime +end +# ----------------------------- +# Simulation Function +# ----------------------------- +function simulate(s, logger) + dt = s.set.dt + tol = s.set.tol + tspan = (0.0, dt) + time_range = 0:dt:s.set.sim_time-dt + steps = length(time_range) + iter = 0 + for i in 1:steps + next_step!(s; dt=s.set.dt) + u = s.get_state(s.integrator) + x = u[1][1, :] + y = u[1][2, :] + z = u[1][3, :] + iter += s.iter + sys_state = SysState{points(s)}() + sys_state.X .= x + sys_state.Y .= y + sys_state.Z .= z + println("iter: $iter", " steps: $steps") + log!(logger, sys_state) + println(x[end], ", ", sys_state.X[end]) + end + println("iter: $iter", " steps: $steps") + return nothing +end +# ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- +# Generate Getters +# ----------------------------- +function generate_getters!(s) + sys = s.sys + c = collect + get_state = ModelingToolkit.getu(sys, + [c(sys.pos)] + ) + s.get_state = (integ) -> get_state(integ) + return nothing +end +function play(s, lg) + conn = getconnections(s) + sl = lg.syslog + total_segmentsvector = Vector{Int64}[] + for conn_pair in conn + push!(total_segmentsvector, Int64[conn_pair[1], conn_pair[2]]) + end + # Add tether segments + for i in 0:(s.set.segments-2) + push!(total_segmentsvector, [6+i, 6+i+1]) + end + # Add final connection from last tether point to bridle point + push!(total_segmentsvector, [6+s.set.segments-1, 1]) + for step in 1:length(0:s.set.dt:s.set.sim_time)-1 #-s.set.dt + # Get positions at this time step + x = sl.X[step] + y = sl.Y[step] + z = sl.Z[step] + # Create points array for all points in the system + pointsvector = Vector{Float64}[] + for i in 1:points(s) #FIX THIS! + push!(pointsvector, Float64[x[i], y[i], z[i]]) + end + # Calculate appropriate limits for the plot + x_min, x_max = 0, 40 + z_min, z_max = 0, 60 + t = s.set.dt * (step-1) + # Plot the kite system at this time step + plot2d(pointsvector, total_segmentsvector, t; + zoom = false, + xlim = (x_min, x_max), + ylim = (z_min, z_max) + ) + # Add a small delay to control animation speed + sleep(0.05) + end + nothing +end +function plot_front_view3(lg) + display(plotxy(lg.y, lg.z; + xlabel="pos_y [m]", + ylabel="height [m]", + fig="front_view")) + nothing +end +# ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- +# Running the Main Function +# ----------------------------- +function main() + global lg, s + set = Settings2() + s = KPS5(set=set) + time_range = 0:set.dt:set.sim_time-set.dt + steps = length(time_range) + logger = Logger(points(s), steps) + init_sim!(s) + generate_getters!(s) + simulate(s, logger) + save_log(logger, "tmp") + lg = load_log("tmp") + play(s, lg) +end + +main() +# ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- \ No newline at end of file diff --git a/src/KiteModels.jl b/src/KiteModels.jl index c3d020146..14ac39899 100644 --- a/src/KiteModels.jl +++ b/src/KiteModels.jl @@ -56,7 +56,7 @@ using ADTypes: AutoFiniteDiff using UnPack import ModelingToolkit.SciMLBase: successful_retcode -export KPS3, KPS4, RamAirKite, KVec3, SimFloat, Measurement, PointMassSystem, ProfileLaw, EXP, LOG, EXPLOG # constants and types +export KPS3, KPS4, KPS5, RamAirKite, KVec3, SimFloat, Measurement, PointMassSystem, ProfileLaw, EXP, LOG, EXPLOG # constants and types export calc_set_cl_cd!, copy_examples, copy_bin, update_sys_state! # helper functions export clear!, find_steady_state!, residual! # low level workers export init_sim!, init!, reinit!, next_step!, init_pos_vel, init_pos, model! # high level workers @@ -119,6 +119,7 @@ function __init__() end include("KPS4.jl") # include code, specific for the four point kite model +include("KPS5.jl") # include code, specific for the five point kite model include("ram_air_kite.jl") # include code, specific for the four point 3 line kite model include("mtk_model.jl") include("KPS3.jl") # include code, specific for the one point kite model diff --git a/test/runtests.jl b/test/runtests.jl index a64d776bf..25704197e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -17,16 +17,17 @@ end::Bool cd("..") KiteUtils.set_data_path("") @testset verbose = true "Testing KiteModels..." begin - include("test_orientation.jl") - include("test_kps3.jl") - include("test_kps4.jl") - if build_is_production_build - include("bench3.jl") - include("bench4.jl") - end - if ! haskey(ENV, "NO_MTK") - include("test_ram_air_kite.jl") - end - include("test_helpers.jl") - include("test_inertia_calculation.jl") + # include("test_orientation.jl") + # include("test_kps3.jl") + # include("test_kps4.jl") + include("test_kps5.jl") + # if build_is_production_build + # include("bench3.jl") + # include("bench4.jl") + # end + # if ! haskey(ENV, "NO_MTK") + # include("test_ram_air_kite.jl") + # end + # include("test_helpers.jl") + # include("test_inertia_calculation.jl") end \ No newline at end of file diff --git a/test/test_kps5.jl b/test/test_kps5.jl new file mode 100644 index 000000000..7525e0d35 --- /dev/null +++ b/test/test_kps5.jl @@ -0,0 +1,91 @@ +using Pkg +if ! ("PackageCompiler" ∈ keys(Pkg.project().dependencies)) + using TestEnv; TestEnv.activate() + Pkg.update() +end +using Test, BenchmarkTools, StaticArrays, LinearAlgebra, KiteUtils +using KiteModels, KitePodModels + +set_data_path(joinpath(dirname(dirname(pathof(KiteModels))), "data")) +set = deepcopy(load_settings("system.yaml")) +kcu::KCU = KCU(set) +kps5::KPS5 = KPS5(kcu) + +pos, vel = nothing, nothing +poss, vels = nothing, nothing + +@testset verbose = true "KPS5 tests...." begin + +function set_defaults() + KiteModels.clear!(kps5) + kps5.set.l_tether = 150.0 + kps5.set.elevation = 60.0 + kps5.set.area = 20.0 + kps5.set.rel_side_area = 50.0 + kps5.set.v_wind = 8.0 + kps5.set.mass = 11.4 + kps5.set.damping = 2 * 473.0 + kps5.set.alpha = 1.0/7 + kps5.set.c_s = 0.6 + kps5.set.kcu_diameter = 0 +end + +function init_392() + KiteModels.clear!(kps5) + kps5.set.l_tether = 392.0 + kps5.set.elevation = 70.0 + kps5.set.area = 10.0 + kps5.set.rel_side_area = 50.0 + kps5.set.v_wind = 9.1 + kps5.set.mass = 6.2 + kps5.set.c_s = 0.6 +end + +function init_150() + KiteModels.clear!(kps5) + kps5.set.l_tether = 150.0 + kps5.set.elevation = 70.0 + kps5.set.area = 10.18 + kps5.set.rel_side_area = 30.6 + kps5.set.v_wind = 9.1 + kps5.set.mass = 6.21 + kps5.set.c_s = 0.6 + kps5.set.damping = 473.0 # unit damping + kps5.set.c_spring = 614600.0 # unit spring coefficent + kps5.set.width = 4.9622 +end + +function init3() + kps5.set.alpha = 0.08163 + KiteModels.clear!(kps5) + kps5.set.l_tether = 150.0 # - kps5.set.height_k - kps5.set.h_bridle + kps5.set.area = 10.18 + kps5.set.rel_side_area = 30.6 + kps5.set.mass = 6.21 + kps5.set.c_s = 0.6 + kps5.set.damping = 473.0 # unit damping coefficient + kps5.set.c_spring = 614600.0 # unit spring coefficent + kps5.set.width = 4.9622 + kps5.set.elevation = 70.7 + kps5.set.profile_law = Int(EXPLOG) + pos, vel = KiteModels.init_pos_vel(kps5) + posd = copy(vel) + veld = zero(vel) + height = 134.14733504839947 + kps5.v_wind .= kps5.v_wind_gnd * calc_wind_factor(kps5.am, height) + kps5.stiffness_factor = 1.0 + KiteModels.init_springs!(kps5) + return pos, vel, posd, veld +end + +set_defaults() + +@testset "calc_rho " begin + @test isapprox(calc_rho(kps5.am, 0.0), 1.225, atol=1e-5) + @test isapprox(calc_rho(kps5.am, 100.0), 1.210756, atol=1e-5) + kps5_= KPS5(kcu) + @test kps5_ isa KPS5 +end + +end +nothing \ No newline at end of file