Merge pull request #68 from symbiotic-engineering/float-stress

Float stress

mdocean/inputs/parameters.m

@@ -83,8 +83,11 @@
     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)");
+    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","h_{stiff}",{16 * in2m},"structures",true,"float stiffener height (m)");
+    table("w_stiff","w_{stiff}",{1  * in2m},"structures",true,"float stiffener width (m)");
+    table("num_sections","N_{sect}",{12},"structures",false,"number of float sections (-)");
     ...
     ...% Economics
     table("m_scale","m_{scale}",{1.1},"economics",false,...

mdocean/inputs/var_bounds.m

@@ -143,8 +143,9 @@
 b.X_start_struct = cell2struct(num2cell(b.X_starts),b.var_names(1:end-1)',1);
 
 b.constraint_names = {'float_too_heavy','float_too_light','spar_too_heavy','spar_too_light',...
-                      'stability','FOS_float_max','FOS_float_fatigue','FOS_col_max',...
-                      'FOS_col_fatigue','FOS_plate_max','FOS_plate_fatigue',...
+                      'stability','FOS_float_max','FOS_float_fatigue',...
+                      'float_top_ratio','float_radial_ratio','float_circ_ratio',... % placeholders for now
+                      'FOS_col_max','FOS_col_fatigue','FOS_plate_max','FOS_plate_fatigue',...
                       'FOS_col_local_max','FOS_col_local_fatigue',...
                       'pos_power','LCOE_max','irrelevant_max_force',...
                       'spar_height_up','spar_height_down','linear_theory'};

mdocean/optimization/gradient_optim.m

@@ -89,10 +89,11 @@
     T_s = D_s * p.T_s_over_D_s;
     h_f = T_f_2 / p.T_f_2_over_h_f;
     D_d = p.D_d_over_D_s * D_s;
+    D_f_in = D_s * p.D_f_in_over_D_s;
 
     MEEM = pi*p.harmonics / (p.besseli_argmax*2);
     prob.Constraints.linear_spar_natural_freq = D_d >= p.D_d_min;
-    prob.Constraints.linear_float_spar_diam = D_s <= D_f - .01;
+    prob.Constraints.linear_float_spar_diam = D_f_in <= D_f - .01;
     prob.Constraints.linear_float_spar_draft = T_f_2 <= T_s - .01;
     prob.Constraints.linear_float_spar_tops = h_s - T_s >= h_f - T_f_2 + .01;
     prob.Constraints.linear_float_seafloor = p.h - T_f_2 >= MEEM * D_f; % M

mdocean/simulation/modules/structures/float_plate_stress.m

@@ -0,0 +1,81 @@
+
+
+
+function [sigma_vm_bot,sigma_vm_top] = float_plate_stress(D_f, D_f_in, F_heave, num_sections, t_bot, t_top, h_stiff, w_stiff, D_f_tu, nu)
+
+    A = pi/4 * (D_f^2 - D_f_in^2);
+    P_hydrodynamic = F_heave / A;
+    W = F_heave / num_sections;
+
+    b_out = pi * D_f    / num_sections;
+    b_in  = pi * D_f_in / num_sections;
+    b_avg = sqrt(b_out * b_in);
+
+    h = (D_f-D_f_in)/2;
+    if h >= b_avg
+        m = (b_out-b_in)/(2*h);
+        shorter_length = b_in*(m+sqrt(1+m^2));
+        longer_length = h;
+    else
+        shorter_length = h;
+        longer_length  = b_out;
+    end
+
+    length_ratio = longer_length / shorter_length;
+    r0 = 1.02;%D_f_tu / 2;
+    width_plate = b_out; % fixme, should be combo of b_out and b_in? but b_in gives imag number
+
+    sigma_vm_bot = bottom_plate_stress(length_ratio, shorter_length, P_hydrodynamic, t_bot, h_stiff, w_stiff, width_plate);
+    sigma_vm_top = top_plate_stress(length_ratio, shorter_length, W, t_top, h_stiff, w_stiff, r0, width_plate, nu);
+
+end
+
+function sigma_vm_top = top_plate_stress(length_ratio, shorter_length, W, t_top, h_stiff, w_stiff, r0, width_plate, nu)
+
+    length_ratio_data = [1 : 0.2 : 2 1000];
+    beta_1_data = [-0.238 -0.078 0.011 0.053 0.068 0.067 0.067];
+    beta_2_data = [0.7542 0.8940 0.9624 0.9906 1.0000 1.004 1.008];
+    beta_1_fcn  = @(length_ratio) interp1(length_ratio_data, beta_1_data, length_ratio);
+    beta_2_fcn  = @(length_ratio) interp1(length_ratio_data, beta_2_data,  length_ratio);
+    beta_1 = beta_1_fcn(length_ratio);
+    beta_2 = beta_2_fcn(length_ratio);
+
+    % Roark's table 11.4 case 8b page 514
+    % fixme - this is potentially not applicable since it assumes force in
+    % middle, but my force is on the edge and therefore creates ~no bending moment
+    % perhaps check ref 26 and 21 to see if they have a more general (off center) equation 
+    sigma_edge_no_stiff = 3*W/(2*pi*t_top^2) * ( (1+nu)*log(2*shorter_length/(pi*r0)) + beta_1);
+    sigma_cent_no_stiff = -beta_2 * W / t_top^2;
+
+    M_edge = sigma_edge_no_stiff * t_top^2 / 6;
+    M_cent = sigma_cent_no_stiff * t_top^2 / 6;
+
+    M_max = max(abs([M_edge,M_cent]));
+    sigma_max = get_stiffener_equivalent_stress(t_top, h_stiff, width_plate, w_stiff, M_max);
+    sigma_vm_top = sigma_max;
+
+end
+
+function sigma_vm = bottom_plate_stress(length_ratio, shorter_length, P_hydrodynamic, t_bot, h_stiff, w_stiff, width_plate)
+
+    % data from table: Timoshenko table 35 p219
+    length_ratio_data = [1 : 0.1 : 2 1000];
+    alpha_shorter_data = -0.0001*[513 581 639 687 726 757 780 799 812 822 829 833];
+    alpha_longer_data  = -0.0001*[513 538 554 563 568 570 571 571 571 571 571 571];
+    alpha_shorter_fcn = @(length_ratio) interp1(length_ratio_data, alpha_shorter_data, length_ratio);
+    alpha_longer_fcn  = @(length_ratio) interp1(length_ratio_data, alpha_longer_data,  length_ratio);
+
+    M_shorter = alpha_shorter_fcn(length_ratio) * P_hydrodynamic * shorter_length^2;
+    M_longer  = alpha_longer_fcn(length_ratio)  * P_hydrodynamic * shorter_length^2;
+
+    sigma_shorter = get_stiffener_equivalent_stress(t_bot, h_stiff, width_plate, w_stiff, M_shorter);
+    sigma_longer  = get_stiffener_equivalent_stress(t_bot, h_stiff, width_plate, w_stiff, M_longer);
+
+    sigma_zz = P_hydrodynamic;
+
+    s_11 = sigma_shorter; s_22 = 0; s_33 = sigma_zz;
+    sigma_vm = sqrt(1/2 * ( (s_11 - s_22).^2 + (s_22 - s_33).^2 + (s_33 - s_11).^2 ));
+
+end
+
+

mdocean/simulation/modules/structures/get_stiffener_equivalent_stress.m

@@ -0,0 +1,23 @@
+function sigma = get_stiffener_equivalent_stress(t_plate, h_stiff, width_plate, width_stiff, moment_per_length)
+% https://ocw.mit.edu/courses/2-080j-structural-mechanics-fall-2013/f8fd2ad49d100766335b4e129a8a4791_MIT2_080JF13_Lecture7.pdf
+
+    h = t_plate;
+    H = h_stiff;
+    a = width_plate/2;
+    b = width_stiff;
+    M = moment_per_length;
+
+    eta_over_h = 1/2 * (1 - b/a * (H/h)^2) / (1 + b/a * (H/h)); % eq 7.67 MIT 2.080
+    h_eq_over_h_3 = 4 * ( (1-3*eta_over_h+3*eta_over_h^2) + b/(2*a) * ( (H/h)^2 + 3*(H/h)^2*eta_over_h + 3*H/h*eta_over_h^2 )  ); % eq 7.69
+    if h_eq_over_h_3 < 0
+        warning('stiffener equivalence broke')
+        h_eq_over_h_3 = 1;
+    end
+    h_eq_over_h = h_eq_over_h_3^(1/3);
+    h_eq = h_eq_over_h * h;
+
+    y_max = h-eta_over_h*h;
+
+    sigma = 12 * M * y_max / h_eq^3;
+
+end

mdocean/simulation/modules/structures.m → mdocean/simulation/modules/structures/structures.m

@@ -1,21 +1,32 @@
 function [FOS1Y,FOS2Y,FOS3Y,FOS_buckling] = structures(...
-          	F_heave_storm, F_surge_storm, F_heave_op, F_surge_op, ...
-            M, h_s, T_s, rho_w, g, sigma_y, sigma_e, A_c, A_lat_sub, D_s, t_s_r, I, E, nu)
+          	F_heave_storm, F_surge_storm, F_heave_op, F_surge_op, ... % forces
+            h_s, T_s, D_s, D_f, D_f_in, num_sections, D_f_tu, ... % bulk dimensions
+            t_s_r, I, A_c, A_lat_sub, t_bot, t_top, h_stiff, w_stiff, ... % structural dimensions
+            M, rho_w, g, sigma_y, sigma_e, E, nu) % constants
 
     F_heave_peak = max(F_heave_storm,F_heave_op);
     F_surge_peak = max(F_surge_storm,F_surge_op);
 
     % peak
-    [FOS1Y, FOS2Y, FOS3Y, FOS_buckling] = structures_one_case(F_heave_peak, F_surge_peak, ...
-            M, h_s, T_s, rho_w, g, sigma_y(M), A_c, A_lat_sub, D_s, t_s_r, I, E, nu);
-
+    [FOS1Y, FOS2Y, FOS3Y, FOS_buckling] = structures_one_case(...
+            F_heave_peak, F_surge_peak, sigma_y(M), ...
+            h_s, T_s, D_s, D_f, D_f_in, num_sections, D_f_tu, ...
+            t_s_r, I, A_c, A_lat_sub, t_bot, t_top, h_stiff, w_stiff, ...
+            M, rho_w, g, E, nu);
+
     % endurance limit (fatigue)
-    [FOS1Y(2), FOS2Y(2), FOS3Y(2), FOS_buckling(2)] = structures_one_case(F_heave_op, F_surge_op, ...
-            M, h_s, T_s, rho_w, g, sigma_e(M), A_c, A_lat_sub, D_s, t_s_r, I, E, nu);
+    [FOS1Y(2), FOS2Y(2), FOS3Y(2), FOS_buckling(2)] = structures_one_case(...
+            F_heave_op, F_surge_op, sigma_e(M), ...
+            h_s, T_s, D_s, D_f, D_f_in, num_sections, D_f_tu, ...
+            t_s_r, I, A_c, A_lat_sub, t_bot, t_top, h_stiff, w_stiff, ...
+            M, rho_w, g, E, nu);
 end
 
-function [FOS1Y, FOS2Y, FOS3Y, FOS_spar_local] = structures_one_case(F_heave, F_surge, ...
-            M, h_s, T_s, rho_w, g, sigma_max, A_c, A_lat_sub, D_s, t_s_r, I, E, nu)
+function [FOS1Y, FOS2Y, FOS3Y, FOS_spar_local] = structures_one_case(...
+            F_heave, F_surge, sigma_max, ...
+            h_s, T_s, D_s, D_f, D_f_in, num_sections, D_f_tu, ...
+            t_s_r, I, A_c, A_lat_sub, t_bot, t_top, h_stiff, w_stiff, ...
+            M, rho_w, g, E, nu)
 
     %% Stress calculations
     depth = [0 T_s T_s]; % max depth
@@ -41,15 +52,18 @@
     % assume ductile material for now - need to use mohr's circle for concrete
     sigma_vm = von_mises(sigma_rr, sigma_tt, sigma_zz, sigma_rt, sigma_tz, sigma_zr);
 
+    %% Float plate stress
+    [sigma_float_bot,sigma_float_top] = float_plate_stress(D_f, D_f_in, F_heave, num_sections, t_bot, t_top, h_stiff, w_stiff, D_f_tu, nu);
+    % ignoring top and side plates for now
+
     %% Buckling calculation
     [FOS_spar,FOS_spar_local] = spar_combined_buckling(F_heave, E(M), I(2), h_s, D_s, A_c(2), t_s_r, ...
                                         P_hydrostatic(2), sigma_max, nu(M));
 
     %% Factor of Safety (FOS) Calculations
-    FOS_yield = sigma_max ./ sigma_vm;
-    FOS1Y = FOS_yield(1);
+    FOS1Y = sigma_max / sigma_float_bot;
     FOS2Y = FOS_spar;
-    FOS3Y = FOS_yield(3);
+    FOS3Y = sigma_vm(3);
 
 end 
 

mdocean/simulation/simulation.m

@@ -35,7 +35,8 @@
 in.D_f_b = p.D_f_b_over_D_f * in.D_f;
 
 % Space between float and spar
-in.D_f_in = p.D_f_in_over_D_s * in.D_s;
+D_f_in = p.D_f_in_over_D_s * in.D_s;
+in.D_f_in = min(D_f_in, in.D_f - 1e-6); % ensure positive area to avoid sim breaking for infeasible inputs
 
 %% Run modules
 [V_d, m_m, m_f_tot, m_s_tot, ...
@@ -56,14 +57,16 @@
                     X_constraints] = dynamics(in, m_f_tot, m_s_tot, V_d, T);
 
 [FOS_float,FOS_spar,FOS_damping_plate,...
-    FOS_spar_local] = structures(F_heave_storm, F_surge_storm, F_heave_op, F_surge_op, ...
-                                 in.M, in.h_s, in.T_s, in.rho_w, in.g, in.sigma_y, in.sigma_e, ...
-                                 A_c, A_lat_sub, in.D_s, in.t_s_r, I, in.E, in.nu);
+    FOS_spar_local] = structures(...
+          	F_heave_storm, F_surge_storm, F_heave_op, F_surge_op, ... % forces
+            in.h_s, in.T_s, in.D_s, in.D_f, in.D_f_in, in.num_sections, in.D_f_tu, ... % bulk dimensions
+            in.t_s_r, I, A_c, A_lat_sub, in.t_f_b, in.t_f_t, in.h_stiff, in.w_stiff, ... % structural dimensions
+            in.M, in.rho_w, in.g, in.sigma_y, in.sigma_e, in.E, in.nu);
 
 LCOE = econ(m_m, in.M, in.cost_m, in.N_WEC, P_avg_elec, in.FCR, in.cost_perN, in.F_max, in.eff_array);
 
 %% Assemble constraints g(x) >= 0
-num_g = 20+numel(p.JPD);
+num_g = 23+numel(p.JPD);
 g = zeros(1,num_g);
 g(1) = V_f_pct;                         % prevent float too heavy
 g(2) = 1 - V_f_pct;                     % prevent float too light
@@ -72,26 +75,28 @@
 g(5) = GM;                              % pitch stability of float-spar system
 g(6) = FOS_float(1) / p.FOS_min - 1;    % float survives max force
 g(7) = FOS_float(2) / p.FOS_min - 1;    % float survives fatigue
-g(8) = FOS_spar(1) / p.FOS_min - 1;           % spar survives max force
-g(9) = FOS_spar(2) / p.FOS_min - 1;           % spar survives fatigue
-g(10) = FOS_damping_plate(1) / p.FOS_min - 1; % damping plate survives max force
-g(11) = FOS_damping_plate(2) / p.FOS_min - 1; % damping plate survives fatigue
-g(12) = FOS_spar_local(1) / p.FOS_min - 1;    % spar survives max force in local buckling
-g(13) = FOS_spar_local(2) / p.FOS_min - 1;    % spar survives fatigue in local buckling
-g(14) = P_avg_elec;                     % positive power
+g(8) = in.t_f_t / in.t_f_b - 0.5 / 0.56;      % float top thickness ratio placeholder
+g(9) = in.t_f_r / in.t_f_b - 0.44 / 0.56;     % float radial thickness ratio placeholder
+g(10) = in.t_f_c / in.t_f_b - 0.44 / 0.56;    % float circumferential thickness ratio placeholder
+g(11) = FOS_spar(1) / p.FOS_min - 1;           % spar survives max force
+g(12) = FOS_spar(2) / p.FOS_min - 1;           % spar survives fatigue
+g(13) = FOS_damping_plate(1) / p.FOS_min - 1; % damping plate survives max force
+g(14) = FOS_damping_plate(2) / p.FOS_min - 1; % damping plate survives fatigue
+g(15) = FOS_spar_local(1) / p.FOS_min - 1;    % spar survives max force in local buckling
+g(16) = FOS_spar_local(2) / p.FOS_min - 1;    % spar survives fatigue in local buckling
+g(17) = P_avg_elec;                     % positive power
 %1 + min(Kp_over_Ks,[],'all');   % spar heave stability (positive effective stiffness)
-g(15) = p.LCOE_max/LCOE - 1;            % prevent more expensive than threshold
-g(16) = F_ptrain_max/in.F_max - 1;      % prevent irrelevant max force -
+g(18) = p.LCOE_max/LCOE - 1;            % prevent more expensive than threshold
+g(19) = F_ptrain_max/in.F_max - 1;      % prevent irrelevant max force -
                                         % this constraint should always be active
                                         % and is only required when p.cost_perN = 0.
-g(17) = X_constraints(1);               % prevent float rising above top of spar
-g(18) = X_constraints(2);               % prevent float going below bottom of spar
-g(19) = X_constraints(3);               % prevent float support tube (PTO attachment) from hitting spar
-g(20) = X_constraints(4);               % float amplitude obeys linear theory
-g(21:end) = X_constraints(5:end);       % prevent rising out of water/slamming
-
-criteria = all(~isinf(g)) && all(~isnan(g)) && all(isreal(g));
-%assert( criteria )
+g(20) = X_constraints(1);               % prevent float rising above top of spar
+g(21) = X_constraints(2);               % prevent float going below bottom of spar
+g(22) = X_constraints(3);               % prevent float support tube (PTO attachment) from hitting spar
+g(23) = X_constraints(4);               % float amplitude obeys linear theory
+g(24:end) = X_constraints(5:end);       % prevent rising out of water/slamming
+
+criteria = all(~isinf([g LCOE P_var])) && all(~isnan([g LCOE P_var])) && all(isreal([g LCOE P_var]));
 if ~criteria
     warning('Inf, NaN, or imaginary constraint detected')
 end