More writing updates

.gitignore

@@ -45,6 +45,8 @@ test-results/
 # figures
 *.fig
 *.png
+# except intentional figures
+!**/non_matlab_figs/*.png
 
 # mat data files
 *.mat

dev/dynamics/multibody/multibody_fudge_factor.m

@@ -0,0 +1,22 @@
+p = parameters();
+b = var_bounds();
+X = [b.X_noms; 1];
+wecsim_filename = 'wecsim_sparcd0_floatcd0_14e825b_01a25103-1490-4189-b376-76cbdf616634';
+report = true;
+override = false;
+power_matrix_compare(X, p, wecsim_filename, report, override)
+
+% add breakpoint inside power_matrix_compare and do the following in
+% command window: (see 1/27/25 slides for screenshots)
+sp = sim.power_mech_unsat;
+ap = actual.power_mech_unsat;
+ratio = ap ./ sp
+
+figure; plot(T,ratio)
+ylabel('RM3 report power / wecsim singlebody power')
+xlabel('T')
+
+y=ratio(H==5.75,:);
+figure; plot(T,y)
+
+% then do apps > curve fitter with rational function: num degree 0, den degree 2

dev/econ/design_cost_scaling.m

@@ -0,0 +1,50 @@
+%%%% PTO cost %%%%%%
+% PTO cost scale with N_WEC (from curve fit)
+% C_PTO_per_wec = C_PTO_const + C_PTO_coeff * N_WEC ^ -C_PTO_alpha
+C_PTO_const = 0.280e6;
+C_PTO_coeff = 0.344e6;
+C_PTO_alpha = 0.206;
+
+% cost that scales with force: generator + bearing + mounting
+cost_force_N1  = 25740+88411+28600;
+cost_force_N100 = 20037+66308+5720;
+
+F_max_nom = 0.926e6;
+C_force_coeff = (cost_force_N1 - cost_force_N100) / (F_max_nom * (1 - 100^(-C_PTO_alpha)))
+C_force_const = (cost_force_N1 / F_max_nom) - C_force_coeff
+
+% cost that scales with neither force nor power: riser cable and controller
+cost_neither_N1  = 88000+5644;
+cost_neither_N100 = 88000+5000;
+
+C_neither_coeff = (cost_neither_N1 - cost_neither_N100) / ((1 - 100^(-C_PTO_alpha)))
+C_neither_const = (cost_neither_N1) - C_neither_coeff
+
+P_max_nom = 286000;
+C_power_coeff = (C_PTO_coeff - C_neither_coeff - C_force_coeff*F_max_nom) / P_max_nom
+C_power_const = (C_PTO_const - C_neither_const - C_force_const*F_max_nom) / P_max_nom
+
+
+cost_force_pred = (C_force_const + C_force_coeff*[1 100].^-C_PTO_alpha) * F_max_nom;
+cost_force_act = [cost_force_N1 cost_force_N100];
+
+cost_neither_pred = C_neither_const + C_neither_coeff*[1 100].^-C_PTO_alpha;
+cost_neither_act = [cost_neither_N1 cost_neither_N100];
+
+cost_pto_pred = C_PTO_const + C_PTO_coeff*[1 100].^-C_PTO_alpha
+cost_power_pred = cost_pto_pred - cost_neither_pred - cost_neither_pred
+
+%%%%%% non design cost %%%%%%%%
+clc
+N_WEC = [1 10 50 100];
+development     = [4553389 8773812 11003159 10820060] ./ N_WEC
+infrastructure  = [990000 4860000 7566000 17310000] ./ N_WEC
+mooring         = [524775 4722975 23614875 47229750] ./ N_WEC
+profitmargin    = [356252 2561152 11323295 21921723] ./ N_WEC
+installation    = [5908552 9081973 21531225 37859591] ./ N_WEC
+decommissioning = 0;%installation; % decomissioning cost not used in LCOE
+contingency     = [1589545 5561144 18827150 35435836] ./ N_WEC
+
+non_design_capex = development + infrastructure + mooring + profitmargin + installation + decommissioning + contingency;
+
+

dev/force_sat/quadratic.m

@@ -68,7 +68,8 @@
 matrix = simplify(matrix,'IgnoreAnalyticConstraints',true)
 
 latex(matrix)
-bases_rk = bases .* [r_k^2, r_k^2*r_b, r_k^2*r_b^2, r_k*r_b^2, r_b^2];
+rkb = [r_k^2, r_k^2*r_b, r_k^2*r_b^2, r_k*r_b^2, r_b^2];
+bases_rk = bases .* rkb;
 pretty(bases_rk')
 latex(bases_rk)
 
@@ -77,6 +78,51 @@
 soln_f = solve(eqn_m_1,f);
 assert(soln_f == 1)
 
+%% replace rk with other variables
+% notebook p97 2-1-25
+syms omega_u_star zeta_u star real positive
+
+rk_new = star^2/(star^2-1);
+rb_new = 1/(1-zeta_u/(zeta*star));
+w_rat_new = omega_u_star / star;
+
+rkb_new = subs(rkb,[r_k,r_b],[rk_new, rb_new]);
+pretty(rkb_new.')
+bases_new = subs(bases,w_ratio,w_rat_new);
+bases_rk_new = bases_new .* rkb_new;
+scale = (star^2-1)^2 * (star*zeta)^-2 * (zeta_u - star*zeta)^2;
+pretty(simplify(bases_rk_new.' * scale,'IgnoreAnalyticConstraints',true))
+
+%% check that parameterization used in paper equals "new" just above
+% notebook p99 2-1-25
+matrix_check = [-4  4/f^2 1/f^2 -1  0  0;
+                 8  0     0      2  8  2;
+                -4 -4    -1     -1 -8 -2];
+bases_unsquared_check = [omega_u_star * (zeta_u - star*zeta),...
+                         omega_u_star*star*zeta,...
+                         omega_u_star^2 - star^2,...
+                         star^2-1];
+bases_squared_check = [bases_unsquared_check(1)^2;...
+                       bases_unsquared_check(2)^2;...
+                       bases_unsquared_check(3)^2;...
+                       bases_unsquared_check(4)^2;...
+                       bases_unsquared_check(1)*bases_unsquared_check(2);...
+                       bases_unsquared_check(3)*bases_unsquared_check(4)];
+
+bases_should_match = [bases_unsquared_check(1)^2;... % a^2
+                     -bases_unsquared_check(1) * bases_unsquared_check(2);... % -ab
+                      bases_unsquared_check(3)^2 + 4*bases_unsquared_check(2)^2;... % c^2+4b^2
+                      bases_unsquared_check(3) * bases_unsquared_check(4);... % cd
+                      bases_unsquared_check(4)^2]; % d^2
+
+should_be_1 = simplify( bases_should_match ./ (bases_rk_new.' * scale) )
+assert(all(should_be_1 == 1))
+
+%% and check that it equals the original (r_k, r_b, quad)
+quad_new = subs(quad,[r_k,r_b,w_ratio],[rk_new,rb_new,w_rat_new]);
+should_be_1_also = simplify( quad_new.' ./ (matrix_check*bases_squared_check/scale) )
+assert(all(should_be_1_also == 1))
+
 %% optimal powertrain case
 bases_rk = subs(bases_rk,{w_ratio,r_b},{1,2}).';
 

mdocean/inputs/parameters.m

@@ -42,7 +42,7 @@
     table("rho_w","\rho_w",{1000},"site",false,"water density (kg/m3)");
     ...%table("mu","\mu",1e-3,"site",false,"dynamic viscosity of water (Pa s)");
     table("g","g",{g},"site",false,"acceleration of gravity (m/s2)");
-    table("h","h",{100},"site",true,"water depth (m)");
+    table("h","h",{45},"site",true,"water depth (m)");
     table("JPD","JPD",{jpd(2:end,2:end)},"site",false,...
         "joint probability distribution of wave (%)");
     table("Hs","H_s",{jpd(2:end,1)},"site",true,"significant wave height (m)");
@@ -61,36 +61,41 @@
     table("rho_m","\rho_m",{[7850 2400 7900]},"structures",true,"material density (kg/m3)");
     table("E","E",{[200e9, 5000*sqrt(4.5*ksi2pa), 200e9]},"structures",true,...
         "young's modulus (Pa)");
-    table("cost_m","cost_m",{[4.28, 125/yd2m^3/2400, 1.84/lb2kg]},"economics",true,"material cost ($/kg)");
+    table("cost_perkg_mult","cost_m",{[4.28, 125/yd2m^3/2400, 1.84/lb2kg]/4.28},"economics",true,"material cost ($/kg)");
         ...% RM3 CBS sheet 1.4 average of cells F21, F34, F46, F56
         ...% https://www.concretenetwork.com/concrete-prices.html
         ...% https://agmetalminer.com/metal-prices/
     table("nu","\nu",{[0.36 0 0.29]},"structures",true,"Poisson's ratio (-)");
+    table("FOS_min","FOS_{min}",{1.5},"structures",true,"minimum FOS (-)");
+    ...
+    ...% Thicknesses and structures: float
+    table("t_f_t_over_t_f_b","t_{f,t}/t_{f,b}",{0.50/0.56},"structures",true,...
+        "float top to bottom thickness ratio (-)");
+    table("t_f_r_over_t_f_b","t_{f,r}/t_{f,b}",{0.44/0.56},"structures",true,...
+        "float radial to bottom thickness ratio (-)");
+    table("t_f_c_over_t_f_b","t_{f,c}/t_{f,b}",{0.44/0.56},"structures",true,...
+        "float circumferential to bottom thickness ratio (-)");
+    table("D_f_tu","D_{f,tu}",{20 * in2m},"structures",true,...
+        "float support tube diameter (m)"); % 24 in p156 report, 20 in cad
+    table("t_f_tu","t_{f,tu}",{.5 * in2m},"structures",true,...
+        "float support tube thickness (m)");
+    table("w_over_h_stiff_f","w_{stiff,f}/h_{stiff,f}",{1/16},"structures",...
+        true,"float stiffener width to height ratio (-)");
+    table("num_sections_f","N_{sect}",{12},"structures",false,"number of float sections (-)");
     ...
-    ...% Thicknesses and structures
-...%     table("t_f_t","t_{f,t}",{0.50 * in2m},"structures",true,"float top thickness (m)");
-...%     table("t_f_r","t_{f,r}",{0.44 * in2m},"structures",true,"float radial wall thickness (m)");
-...%     table("t_f_c","t_{f,c}",{0.44 * in2m},"structures",true,...
-...%         "float circumferential gusset thickness (m)");
-...%     table("t_f_b","t_{f,b}",{0.56 * in2m},"structures",true,"float bottom thickness (m)");
-...%     table("t_s_r","t_{s,r}",{1.00 * in2m},"structures",true,"vertical column thickness (m)");
-    table("t_d_max","t_{d,max}",{1.00 * in2m},"structures",true,...
-        "max thickness of damping plate before making it hollow (m)");
-...%     table("t_d_tu","t_{d,tu}",{1.00 * in2m},"structures",true,...
-...%         "damping plate support tube radial wall thickness (m)");
+    ... Thicknesses and structures: damping plate
+    table("t_d_tu","t_{d,tu}",{1.00 * in2m},"structures",true,...
+        "damping plate support tube radial wall thickness (m)");
     table("D_d_tu","D_{d,tu}",{48.00 * in2m},"structures",true,...
         "damping plate support tube diameter (m)");
     table("theta_d_tu","\theta_{d,tu}",{atan(17.5/15)},"structures",true,...
         "angle from horizontal of damping plate support tubes (rad)");
-    table("FOS_min","FOS_{min}",{1.5},"structures",true,"minimum FOS (-)");
-    table("D_f_tu","D_{f,tu}",{20 * in2m},"structures",true,"float support tube diameter (m)"); % 24 in p156 report, 20 in cad
-    table("t_f_tu","t_{f,tu}",{.5 * in2m},"structures",true,"float support tube thickness (m)");
-    table("h_stiff_f","h_{stiff,f}",{16 * in2m},"structures",true,"float stiffener height (m)");
-    table("w_stiff_f","w_{stiff,f}",{1  * in2m},"structures",true,"float stiffener width (m)");
-    table("h_stiff_d","h_{stiff,d}",{[12.5 .5 22  1]*in2m},"structures",true,"damping plate stiffener height (m)");
-    table("w_stiff_d","w_{stiff,d}",{[.5 10 1 12]  * in2m},"structures",true,"damping plate stiffener width (m)");
-    table("num_sections","N_{sect}",{12},"structures",false,"number of float sections (-)");
-    table("FOS_mult_d","FOS_{mult,d}",{7.5},"structures",true,"damping plate factor of safety multiplier (-)");
+    table("h_over_h1_stiff_d","h_{stiff,d}/h_{1,stiff,d}",{[12.5 0.5 22  1]/22},...
+        "structures",true,"damping plate stiffener height (m)");
+    table("w_over_h1_stiff_d","w_{stiff,d}/h_{1,stiff,d}",{[.5 10 1 12]/22},...
+        "structures",true,"damping plate stiffener width (m)");
+    table("FOS_mult_d","FOS_{mult,d}",{7.5},"structures",true,...
+        "damping plate factor of safety multiplier (-)");
     table("num_terms_plate","N_{plate}",{100},"structures",false,...
         "number of terms for damping plate concentrated load (-)");
     table("radial_mesh_plate","N_{r,plate}",{20},"structures",false,...
@@ -99,22 +104,24 @@
         "number of damping plate stiffeners (-)");
     ...
     ...% Economics
-    table("m_scale","m_{scale}",{1.1},"economics",false,...
-        "factor to account for mass of neglected stiffeners (-)");
+    table("m_scale","m_{scale}",{1.10},"economics",false,...
+        "factor to account for mass of neglected stiffeners and device access (-)");
+    table("power_scale","P_{scale}",{85.9/117.7},"economics",false,...
+        "factor to scale power for validation tuning (-)")
     table("FCR","FCR",{0.113},"economics",true,...
         "fixed charge rate (-), see RM3 report p63");
     table("N_WEC","N_{WEC}",{100},"economics",true,"number of WECs in array (-)");
-    table("LCOE_max","LCOE_{max}",{.5},"economics",true,"maximum LCOE ($/kWh)");
-    table("avg_power_min","P_{avg,elec,min}",{100},"economics",true,...
-        "minimum average electrical power (W)"); % set to a negative number (not zero) to disable constraint
+    table("LCOE_max","LCOE_{max}",{.75},"economics",true,"maximum LCOE ($/kWh)");
+...%     table("avg_power_min","P_{avg,elec,min}",{100},"economics",true,...
+...%         "minimum average electrical power (W)"); % set to a negative number (not zero) to disable constraint
     table("eff_array","\eta_{array}",{0.95*0.98},"economics",true,...
         "array availability and transmission efficiency (-)");
-    table("cost_perN", "$/N",{0.1656},"economics",true,...
-        "cost per Newton ($/N), from https://doi.org/10.1016/j.ifacol.2022.10.531");
-    table("cost_perW", "$/W",{1},"economics",true,...
-        "cost per Watt ($/W)");
+    table("cost_perN_mult", "$/N",{1},"economics",true,...
+        "cost per Newton multiplier (-), 0.1656 from https://doi.org/10.1016/j.ifacol.2022.10.531");
+    table("cost_perW_mult", "$/W",{1},"economics",true,...
+        "cost per Watt multiplier (-)");
     ...
-    ...% Geometric ratios
+    ...% Geometric ratios of bulk dimensions
     table("D_d_min","D_{d,min}",{30},"geometry",true,...
         "minimum damping plate diameter");
     table("D_d_over_D_s","D_d/D_s",{30/6},"geometry",true,...
@@ -133,14 +140,13 @@
         true,"ratio of float inner diameter to spar diameter (-)")
     ...
     ...% Dynamics: device parameters
-    table("C_d_float","C_{d,float}",{0},"dynamics",true,"coefficient of drag for float");
-    table("C_d_spar","C_{d,spar}",{0},"dynamics",true,"spar coefficient of drag");
-    table("power_max","power_{max}",{Inf},"dynamics",true,"maximum power (W)");
+    table("C_d_float","C_{d,float}",{1},"dynamics",true,"coefficient of drag for float");
+    table("C_d_spar","C_{d,spar}",{1},"dynamics",true,"spar coefficient of drag");
     table("eff_pto","\eta_{pto}",{0.80},"dynamics",true,"PTO efficiency (-)");
     ...
     ...% Dynamics: simulation type
     table("control_type","control type",{'damping'},"dynamics",false,...
-        "reactive or constant impedance or damping");
+        "reactive or damping");
     table("use_MEEM","use_MEEM",{true},"dynamics",false,...
         "whether to use MEEM for hydro coeffs (boolean)");
     table("use_multibody","use_multibody",{false},"dynamics",false,...