Merge pull request #62 from ShipMMG/dev

Considering wind force and moment in MMG 3DOF model

Author: taiga4112

.JuliaFormatter.toml

@@ -0,0 +1,7 @@
+remove_extra_newlines=true
+join_lines_based_on_source=true
+whitespace_in_kwargs=false
+short_to_long_function_def=true
+always_for_in=true
+verbose=true
+margin=88

.github/workflows/Test.yml

@@ -0,0 +1,30 @@
+name: Run tests
+
+on: 
+  push:
+    branches:
+      - main
+  pull_request:
+    branches:
+      - main
+      - dev
+
+jobs:
+  test:
+    runs-on: ${{ matrix.os }}
+    strategy:
+      matrix:
+        julia-version: ['1']
+        julia-arch: [x64]
+        os: [ubuntu-latest]
+
+    steps:
+      - uses: actions/checkout@v2
+      - uses: julia-actions/[email protected]
+        with:
+          version: ${{ matrix.julia-version }}
+          arch: ${{ matrix.julia-arch }}
+      - uses: julia-actions/julia-buildpkg@v1
+      - uses: julia-actions/julia-runtest@v1
+        with:
+          annotate: true

Project.toml

@@ -1,17 +1,18 @@
 name = "ShipMMG"
 uuid = "37f2b0bf-0c13-4883-8808-e75eb56597e7"
 authors = ["Taiga MITSUYUKI <[email protected]>"]
-version = "0.0.6"
+version = "0.0.7"
 
 [deps]
 AdvancedPS = "576499cb-2369-40b2-a588-c64705576edc"
 Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+KernelDensity = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b"
 Libtask = "6f1fad26-d15e-5dc8-ae53-837a1d7b8c9f"
 ParameterizedFunctions = "65888b18-ceab-5e60-b2b9-181511a3b968"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Turing = "fce5fe82-541a-59a6-adf8-730c64b5f9a0"
 
 [compat]
@@ -22,7 +23,6 @@ Distributions = "^0.25"
 Libtask = "^0.8"
 ParameterizedFunctions = "^5.12"
 Parameters = "0.12"
-Plots = "^1"
 Turing = "^0.24"
 julia = "^1.6"
 

src/ShipMMG.jl

@@ -3,17 +3,19 @@ module ShipMMG
 using DifferentialEquations
 using ParameterizedFunctions
 using Dierckx
-using Plots
 using Parameters
 using Distributions
 using Turing
+using ForwardDiff
 
 include("data.jl")
 export calc_position,
     ShipData,
+    EnvironmentData,
     get_KVLCC2_L7_basic_params,
     get_KVLCC2_L7_maneuvering_params,
     get_KVLCC2_L7_params,
+    get_example_ship_wind_force_moment_params,
     nuts_sampling_single_thread,
     nuts_sampling_multi_threads
 
@@ -35,7 +37,17 @@ export Mmg3DofBasicParams,
     estimate_mmg_approx_lsm_time_window_sampling,
     create_model_for_mcmc_sample_mmg
 
-include("draw.jl")
-export draw_gif_result, calc_position
+include("mmg_wind.jl")
+export Mmg3DofBasicParams,
+    Mmg3DofManeuveringParams,
+    Mmg3DofWindForceMomentParams,
+    apparent_wind_speed_and_angle,
+    mmg_3dof_simulate,
+    mmg_3dof_wind_model!,
+    mmg_3dof_zigzag_test,
+    estimate_mmg_approx_lsm,
+    estimate_mmg_approx_lsm_time_window_sampling,
+    create_model_for_mcmc_sample_mmg,
+    wind_force_and_moment_coefficients
 
 end

src/data.jl

@@ -10,7 +10,13 @@
     n_p::Tn_p
 end
 
-function get_KVLCC2_L7_basic_params(ρ = 1025.0)
+@with_kw struct EnvironmentData{Tt,TU_W,Tψ_W}
+    time::Tt
+    U_W::TU_W
+    ψ_W::Tψ_W
+end
+
+function get_KVLCC2_L7_basic_params(ρ=1025.0)
     L_pp = 7.00  # 船長Lpp[m]
     B = 1.27  # 船幅[m]
     d = 0.46  # 喫水[m]
@@ -144,19 +150,140 @@ function get_KVLCC2_L7_params()
     basic_params, maneuvering_params
 end
 
-function calc_position(time_vec, u_vec, v_vec, r_vec; x0 = 0.0, y0 = 0.0, ψ0 = 0.0)
+"""
+    wind_force_and_moment_coefficients(ψ_A,p)
+
+compute the parameter C_X,C_Y,C_N from Fujiwara's estimation method (1994).
+
+# Arguments
+-`ψ_A`: Wind attack of angle [rad]
+- L_pp
+- B
+- A_OD
+- A_F
+- A_L
+- H_BR
+- H_C
+- C
+"""
+function wind_force_and_moment_coefficients(
+    ψ_A,
+    L_pp,
+    B,
+    A_OD,
+    A_F,
+    A_L,
+    H_BR,
+    H_C,
+    C,
+)
+    #C_LF1の場合で調整
+    C_CF = 0.404 + 0.368 * A_F / (B * H_BR) + 0.902 * H_BR / L_pp
+
+    if deg2rad(0) <= ψ_A <= deg2rad(90)
+        C_LF = -0.992 + 0.507 * A_L / (L_pp * B) + 1.162 * C / L_pp
+        C_XLI = 0.458 + 3.245 * A_L / (L_pp * H_BR) - 2.313 * A_F / (B * H_BR)
+        C_ALF = -0.585 - 0.906 * A_OD / A_L + 3.239 * B / L_pp
+        C_YLI = pi * A_L / L_pp^2 + 0.116 + 3.345 * A_F / (L_pp * B)
+
+    elseif deg2rad(90) < ψ_A <= deg2rad(180)
+        C_LF =
+            0.018 - 5.091 * B / L_pp + 10.367 * H_C / L_pp - 3.011 * A_OD / L_pp^2 -
+            0.341 * A_F / B^2
+        C_XLI =
+            -1.901 + 12.727 * A_L / (L_pp * H_BR) + 24.407 * A_F / A_L -
+            40.310 * B / L_pp - 0.341 * A_F / (B * H_BR)
+        C_ALF = -0.314 - 1.117 * A_OD / A_L
+        C_YLI = pi * A_L / L_pp^2 + 0.446 + 2.192 * A_F / L_pp^2
+
+    elseif deg2rad(180) < ψ_A <= deg2rad(270)
+        C_LF =
+            0.018 - 5.091 * B / L_pp + 10.367 * H_C / L_pp - 3.011 * A_OD / L_pp^2 -
+            0.341 * A_F / B^2
+        C_XLI =
+            -1.901 + 12.727 * A_L / (L_pp * H_BR) + 24.407 * A_F / A_L -
+            40.310 * B / L_pp - 0.341 * A_F / (B * H_BR)
+        C_ALF = -(-0.314 - 1.117 * A_OD / A_L)
+        C_YLI = -(pi * A_L / L_pp^2 + 0.446 + 2.192 * A_F / L_pp^2)
+
+    elseif deg2rad(270) < ψ_A <= deg2rad(360)
+        C_LF = -0.992 + 0.507 * A_L / (L_pp * B) + 1.162 * C / L_pp
+        C_XLI = 0.458 + 3.245 * A_L / (L_pp * H_BR) - 2.313 * A_F / (B * H_BR)
+        C_ALF = -(-0.585 - 0.906 * A_OD / A_L + 3.239 * B / L_pp)
+        C_YLI = -(pi * A_L / L_pp^2 + 0.116 + 3.345 * A_F / (L_pp * B))
+    end
+
+    C_X =
+        C_LF * cos(ψ_A) +
+        C_XLI * (sin(ψ_A) - sin(ψ_A) * cos(ψ_A)^2 / 2) * sin(ψ_A) * cos(ψ_A) +
+        C_ALF * sin(ψ_A) * cos(ψ_A)^3
+    C_Y =
+        C_CF * sin(ψ_A)^2 +
+        C_YLI * (cos(ψ_A) + sin(ψ_A)^2 * cos(ψ_A) / 2) * sin(ψ_A) * cos(ψ_A)
+    C_N = C_Y * (0.297 * C / L_pp - 0.149 * (ψ_A - deg2rad(90)))
+
+    C_X, C_Y, C_N
+end
+
+function get_example_ship_wind_force_moment_params()
+    L_pp = 7.00  # 船長Lpp[m]
+    B = 1.27  # 船幅[m]
+    d = 0.46  # 喫水[m]
+    D = 0.6563 # 深さ[m]
+    A_OD = 0.65 # デッキ上の構造物の側面投影面積[m^2]
+    H_BR = 0.85 # 喫水からブリッジ主要構造物の最高位[m]
+    H_C = 0.235 # 喫水から側面積中心までの高さ[m]
+    C = 0.0 # 船体中心から側面積中心までの前後方向座標(船首方向を正)[m]
+
+    A_OD = A_OD # デッキ上の構造物の側面投影面積[m^2]
+    A_F = (D - d) * B  # 船体の正面投影面積[m^2]
+    A_L = (D - d) * L_pp # 船体の側面投影面積[m^2]
+    H_BR = H_BR # 喫水からブリッジ主要構造物の最高位[m]
+    H_C = H_C # 喫水から側面積中心までの高さ[m]
+    C = C # 船体中心から側面積中心までの前後方向座標[m]
+
+    ψ_A_vec = deg2rad.(collect(0:10:360))
+    C_X_vec = Array{Float64}(undef, length(ψ_A_vec))
+    C_Y_vec = Array{Float64}(undef, length(ψ_A_vec))
+    C_N_vec = Array{Float64}(undef, length(ψ_A_vec))
+    for (index, ψ_A) in enumerate(ψ_A_vec)
+        C_X, C_Y, C_N = wind_force_and_moment_coefficients(
+            ψ_A,
+            L_pp,
+            B,
+            A_OD,
+            A_F,
+            A_L,
+            H_BR,
+            H_C,
+            C,
+        )
+        C_X_vec[index] = C_X
+        C_Y_vec[index] = C_Y
+        C_N_vec[index] = C_N
+    end
+    spl_C_X = Spline1D(ψ_A_vec, C_X_vec)
+    spl_C_Y = Spline1D(ψ_A_vec, C_Y_vec)
+    spl_C_N = Spline1D(ψ_A_vec, C_N_vec)
+
+    Mmg3DofWindForceMomentParams(A_F, A_L, spl_C_X, spl_C_Y, spl_C_N)
+end
+
+function calc_position(time_vec, u_vec, v_vec, r_vec; x0=0.0, y0=0.0, ψ0=0.0)
     dim = length(time_vec)
     x_vec = zeros(Float64, dim)
     y_vec = zeros(Float64, dim)
     ψ_vec = zeros(Float64, dim)
     x_vec[1] = x0
     y_vec[1] = y0
     ψ_vec[1] = ψ0
-    for i = 2:dim
+    for i in 2:dim
         Δt = time_vec[i] - time_vec[i-1]
         ψ_vec[i] = ψ_vec[i-1] + r_vec[i] * Δt
-        x_vec[i] = x_vec[i-1] + (u_vec[i] * cos(ψ_vec[i]) - v_vec[i] * sin(ψ_vec[i])) * Δt
-        y_vec[i] = y_vec[i-1] + (u_vec[i] * sin(ψ_vec[i]) + v_vec[i] * cos(ψ_vec[i])) * Δt
+        x_vec[i] =
+            x_vec[i-1] + (u_vec[i] * cos(ψ_vec[i]) - v_vec[i] * sin(ψ_vec[i])) * Δt
+        y_vec[i] =
+            y_vec[i-1] + (u_vec[i] * sin(ψ_vec[i]) + v_vec[i] * cos(ψ_vec[i])) * Δt
     end
     x_vec, y_vec, ψ_vec
 end
@@ -165,12 +292,12 @@ function nuts_sampling_single_thread(
     model,
     n_samples::Int,
     n_chains::Int;
-    target_acceptance::Float64 = 0.65,
-    progress = true,
+    target_acceptance::Float64=0.65,
+    progress=true,
 )
     sampler = NUTS(target_acceptance)
     mapreduce(
-        c -> sample(model, sampler, n_samples, progress = progress),
+        c -> sample(model, sampler, n_samples, progress=progress),
         chainscat,
         1:n_chains,
     )
@@ -180,9 +307,9 @@ function nuts_sampling_multi_threads(
     model,
     n_samples::Int,
     n_chains::Int;
-    target_acceptance::Float64 = 0.65,
-    progress = false,
+    target_acceptance::Float64=0.65,
+    progress=false,
 )
     sampler = NUTS(target_acceptance)
-    sample(model, sampler, MCMCThreads(), n_samples, n_chains, progress = progress)
+    sample(model, sampler, MCMCThreads(), n_samples, n_chains, progress=progress)
 end