Merge pull request #67 from symbiotic-engineering/spar-stress

spar buckling runs when set to compact

mdocean/inputs/parameters.m

@@ -52,15 +52,17 @@
     ...
     ...% Materials: [  Structural Steel ASTM-A36 (Ductile)
     ...%               Concrete (Brittle)
-    ...%               Stainless Steel 304 (Ductile) ]
+    ...%               Stainless Steel 316 (Ductile) ASTM-A240 ]
     table("sigma_y","\sigma_y",{[36,4.5,30]* ksi2pa},"structures",true,"yield strength (Pa)");
+    table("sigma_e","\sigma_e",{[58*.45, 0, 75*.45]*ksi2pa},"structures",true,"endurance limit (Pa)")
     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)");
         ...% 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 (-)");
     ...
     ...% Thicknesses and structures
     table("t_f_t","t_{f,t}",{0.50 * in2m},"structures",true,"float top thickness (m)");

mdocean/inputs/validation/validation_inputs.m

@@ -17,7 +17,9 @@
                 'power_avg',    85.9e3, ... % S14 in CBS Performance & Economics
                 'power_max',    286e3, ...  % S15 in CBS Performance & Economics
                 'force_heave',  8500e3,...  % p156 RM3 report
-                'FOS_b',        3,...       % p158 RM3 report
+                'FOS_spar',     11.1,...    % p158 RM3 report, but modified 
+                ...                         % assuming that RM3 report used r 
+                ...                         % instead of D when calculating I
 				'c_v',			nominal_c_v()); % calculated from data in CBS Performance & Economics
 
 RM3_wecsim = struct(...

mdocean/inputs/var_bounds.m

@@ -83,15 +83,17 @@
 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_yield','FOS_col_yield','FOS_plate','FOS_col_buckling',...
+                      'stability','FOS_float_max','FOS_float_fatigue','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'};
 i1 = length(b.constraint_names);
 for i = (i1+1):(i1+14*15)
     b.constraint_names{i} = strcat('prevent_slamming',num2str(i-i1));
 end
 
-b.lin_constraint_names = {'spar_natural_freq','float_spar_diam','float_spar_draft',
+b.lin_constraint_names = {'spar_natural_freq','float_spar_diam','float_spar_draft',...
                           'float_spar_tops','float_seafloor','spar_seafloor'};
 
 [~,idxs_sort] = sort(b.var_names(1:end-1)); % alphabetical design variable indices

mdocean/simulation/modules/dynamics/dynamics.m

@@ -1,9 +1,12 @@
-function [F_heave_max, F_surge_max, F_ptrain_max, ...
+function [F_heave_storm, F_surge_storm, F_heave_op, F_surge_op, F_ptrain_max, ...
     P_var, P_avg_elec, P_matrix_elec, X_constraints, B_p, X_u, X_f, X_s, P_matrix_mech] = dynamics(in,m_float,m_spar,V_d,draft)
 
     % use probabilistic sea states for power and PTO force and max amplitude
     [T,Hs] = meshgrid(in.T,in.Hs);
-    [P_matrix_mech,X_constraints,B_p,X_u,X_f,X_s,~,~,F_ptrain_max] = get_power_force(in,T,Hs,m_float,m_spar,V_d,draft);
+    [P_matrix_mech,X_constraints,B_p,...
+        X_u,X_f,X_s,F_heave_op,...
+        F_surge_op,F_ptrain_max] = get_power_force(in,T,Hs,m_float,m_spar,...
+                                                        V_d,draft,in.F_max);
 
     % account for powertrain electrical losses
     P_matrix_elec = P_matrix_mech * in.eff_pto;
@@ -18,8 +21,8 @@
     assert(isreal(P_avg_elec))
 
     % use max sea states for structural forces
-    [~,~,~,~,~,~,F_heave_max,F_surge_max,~] = get_power_force(in, ...
-                                in.T_struct, in.Hs_struct, m_float, m_spar, V_d, draft);
+    [~,~,~,~,~,~,F_heave_storm,F_surge_storm,~] = get_power_force(in, ...
+                                in.T_struct, in.Hs_struct, m_float, m_spar, V_d, draft, in.F_max);
 
     % coefficient of variance (normalized standard deviation) of power
     P_var = std(P_matrix_elec(:), in.JPD(:)) / P_avg_elec;
@@ -29,7 +32,7 @@
 end
 
 function [P_matrix, X_constraints, B_p, mag_X_u, mag_X_f, mag_X_s,...
-          F_heave_f, F_surge, F_ptrain_max] = get_power_force(in,T,Hs, m_float, m_spar, V_d, draft)
+          F_heave_f, F_surge, F_ptrain_max] = get_power_force(in,T,Hs, m_float, m_spar, V_d, draft, F_max)
 
     % get dynamic coefficients for float and spar
     % fixme: eventually should use in.D_f_in to allow a radial gap between float and spar
@@ -52,7 +55,7 @@
      mag_X_f,phase_X_f,...
      mag_X_s,phase_X_s,...
      B_p,K_p] = get_response_drag(w,m_f,m_s,m_c,B_h_f,B_h_s,B_c,K_h_f,K_h_s,...
-                                            F_f_mag,F_f_phase,F_s_mag,F_s_phase,in.F_max,...
+                                            F_f_mag,F_f_phase,F_s_mag,F_s_phase,F_max,...
                                             drag_const_f,drag_const_s,mag_v0_f,mag_v0_s, ...
                                             X_max,in.control_type,in.use_multibody,...
                                             in.X_tol,in.phase_X_tol,in.max_drag_iters);
@@ -87,13 +90,14 @@
     if nargout > 3
         % powertrain force
         F_ptrain_max = max(mag_U,[],'all');
-        F_ptrain_max = min(F_ptrain_max, in.F_max);
+        F_ptrain_max = min(F_ptrain_max, F_max);
 
-        % heave force: includes powertrain force and D'Alembert force
-        %F_heave_fund = sqrt( (mult .* B_p .* w).^2 + (mult .* K_p - m_float .* w.^2).^2 ) .* X;
-
-        F_heave_f = combine_ptrain_dalembert_forces(m_float, w, mag_X_f, phase_X_f, mag_U, phase_U, in.F_max);
-        F_heave_s = combine_ptrain_dalembert_forces(m_spar,  w, mag_X_s, phase_X_s, mag_U, phase_U, in.F_max);
+        % heave force: includes powertrain force and D'Alembert force        
+        F_heave_f = combine_ptrain_dalembert_forces(m_float, w, mag_X_f, phase_X_f, mag_U, phase_U, F_max);
+        F_heave_s = combine_ptrain_dalembert_forces(m_spar,  w, mag_X_s, phase_X_s, mag_U, phase_U, F_max);
+
+        F_heave_f = max(F_heave_f,[],'all');
+        F_heave_s = max(F_heave_s,[],'all');
 
         % surge force
         F_surge = max(Hs,[],'all') * in.rho_w * in.g * V_d .* (1 - exp(-max(k_wvn,[],'all')*draft));
@@ -266,7 +270,7 @@
     end
 
     % get force-saturated response
-    r = min(F_max ./ mag_U_unsat, 1);%fcn2optimexpr(@min, in.F_max ./ F_ptrain_unsat, 1);
+    r = min(F_max ./ mag_U_unsat, 1);%fcn2optimexpr(@min, F_max ./ F_ptrain_unsat, 1);
     alpha = 2/pi * ( 1./r .* asin(r) + sqrt(1 - r.^2) );
     f_sat = alpha .* r;
     mult = get_multiplier(f_sat,m_f,B_f,K_f,w, B_f./B_p, K_f./K_p); % fixme this is wrong for multibody

mdocean/simulation/modules/geometry.m

@@ -1,5 +1,5 @@
 function [V_d, m_m, m_f_tot, m_s_tot,...
-         A_c, A_lat_sub, r_over_t, ...
+         A_c, A_lat_sub, ...
          I, T, V_f_pct, V_s_pct, GM, mass,...
          CB_f_from_waterline,CG_f_from_waterline] = geometry(D_s, D_f, D_f_in, D_f_b, ...
                                             T_f_1, T_f_2, h_f, h_s, h_fs_clear, D_f_tu, ...
@@ -144,9 +144,6 @@
 
 A_c = [A_f_cross_top, A_vc_c, A_d_c];
 A_lat_sub = [A_f_l A_vc_l A_d_l];
-r_over_t = [0,... % D_sft/(2*t_sf) 
-            D_s/(2*t_s_r),...
-            0];%D_or/(2*t_r)];
 I = [I_f, I_vc, I_rp];
 T = [T_f_2, T_s, h_d];          % drafts: used to calculated F_surge in dynamics.m
 m_m = m_f_m + m_s_m;            % total mass of material

mdocean/simulation/modules/structures.m

@@ -1,6 +1,21 @@
 function [FOS1Y,FOS2Y,FOS3Y,FOS_buckling] = structures(...
-          	F_heave, F_surge, M, h_s, T_s, rho_w, g, ...
-            sigma_y, A_c, A_lat_sub, r_over_t, I, E)
+          	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_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);
+
+    % 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);
+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)
 
     %% Stress calculations
     depth = [0 T_s T_s]; % max depth
@@ -9,7 +24,7 @@
     sigma_surge = F_surge ./ A_lat_sub;
 
     sigma_rr = P_hydrostatic + sigma_surge;     % radial compression
-    sigma_tt = P_hydrostatic .* r_over_t;       % hoop stress
+    sigma_tt = 0;%P_hydrostatic .* r_over_t;       % hoop stress
     sigma_zz = F_heave ./ A_c;                  % axial compression
     sigma_rt = sigma_surge;                     % shear
     sigma_tz = [0 0 0];
@@ -27,19 +42,65 @@
     sigma_vm = von_mises(sigma_rr, sigma_tt, sigma_zz, sigma_rt, sigma_tz, sigma_zr);
 
     %% Buckling calculation
-    K = 2; % fixed-free - top is fixed by float angular stiffness, bottom is free
-    L = h_s;
-    F_buckling = pi^2 * E(M) * I(2) / (K*L)^2;
-
+    [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_y(M) ./ sigma_vm;
+    FOS_yield = sigma_max ./ sigma_vm;
     FOS1Y = FOS_yield(1);
-    FOS2Y = FOS_yield(2);
+    FOS2Y = FOS_spar;
     FOS3Y = FOS_yield(3);
-    FOS_buckling = F_buckling ./ F_heave;
 
 end 
 
+function [FOS_spar, FOS_spar_local] = spar_combined_buckling(F, E, I, L, D, A, t, q, sigma_0, nu)
+    % euler buckling
+    K = 2; % fixed-free - top is fixed by float angular stiffness, bottom is free
+    F_buckling = pi^2 * E * I / (K*L)^2;
+    sigma_EA = F_buckling / A;
+
+    % hoop stress
+    sigma_theta = q*D/(2*t);
+
+    % local buckling of plate element: 2/9.7
+    k_s = 1.33; % uniform compression, fixed-free, from table 3 on page 30
+    s = D; % fixme - not sure if this is correct
+    P_r = 0.6; % steel proportional linear elastic limit
+    sigma_Ex = k_s * pi^2 * E / (12 * (1-nu^2)) * (t/s)^2;
+    if sigma_Ex <= P_r * sigma_0
+        sigma_Cx = sigma_Ex;
+    else
+        sigma_Cx = sigma_0 * (1 - P_r * (1-P_r) * sigma_0/sigma_Ex);
+    end
+
+    % stress at failure: 2/7.3
+    compact = true; % assumption, eventually this should be a calculation
+    if compact
+        sigma_F = sigma_0; % yield
+    else
+        sigma_F = sigma_Cx; % local buckling of plate element
+    end
+
+    % whether the failure mode is pure euler buckling or combined loading
+    sigma_EA_thresh = P_r * sigma_F * (1 - sigma_theta / sigma_F);
+    pure_buckling = sigma_EA <= sigma_EA_thresh; 
+    if pure_buckling
+        sigma_C_A_theta = sigma_EA;
+    else
+        zeta = 1 - P_r*(1 - P_r)* sigma_F/sigma_EA - sigma_theta/sigma_F;
+        omega = 0.5 * sigma_theta/sigma_F * (1 - .5 * sigma_theta/sigma_F);
+        Lambda = 1/2 * (zeta + sqrt(zeta^2 + 4*omega));
+        sigma_C_A_theta = sigma_F * Lambda;
+    end
+
+    sigma_ac = F ./ A + q;
+    FOS_spar = sigma_C_A_theta / sigma_ac;
+
+    % local buckling of tube: final part of 2/7.3, which combines 2/9.1
+    % and 2/9.5
+    FOS_spar_local = 3; % fixme: need to implement
+end
+
 function s_vm = von_mises(s_11, s_22, s_33, s_12, s_23, s_31)
 
     principal_term = 1/2 * ( (s_11 - s_22).^2 + (s_22 - s_33).^2 + (s_33 - s_11).^2 );

mdocean/simulation/simulation.m

@@ -4,7 +4,7 @@
 
 %% Assemble inputs
 in = p;
-  
+
 in.D_s   = X(1);      % inner diameter of float (m)
 in.D_f   = X(2);      % outer diameter of float (m)
 in.T_f_2 = X(3);      % draft of float (m)
@@ -32,7 +32,7 @@
 
 %% Run modules
 [V_d, m_m, m_f_tot, m_s_tot, ...
- A_c, A_lat_sub, r_over_t, ...
+ A_c, A_lat_sub, ...
  I, T, V_f_pct, V_s_pct, GM] = geometry(in.D_s, in.D_f, in.D_f_in, in.D_f_b, ...
                                         in.T_f_1, in.T_f_2, in.h_f, in.h_s, ...
                                         in.h_fs_clear, in.D_f_tu, in.t_f_t, ...
@@ -43,39 +43,44 @@
 
 m_f_tot = max(m_f_tot,1e-3); % zero out negative mass produced by infeasible inputs
 
-[F_heave_max, F_surge_max, F_ptrain_max, ...
-	    P_var, P_avg_elec, P_matrix_elec, ...
-        X_constraints] = dynamics(in, m_f_tot, m_s_tot, V_d, T);
+[F_heave_storm, F_surge_storm, ...
+ F_heave_op, F_surge_op, F_ptrain_max, ...
+ P_var, P_avg_elec, P_matrix_elec, ...
+                    X_constraints] = dynamics(in, m_f_tot, m_s_tot, V_d, T);
 
-[FOS1Y,FOS2Y,FOS3Y,FOS_buckling] = structures(...
-                                    F_heave_max, F_surge_max,...
-                                    in.M, in.h_s, in.T_s, in.rho_w, in.g, in.sigma_y, A_c, ...
-                                    A_lat_sub, r_over_t, I, in.E);
+[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);
 
 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 = 15+numel(p.JPD);
+num_g = 19+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
 g(3) = V_s_pct;                         % prevent spar too heavy
 g(4) = 1 - V_s_pct;                     % prevent spar too light
 g(5) = GM;                              % pitch stability of float-spar system
-g(6) = FOS1Y / p.FOS_min - 1;           % float survives max force
-g(7) = FOS2Y / p.FOS_min - 1;           % spar survives max force
-g(8) = FOS3Y / p.FOS_min - 1;           % damping plate survives max force
-g(9) = FOS_buckling / p.FOS_min - 1;    % spar survives max force in buckling
-g(10) = P_avg_elec;                     % positive power
+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
 %1 + min(Kp_over_Ks,[],'all');   % spar heave stability (positive effective stiffness)
-g(11) = p.LCOE_max/LCOE - 1;            % prevent more expensive than threshold
-g(12) = F_ptrain_max/in.F_max - 1;      % prevent irrelevant max force -
+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 -
                                         % this constraint should always be active
                                         % and is only required when p.cost_perN = 0.
-g(13) = X_constraints(1);               % prevent float rising above top of spar
-g(14) = X_constraints(2);               % prevent float going below bottom of spar
-g(15) = X_constraints(3);               % float amplitude obeys linear theory
-g(16:end) = X_constraints(4:end);       % prevent rising out of water/slamming
+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);               % float amplitude obeys linear theory
+g(20:end) = X_constraints(4:end);       % prevent rising out of water/slamming
 % fixme: add another X_constraint that h_fs_clear is greater than X_u
 
 criteria = all(~isinf(g)) && all(~isnan(g)) && all(isreal(g));
@@ -85,7 +90,7 @@
 end
 
 if nargout > 4 % if returning extra struct output for validation
-    [~, ~, ~, ~, ~, ~, ~, ~, ~, ~, ~, ~, ...
+    [~, ~, ~, ~, ~, ~, ~, ~, ~, ~, ~, ...
      mass, CB_f,CG_f] = geometry(in.D_s, in.D_f, in.D_f_in, in.D_f_b, ...
                                  in.T_f_1, in.T_f_2, in.h_f, in.h_s, ...
                                  in.h_fs_clear, in.D_f_tu, in.t_f_t, ...
@@ -104,9 +109,9 @@
     val.LCOE = LCOE;
     val.power_avg = P_avg_elec;
     val.power_max = max(P_matrix_elec,[],'all');
-    val.force_heave = F_heave_max;
+    val.force_heave = F_heave_storm;
     val.force_ptrain = F_ptrain_max;
-    val.FOS_b = FOS_buckling;
+    val.FOS_spar = FOS_spar(1);
 	val.c_v = P_var;
     val.B_p = B_p;
     val.X_u = X_u;