Add a 3 line kite model. (#60)

* serialization of initial integrator * even faster by serializing only once * stable data type * update documentation * add check for model * less memory * remove float64 * move to simfloat * two point aero model * working two point aero model * xournal notes with formulas * implement new aero model * add comment to explain point numbers * implement kps4_3l functions * plot lift at different wing positions * add kps4 3 line settings * plot * bug free init * implement 3 lines * runs succesfully, not tested * find steady state runs but doesnt diverge * converging but not solving steady state * not solving steady state * notes * updated notes * change page order * change page order * Fix precompilation problem * add file to .gitignore * different tether length * remove Plots from Project * use package directly * correct bridle center distance * shorter name * update bridle center distance * add kps4_3l settings * solve forward slash in branch name * succesful but very slow step * working simulation, just slow * still really slow * benchmark test showing big mem alloc * add 3l test * a lot less memory in calc aero forces * magic memory alloc * remove a lot of heap alloc * working but slow * a lot less allocations * zero heap alloc in calc_aero_forces * no heap alloc * working no heap alloc version * fix E distance * flying makes sense and is fast * add kps4_3l to precompile workload * add orient calc func * export function * remove unnecesary print * remove controlplots dependency * remove unused import * remove unnecesary imports * calc azimuth elevation * need to test int history * fix keyword error * remove precompiling for 3 and 4 * fix issues after merge * fix benchmark mem alloc * add test cases for 3 line model * one 4p model test failing * succesful tests * fix precompile * fix typo --------- Co-authored-by: Uwe Fechner <[email protected]> Co-authored-by: Uwe Fechner <[email protected]>
ufechner7 · Aug 9, 2024 · 52a15a9 · 52a15a9
1 parent e0921f9
commit 52a15a9
Show file tree

Hide file tree

Showing 20 changed files with 1,956 additions and 330 deletions.
diff --git a/.gitignore b/.gitignore
@@ -1,3 +1,4 @@
+*.xopp
 /Manifest.toml
 /Manifest.toml.1.7
 bin/kps-image-1.8-initial.so
@@ -21,6 +22,9 @@ bin/kps-image-1.10-3line.so
 bin/kps-image-1.10-controlplots.so
 bin/kps-image-1.11-main.so
 bin/kps-image-1.10-torque.so
+bin/kps-image-1.10-feat-KPS4_3L.so
+data/.steady_state_history.bin
+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
diff --git a/Project.toml b/Project.toml
@@ -5,6 +5,7 @@ version = "0.6.3"
 
 [deps]
 AtmosphericModels = "c59cac55-771d-4f45-b14d-1c681463a295"
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
@@ -17,8 +18,10 @@ Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc"
+Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
+TestEnv = "1e6cf692-eddd-4d53-88a5-2d735e33781b"
 WinchModels = "7dcfa46b-7979-4771-bbf4-0aee0da42e1f"
 
 [compat]

diff --git a/bin/create_sys_image b/bin/create_sys_image
@@ -22,7 +22,7 @@ if [[ $julia_major == "1.1" ]]; then
     julia_major=${julia_version:0:4} 
 fi
 if [ -d .git ] || git rev-parse --git-dir > /dev/null 2>&1 ; then
-    branch=$(git rev-parse --abbrev-ref HEAD)
+    branch=$(git rev-parse --abbrev-ref HEAD | sed 's/\//-/g')
 else
     branch=""
 fi

diff --git a/data/3l_settings.yaml b/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]               
diff --git a/data/settings.yaml b/data/settings.yaml
@@ -12,7 +12,7 @@ system:
 
 initial:
     l_tether: 150.0        # initial tether length       [m]
-    elevation: 70.7        # initial elevation angle   [deg]
+    elevation: 70.8        # initial elevation angle   [deg]
     v_reel_out: 0.0        # initial reel out speed    [m/s]
     depower:   25.0        # initial depower settings    [%]
 

diff --git a/data/system.yaml b/data/system.yaml
@@ -1,3 +1,3 @@
 system:
-    sim_settings: "settings.yaml"  # simulator settings
+    sim_settings: "3l_settings.yaml"  # simulator settings
 
diff --git a/docs/3-line_kite_model.pdf b/docs/3-line_kite_model.pdf
diff --git a/docs/3-line_kite_model.xopp b/docs/3-line_kite_model.xopp
diff --git a/src/KPS3.jl b/src/KPS3.jl
@@ -482,7 +482,7 @@ end
 function init(s::KPS3, X=zeros(2 * (s.set.segments)); old=false, delta = 0.0)
     res1_, res2_ = init_inner(s, X; old=old, delta = delta)
     res1, res2  = vcat(reduce(vcat, res1_), [s.l_tether, 0]), vcat(reduce(vcat, res2_),[0,0])
-    MVector{6*(s.set.segments)+2, Float64}(res1), MVector{6*(s.set.segments)+2, Float64}(res2)
+    MVector{6*(s.set.segments)+2, SimFloat}(res1), MVector{6*(s.set.segments)+2, SimFloat}(res2)
 end
 
 """

diff --git a/src/KPS4.jl b/src/KPS4.jl
@@ -44,6 +44,8 @@ const SPRINGS_INPUT = [0.    1.  150.
                        5.    2.   -1. # s8, p11, p8
                        2.    3.   -1.] # s9, p8, p9
 
+# KCU = p7, A = p8, B = p9, C = p10, D = p11
+
 # struct, defining the phyical parameters of one spring
 @with_kw struct Spring{I, S}
     p1::I = 1         # number of the first point
@@ -53,7 +55,7 @@ const SPRINGS_INPUT = [0.    1.  150.
     damping::S  = 0.1 # damping coefficent [Ns/m]
 end
 
-const SP = Spring{Int16, Float64}
+const SP = Spring{Int16, SimFloat}
 const KITE_PARTICLES = 4
 const KITE_SPRINGS = 9
 const KITE_ANGLE = 3.83 # angle between the kite and the last tether segment due to the mass of the control pod
@@ -164,10 +166,10 @@ $(TYPEDFIELDS)
     z::T =                 zeros(S, 3)
 end
 
-@inline @inbounds function norm(vec::SVector{3, Float64})
+@inline @inbounds function norm(vec::SVector{3, SimFloat})
     sqrt(vec[1]*vec[1]+vec[2]*vec[2]+vec[3]*vec[3])
 end
-@inline @inbounds function norm(vec::MVector{3, Float64})
+@inline @inbounds function norm(vec::MVector{3, SimFloat})
     sqrt(vec[1]*vec[1]+vec[2]*vec[2]+vec[3]*vec[3])
 end
 
@@ -223,7 +225,7 @@ Calculate the drag force of the tether segment, defined by the parameters pos1,
 and distribute it equally on the two particles, that are attached to the segment.
 The result is stored in the array s.forces. 
 """
-@inline function calc_particle_forces!(s, pos1, pos2, vel1, vel2, spring, segments, d_tether, rho, i)
+@inline function calc_particle_forces!(s::KPS4, pos1, pos2, vel1, vel2, spring, segments, d_tether, rho, i)
     l_0 = spring.length # Unstressed length
     k = spring.c_spring * s.stiffness_factor  # Spring constant
     c = spring.damping                        # Damping coefficient    
@@ -424,7 +426,7 @@ function residual!(res, yd, y::MVector{S, SimFloat}, s::KPS4, time) where S
     posd1, veld1 = partd[:,:,1], partd[:,:,2]
     posd = SVector{div(T,6)+1}(if i==1 SVector(0.0,0,0) else SVector(posd1[:,i-1]) end for i in 1:div(T,6)+1)
     veld = SVector{div(T,6)+1}(if i==1 SVector(0.0,0,0) else SVector(veld1[:,i-1]) end for i in 1:div(T,6)+1)
-    @assert ! isnan(pos[2][3])
+    @assert isfinite(pos[2][3])
 
     # core calculations
     s.l_tether = length