Foil sjekte

foilsjekte.jl

using Plots
using Optim

MS_KNOTS = 1.943844
KG_LBS = 2.204623
KW_HP = 1.341022
M_FEET = 3.28084

# Dimensjoner på båten (i meter, kg etc)
# https://marex.no/wp-content/uploads/2016/12/Duckie-Catalog.pdf
l = 6.5
b = 2.46
m = 1400
m_pass = 6 * 80
m_tot = m + m_pass
g = 9.81

# hastigheter
vv = 0.5:0.5:10.0 # ca 20 knop

# foil parametre
# http://airfoiltools.com/airfoil/details?airfoil=clarky-il
rho = 1000
wingspan = b
n_wings = 3.0
aspect = 15
area = wingspan ^2 / aspect * n_wings
c_l = 1.2 # Reynolds number approx 200.000 .. 2.000.000
g_e = 60 # at reynolds 200.000
c_d = c_l / 60.0 # skin friction?


POLAR_FILE = [
  -3.500   0.0091   0.01659   0.00877  -0.0916   0.9112   0.0772
  -3.250   0.0342   0.01609   0.00829  -0.0910   0.9027   0.0811
  -3.000   0.0668   0.01553   0.00768  -0.0915   0.8975   0.0871
  -2.750   0.0940   0.01485   0.00710  -0.0912   0.8900   0.0985
  -2.500   0.1221   0.01403   0.00642  -0.0910   0.8825   0.1304
  -2.250   0.1479   0.01328   0.00605  -0.0905   0.8741   0.2081
  -2.000   0.1751   0.01271   0.00578  -0.0901   0.8662   0.2941
  -1.750   0.1998   0.01226   0.00560  -0.0894   0.8568   0.3729
  -1.500   0.2249   0.01156   0.00536  -0.0885   0.8491   0.4994
  -1.250   0.2429   0.01087   0.00534  -0.0860   0.8383   0.6774
  -1.000   0.2709   0.01036   0.00534  -0.0841   0.8309   0.8768
  -0.750   0.3315   0.01030   0.00523  -0.0899   0.8231   0.9730
  -0.500   0.3953   0.01014   0.00489  -0.0970   0.8170   1.0000
  -0.250   0.4182   0.01015   0.00480  -0.0960   0.8055   1.0000
   0.000   0.4424   0.01016   0.00470  -0.0952   0.7954   1.0000
   0.250   0.4680   0.01015   0.00455  -0.0945   0.7864   1.0000
   0.500   0.4915   0.01020   0.00453  -0.0935   0.7750   1.0000
   0.750   0.5161   0.01024   0.00447  -0.0927   0.7644   1.0000
   1.000   0.5419   0.01026   0.00437  -0.0920   0.7544   1.0000
   1.250   0.5656   0.01032   0.00437  -0.0911   0.7423   1.0000
   1.500   0.5899   0.01039   0.00437  -0.0902   0.7308   1.0000
   1.750   0.6149   0.01046   0.00436  -0.0895   0.7196   1.0000
   2.000   0.6402   0.01053   0.00433  -0.0887   0.7080   1.0000
   2.250   0.6641   0.01061   0.00436  -0.0877   0.6945   1.0000
   2.500   0.6883   0.01071   0.00440  -0.0868   0.6807   1.0000
   2.750   0.7126   0.01081   0.00444  -0.0859   0.6666   1.0000
   3.000   0.7370   0.01092   0.00449  -0.0851   0.6519   1.0000
   3.250   0.7612   0.01104   0.00456  -0.0842   0.6360   1.0000
   3.500   0.7852   0.01117   0.00462  -0.0832   0.6187   1.0000
   3.750   0.8091   0.01133   0.00468  -0.0823   0.6003   1.0000
   4.000   0.8324   0.01152   0.00480  -0.0812   0.5801   1.0000
   4.250   0.8553   0.01174   0.00494  -0.0801   0.5580   1.0000
   4.500   0.8781   0.01201   0.00508  -0.0790   0.5358   1.0000
   4.750   0.9005   0.01231   0.00530  -0.0779   0.5124   1.0000
   5.000   0.9227   0.01266   0.00551  -0.0768   0.4911   1.0000
   5.250   0.9450   0.01300   0.00576  -0.0758   0.4700   1.0000
   5.500   0.9671   0.01335   0.00602  -0.0747   0.4502   1.0000
   5.750   0.9892   0.01371   0.00630  -0.0737   0.4321   1.0000
   6.000   1.0116   0.01408   0.00660  -0.0727   0.4164   1.0000
   6.250   1.0345   0.01444   0.00693  -0.0719   0.4033   1.0000
   6.500   1.0576   0.01483   0.00730  -0.0711   0.3917   1.0000
   6.750   1.0804   0.01525   0.00766  -0.0703   0.3808   1.0000
   7.000   1.1032   0.01563   0.00805  -0.0695   0.3699   1.0000
   7.250   1.1262   0.01602   0.00848  -0.0687   0.3598   1.0000
   7.500   1.1486   0.01648   0.00890  -0.0679   0.3499   1.0000
   7.750   1.1705   0.01686   0.00933  -0.0669   0.3391   1.0000
   8.000   1.1918   0.01726   0.00978  -0.0659   0.3278   1.0000
   8.250   1.2120   0.01769   0.01023  -0.0647   0.3155   1.0000
   8.500   1.2310   0.01811   0.01066  -0.0634   0.3022   1.0000
   8.750   1.2490   0.01852   0.01110  -0.0618   0.2876   1.0000
   9.000   1.2663   0.01891   0.01155  -0.0602   0.2720   1.0000
   9.250   1.2826   0.01932   0.01203  -0.0584   0.2533   1.0000
   9.500   1.2960   0.01984   0.01253  -0.0563   0.2317   1.0000
   9.750   1.3072   0.02048   0.01314  -0.0538   0.2074   1.0000
  10.000   1.3158   0.02138   0.01392  -0.0512   0.1848   1.0000
  10.250   1.3246   0.02239   0.01486  -0.0487   0.1663   1.0000
  10.500   1.3333   0.02347   0.01590  -0.0464   0.1524   1.0000
  10.750   1.3420   0.02458   0.01699  -0.0442   0.1408   1.0000
  11.000   1.3494   0.02581   0.01818  -0.0421   0.1314   1.0000
  11.250   1.3609   0.02682   0.01927  -0.0404   0.1222   1.0000
  11.500   1.3704   0.02799   0.02047  -0.0387   0.1126   1.0000
  11.750   1.3775   0.02935   0.02184  -0.0370   0.1040   1.0000
  12.000   1.3859   0.03069   0.02323  -0.0355   0.0937   1.0000
  12.250   1.3920   0.03227   0.02483  -0.0339   0.0811   1.0000
	]

clark_cls = POLAR_FILE[:,2]
clark_cds = POLAR_FILE[:,3]

function calc_c_d(c_l, coeff)
	coeff[1] + coeff[2] * c_l + coeff[3] * c_l^2 +  coeff[4] * c_l ^3
end

function cost_fun(coeffs)
	sum((clark_cds .- map(c_l -> calc_c_d(c_l, coeffs), clark_cls)).^2)
end

opt = optimize(cost_fun, [0.0, 0.0, 0.0, 0.0])
minimizer = Optim.minimizer(opt)

function clark_c_d(c_l)
	calc_c_d(c_l, minimizer)
end

function adjust_cl_for_aspect_ratio(aspect, c_l, c_d)
	c_d + c_l^2  / pi / aspect
end

c_l_max = maximum(clark_cls)
lift_max_foil = 1 ./ 2 .* rho .* vv.^2 .* c_l_max .* area
actual_c_l = c_l_max .* (min.(g * m_tot, lift_max_foil) ./ lift_max_foil)
actual_c_d = adjust_cl_for_aspect_ratio.(aspect, actual_c_l, clark_c_d.(actual_c_l)) 
actual_lift = 1 ./ 2 .* rho .* vv.^2 .* actual_c_l .* area
actual_foil_drag = 1 ./ 2 .* rho .* vv.^2 .* actual_c_d .* area
power_foil = actual_foil_drag .* vv / 1000.0




# displacement hull
# bruker regnearket på https://www.boatdesign.net/threads/calculating-power-for-displacement-boats.38487/
slr = vv ./ MS_KNOTS ./ sqrt(l / M_FEET)
# Gerr Method A
shp = max.(0.0, (m_tot .- actual_lift ./ 9.81)) ./ KG_LBS ./ (10.665 ./ slr).^3
power_displacement = shp .* KW_HP


# https://en.wikipedia.org/wiki/Hull_speed
v_hull = 1.34 * sqrt(l / M_FEET) * MS_KNOTS

p0 = plot([v_hull, v_hull], [0.0, 30.0], labels = "")
p0 = plot!(p0, vv, [power_foil, power_displacement, power_foil + power_displacement], labels=["P Foil" "P Displacement" "Sum P"], lw = 2.0, ylims = (0, 30), xlims = (0,maximum(vv)), ylabel = "kW", xlabel = "m/s")
p1 = plot([v_hull, v_hull], [0.0, m_tot * g / 10], labels = "")
p1 = plot!(p1, vv, actual_lift / 10, lw = 2.0,  labels = "Lift of foil [daN]")
p2 = plot(vv, [actual_c_d actual_c_l], labels = ["C_L" "C_D"], lw = 2.0)
p3 = plot(vv, (power_foil + power_displacement) * 100 ./ vv, lw = 2.0,  labels = "daN propulsion")
plot(p0, p1, p2, p3, layout = (2,2), size = (1200, 1000))