Merge pull request #55 from ShipMMG/feature/wind_force

Add wind force and moment

.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

src/ShipMMG.jl

@@ -13,6 +13,7 @@ export calc_position,
     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
 
@@ -27,11 +28,13 @@ export kt_simulate,
 include("mmg.jl")
 export Mmg3DofBasicParams,
     Mmg3DofManeuveringParams,
+    Mmg3DofWindForceMomentParams,
     mmg_3dof_simulate,
     mmg_3dof_model!,
     mmg_3dof_zigzag_test,
     estimate_mmg_approx_lsm,
     estimate_mmg_approx_lsm_time_window_sampling,
-    create_model_for_mcmc_sample_mmg
+    create_model_for_mcmc_sample_mmg,
+    wind_force_and_moment_coefficients
 
 end

src/data.jl

@@ -10,7 +10,7 @@
     n_p::Tn_p
 end
 
-function get_KVLCC2_L7_basic_params(ρ = 1025.0)
+function get_KVLCC2_L7_basic_params(ρ=1025.0)
     L_pp = 7.00  # 船長Lpp[m]
     B = 1.27  # 船幅[m]
     d = 0.46  # 喫水[m]
@@ -144,19 +144,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 +286,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 +301,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

src/kt.jl

@@ -85,9 +85,8 @@ function kt_simulate(
     Ψ0=0.0,
     algorithm=Tsit5(),
     reltol=1e-8,
-    abstol=1e-8
+    abstol=1e-8,
 )
-
     spl_δ = Spline1D(time_list, δ_list)
 
     X0 = [u0; v0; r0; x0; y0; Ψ0; δ_list[1]]
@@ -168,7 +167,7 @@ function kt_zigzag_test(
     δ_rad_rate=10.0 * π / 180,
     algorithm=Tsit5(),
     reltol=1e-8,
-    abstol=1e-8
+    abstol=1e-8,
 )
     target_Ψ_rad_deviation = abs(target_Ψ_rad_deviation)
 
@@ -205,7 +204,7 @@ function kt_zigzag_test(
             y0 = final_y_list[start_index-1]
         end
 
-        for i = (start_index+1):length(time_list)
+        for i in (start_index+1):length(time_list)
             Δt = time_list[i] - time_list[i-1]
             if target_δ_rad > 0
                 δ = δ_list[i-start_index] + δ_rad_rate * Δt
@@ -306,7 +305,7 @@ function estimate_kt_lsm_time_window_sampling(data::ShipData, window_size::Int)
     n_samples = length(time_vec) - window_size
     K_samples = zeros(n_samples)
     T_samples = zeros(n_samples)
-    for i = 1:n_samples
+    for i in 1:n_samples
         K_samples[i], T_samples[i] =
             estimate_kt_lsm(r[i:i+window_size], dr[i:i+window_size], δ[i:i+window_size])
     end
@@ -317,7 +316,7 @@ function create_model_for_mcmc_sample_kt(
     data::ShipData;
     σ_r_prior_dist::Distribution=Chi(5),
     K_prior_dist::Distribution=Uniform(0.01, 10.0),
-    T_prior_dist::Distribution=truncated(Normal(100.0, 50.0), 10.0, 200.0)
+    T_prior_dist::Distribution=truncated(Normal(100.0, 50.0), 10.0, 200.0),
 )
     time_obs = data.time
     r_obs = data.r
@@ -347,7 +346,7 @@ function create_model_for_mcmc_sample_kt(
         prob = remake(prob1, p=p)
         sol = solve(prob, Tsit5())
         predicted = sol(time_obs)
-        for i = 1:length(predicted)
+        for i in 1:length(predicted)
             r_obs[i] ~ Normal(predicted[i][1], σ_r) # index number of r is 1
         end
     end