Add kps4 3 line example (#61)

* add 3 line example


---------

Co-authored-by: Uwe Fechner <[email protected]>
Co-authored-by: Uwe Fechner <[email protected]>

.gitignore

@@ -28,3 +28,4 @@ data/sim_log.bin
 docs/src/winch.png
 bin/kps-image-1.10-main.so.bak
 bin/kps-image-1.10-1-Bart-1-main.so
+bin/kps-image-1.10-kps4_3l_example.so

data/3l_settings.yaml

@@ -0,0 +1,117 @@
+system:
+    log_file: "data/log_8700W_8ms" # filename without extension  [replay only]
+                                   #   use / as path delimiter, even on Windows 
+    time_lapse:   1.0              # relative replay speed
+    sim_time:   100.0              # 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
+
+initial:
+    l_tether: 50.0        # initial tether length       [m]
+    elevation: 70.8        # initial elevation angle   [deg]
+    v_reel_out: 0.0        # initial reel out speed    [m/s]
+    depower:   25.0        # initial depower settings    [%]
+
+solver:
+    abs_tol: 0.006        # absolute tolerance of the DAE solver [m, m/s]
+    rel_tol: 0.01         # relative tolerance of the DAE solver [-]
+    solver: "DFBDF"        # DAE solver, IDA or DImplicitEuler, DFBDF
+    linear_solver: "GMRES" # can be GMRES or Dense or LapackDense (only for IDA)
+    max_order: 4           # maximal order, usually between 3 and 5
+    max_iter:  10000         # 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.83    # max. steering angle of the side planes for four point model [degrees]
+
+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_3L" # name of the kite model to use (KPS3 or KPS4)
+    version: 2             # version of the model to use
+    mass:  0.9             # kite mass [kg]
+    area: 10.18            # projected kite area         [m²]
+    rel_side_area: 30.6    # relative side area           [%]
+    height: 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: 
+    alpha_zero:   10.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
+
+kps4_3l:
+    radius: 2.0
+    bridle_center_distance: 4.0
+    middle_length: 1.5
+    tip_length: 0.62
+    min_steering_line_distance: 1.0
+    width:         4.1     # width of the kite                          [m]
+    aero_surfaces: 3      # number of aerodynamic surfaces in 3 line model
+
+kcu:
+    kcu_mass: 8.4                # mass of the kite control unit   [kg]
+    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
+
+bridle:
+    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)
+
+tether:
+    d_tether:  1           # 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:
+    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]
+
+environment:
+    v_wind: 15.51             # wind speed at reference height          [m/s]
+    v_wind_ref: [15.51, 0.0]  # wind speed vector at reference height   [m/s]
+    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 Cabau, 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]
+                             # Cabau: 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]               

data/system.yaml

@@ -1,3 +1,3 @@
 system:
-    sim_settings: "settings.yaml"  # simulator settings
+    sim_settings: "3l_settings.yaml"  # simulator settings
 

examples/reel_out_4p_3l.jl

@@ -0,0 +1,105 @@
+using Printf
+using KiteModels, KitePodModels, KiteUtils
+
+set = deepcopy(se())
+
+# the following values can be changed to match your interest
+dt = 0.05
+set.solver="DFBDF" # IDA or DFBDF
+STEPS = 100
+PLOT = true
+FRONT_VIEW = false
+ZOOM = false
+PRINT = true
+STATISTIC = true
+ALPHA_ZERO = 8.8 
+# end of user parameter section #
+
+set.alpha_zero = ALPHA_ZERO
+set.version = 2
+
+kcu::KCU = KCU(set)
+kps4_3l::KPS4_3L = KPS4_3L(kcu)
+
+# if PLOT
+    using Pkg
+    if ! ("ControlPlots" ∈ keys(Pkg.project().dependencies))
+        using TestEnv; TestEnv.activate()
+    end
+    using ControlPlots
+# end
+
+v_time = zeros(STEPS)
+v_speed = zeros(STEPS)
+v_force = zeros(STEPS)
+
+function simulate(integrator, steps, plot=false)
+    start = integrator.p.iter
+    for i in 1:steps
+        if PRINT
+            lift, drag = KiteModels.lift_drag(kps4_3l)
+            @printf "%.2f: " round(integrator.t, digits=2)
+            println("lift, drag  [N]: $(round(lift, digits=2)), $(round(drag, digits=2))")
+        end
+        acc = 0.0
+        if kps4_3l.t_0 > 1.0
+            acc = 0.1
+        end
+        set_speeds = kps4_3l.set_speeds.+acc*dt
+        v_time[i] = kps4_3l.t_0
+        v_speed[i] = kps4_3l.reel_out_speeds[1]
+        v_force[i] = winch_force(kps4_3l)[1]
+        KiteModels.next_step!(kps4_3l, integrator; set_values=set_speeds, dt=dt)
+
+        if plot
+            reltime = i*dt-dt
+            if mod(i, 5) == 1
+                plot2d(kps4_3l.pos, reltime; zoom=ZOOM, front=FRONT_VIEW, 
+                                        segments=set.segments, fig="side_view")            
+            end
+        end
+    end
+    (integrator.p.iter - start) / steps
+end
+
+integrator = KiteModels.init_sim!(kps4_3l, stiffness_factor=0.3, prn=STATISTIC)
+kps4_3l.set_speeds = [0.0, 0.0, 0.0]
+
+if PLOT
+    av_steps = simulate(integrator, STEPS, true)
+else
+    println("\nStarting simulation...")
+    simulate(integrator, 100)
+    runtime = @elapsed av_steps = simulate(integrator, STEPS-100)
+    println("\nSolver: $(set.solver)")
+    println("Total simulation time: $(round(runtime, digits=3)) s")
+    speed = (STEPS-100) / runtime * dt
+    println("Simulation speed: $(round(speed, digits=2)) times realtime.")
+end
+lift, drag = KiteModels.lift_drag(kps4_3l)
+println("lift, drag  [N]: $(round(lift, digits=2)), $(round(drag, digits=2))")
+println("Average number of callbacks per time step: $av_steps")
+
+if PLOT
+    p = plotx(v_time, v_speed, v_force; ylabels=["v_reelout  [m/s]","tether_force [N]"], fig="winch")
+    display(p)
+end
+# savefig("docs/src/reelout_force_4p.png")
+
+# Solver: DFBDF, reltol=0.000001
+# Total simulation time: 1.049 s
+# Simulation speed: 23.83 times realtime.
+# lift, drag  [N]: 545.24, 102.55
+# Average number of callbacks per time step: 625.604
+
+# Solver: DFBDF, reltol=0.001
+# Total simulation time: 0.178 s
+# Simulation speed: 140.08 times realtime.
+# lift, drag  [N]: 558.69, 105.3
+# Average number of callbacks per time step: 114.64
+
+# Solver: IDA
+# Total simulation time: 1.385 s
+# Simulation speed: 18.05 times realtime.
+# lift, drag  [N]: 543.58, 102.19
+# Average number of callbacks per time step: 756.074