diff --git a/dev/structures/damping-plate/BoedoPrantilAnnularPlate.mlx b/dev/structures/damping-plate/BoedoPrantilAnnularPlate.mlx index c02332013..28b8f30f0 100644 Binary files a/dev/structures/damping-plate/BoedoPrantilAnnularPlate.mlx and b/dev/structures/damping-plate/BoedoPrantilAnnularPlate.mlx differ diff --git a/dvc.lock b/dvc.lock index 8588602e2..202394847 100644 --- a/dvc.lock +++ b/dvc.lock @@ -621,17 +621,17 @@ stages: size: 5662 - path: ./mdocean/simulation/modules/structures/structures.m hash: md5 - md5: ab0450a9b4e934c7459fbd1e4cc9b11c - size: 15488 + md5: 9305a5e2bd2a23183327f172eefb17e5 + size: 13360 - path: ./mdocean/simulation/simulation.m hash: md5 - md5: ab665c28e85b454ba399f4d87df0ccd5 - size: 10847 + md5: 4ef3e1de6586d04953ef7e98a72dc7f7 + size: 10848 outs: - path: results/Comparison/intermed.mat hash: md5 - md5: af4aa628a0342b4e0876dec1a0bb97df - size: 739536 + md5: 7eecbd31511cbae8b7c008d0b68c4d8c + size: 739828 postpro-Comparison: cmd: matlab -batch "add_mdocean_path(); obj=Comparison; obj.run_all_from_load();" @@ -763,45 +763,45 @@ stages: size: 5662 - path: ./mdocean/simulation/modules/structures/structures.m hash: md5 - md5: ab0450a9b4e934c7459fbd1e4cc9b11c - size: 15488 + md5: 9305a5e2bd2a23183327f172eefb17e5 + size: 13360 - path: ./mdocean/simulation/simulation.m hash: md5 - md5: ab665c28e85b454ba399f4d87df0ccd5 - size: 10847 + md5: 4ef3e1de6586d04953ef7e98a72dc7f7 + size: 10848 - path: results/Comparison/intermed.mat hash: md5 - md5: af4aa628a0342b4e0876dec1a0bb97df - size: 739536 + md5: 7eecbd31511cbae8b7c008d0b68c4d8c + size: 739828 outs: - path: results/Comparison/comparison_power_matrix.pdf hash: md5 - md5: 94c24e7f1a78ea7772bf8e64b17acc5c - size: 24580 + md5: f055c7574735ad2d17261d891d23de1a + size: 24649 - path: results/Comparison/end.mat hash: md5 - md5: 5da0284f80ebbd58427f98e84a2fdcaa + md5: df8d2aed88af9040188ca7806ade4e68 size: 191 - path: results/Comparison/optimal_design_vars.tex hash: md5 - md5: fce8a85daa819df83f86fc2a8f20ecc2 - size: 922 + md5: 9d2f86c3cf5883756b57a50da688d28a + size: 949 - path: results/Comparison/optimal_outputs.tex hash: md5 - md5: 822413ae8a07f5d72668286564cf571f - size: 2522 + md5: 63034c9707d55c0611e60580e3130b8e + size: 2494 - path: results/Comparison/overlaid_geometry.pdf hash: md5 - md5: 144a1925c3013731989fdaae735b29f4 - size: 6225 + md5: 306e49e55fd56eef0b22953b9c9f8448 + size: 5902 - path: results/Comparison/overlaid_hydro_coeffs.pdf hash: md5 - md5: 152aca3205587c40f90bbfad3afbe114 - size: 17220 + md5: 9489b4ad52bdc4acf6340abb9c66ed50 + size: 17339 - path: results/Comparison/probability_CDF.pdf hash: md5 - md5: 679cee36650b20ec418f6f7818e2d2cd - size: 18160 + md5: ac03a20fc46f3dc29996eff4dd3e8a0f + size: 18401 analysis-Constraints: cmd: matlab -batch "add_mdocean_path(); obj=Constraints; obj.run_analysis();" diff --git a/dvc.yaml b/dvc.yaml index 143ccad9a..a72dd6008 100644 --- a/dvc.yaml +++ b/dvc.yaml @@ -1,14 +1,4 @@ stages: - _check-env-tex: - cmd: calkit check environment --name tex - deps: [] - outs: - - .calkit/env-locks/tex.json: - cache: false - persist: true - always_changed: true - desc: Automatically generated by Calkit. Changes made here will be - overwritten. analysis-AllFigCompare: cmd: matlab -batch "add_mdocean_path(); obj=AllFigCompare; obj.run_analysis();" @@ -1498,12 +1488,11 @@ stages: desc: Automatically generated from the 'move-results' stage in calkit.yaml. Changes made here will be overwritten. build-RE-paper: - cmd: calkit xenv -n tex --no-check -- latexmk -cd -norc - -interaction=nonstopmode -r pubs/renewable-energy-mdo/.latexmkrc - -silent -synctex=1 -pdf pubs/renewable-energy-mdo/mdocean.tex + cmd: calkit latex build -e tex --no-check -r + pubs/renewable-energy-mdo/.latexmkrc + pubs/renewable-energy-mdo/mdocean.tex deps: - pubs/renewable-energy-mdo/mdocean.tex - - pubs/renewable-energy-mdo/.latexmkrc - pubs/renewable-energy-mdo/figs/from-matlab/ - pubs/renewable-energy-mdo/figs/manual/ - pubs/renewable-energy-mdo/sections/ @@ -1511,7 +1500,8 @@ stages: - pubs/renewable-energy-mdo/references.bib - pubs/renewable-energy-mdo/numbers.tex - pubs/elsarticle-num-names.bst - - .calkit/env-locks/tex.json + - pubs/renewable-energy-mdo/.latexmkrc + - .calkit/env-locks/tex outs: - pubs/renewable-energy-mdo/mdocean.pdf - pubs/renewable-energy-mdo/mdocean.synctex.gz: @@ -1523,13 +1513,11 @@ stages: desc: Automatically generated from the 'build-RE-paper' stage in calkit.yaml. Changes made here will be overwritten. build-AOR-paper: - cmd: calkit xenv -n tex --no-check -- latexmk -cd -norc - -interaction=nonstopmode -r - pubs/applied-ocean-research-model/.latexmkrc -silent -synctex=1 - -pdf pubs/applied-ocean-research-model/main.tex + cmd: calkit latex build -e tex --no-check -r + pubs/applied-ocean-research-model/.latexmkrc + pubs/applied-ocean-research-model/main.tex deps: - pubs/applied-ocean-research-model/main.tex - - pubs/applied-ocean-research-model/.latexmkrc - pubs/applied-ocean-research-model/sections/ - pubs/applied-ocean-research-model/figs/from-matlab/ - pubs/applied-ocean-research-model/figs/manual/ @@ -1538,7 +1526,8 @@ stages: - pubs/applied-ocean-research-model/numbers.tex - pubs/elsarticle-num-names.bst - mdocean/simulation/modules/OpenFLASH/pubs/src/meem-appendix.tex - - .calkit/env-locks/tex.json + - pubs/applied-ocean-research-model/.latexmkrc + - .calkit/env-locks/tex outs: - pubs/applied-ocean-research-model/main.pdf - pubs/applied-ocean-research-model/main.synctex.gz: diff --git a/mdocean/plots/matlab_figs/structures/viz_damping_plate.m b/mdocean/plots/matlab_figs/structures/viz_damping_plate.m index 80fe6877e..c94d03076 100644 --- a/mdocean/plots/matlab_figs/structures/viz_damping_plate.m +++ b/mdocean/plots/matlab_figs/structures/viz_damping_plate.m @@ -21,7 +21,7 @@ rho = linspace(lam,1,15); N = 100; - [w_nondim,Mr_nondim] = calc_plate(lam,nu,theta,rho,N); + [w_nondim,Mr_nondim] = concentrated_plate_nondim(lam,nu,theta,rho,N); Mr_dim = Mr_nondim * P/(2*pi); sigma_dim = Mr_nondim * 3/pi * P/h^2; @@ -101,25 +101,36 @@ end function [w_nondim,Mr_nondim] = calc_plate_superposed(a,b,F_heave,nu,theta,rho,N) -% fixme: this assumes that F_tube = F_heave/4, ie that all force goes into -% the tubes and none goes into the bottom of the spar, which is wrong. -% Should actually use compatibility and tube stiffness to get F_tube, as in -% damping_plate_structures() in structures.m. - lam = b/a; - [w_nondim_conc,Mr_nondim_conc] = calc_plate(lam,nu,theta,rho,N); - [w_nondim_dist,Mr_nondim_dist] = plate_distributed(a,b,F_heave,nu,rho); - w_nondim = w_nondim_conc + w_nondim_dist'; - Mr_nondim = Mr_nondim_conc + Mr_nondim_dist'; -end + % fixme: this assumes that F_tube = F_heave/4, ie that all force goes into + % the tubes and none goes into the bottom of the spar, which is wrong. + % Should actually use compatibility and tube stiffness to get F_tube, as in + % damping_plate_structures() in structures.m. + % lam = b/a; + % [w_nondim_conc,Mr_nondim_conc] = calc_plate(lam,nu,theta,rho,N); + % [w_nondim_dist,Mr_nondim_dist] = plate_distributed(a,b,F_heave,nu,rho); + % w_nondim = w_nondim_conc + w_nondim_dist'; + % Mr_nondim = Mr_nondim_conc + Mr_nondim_dist'; + + vb = var_bounds(); -function [w_nondim,Mr_nondim,Mt_nondim] = plate_distributed(a,b,F_heave,nu,rho) - [w_nondim,Mr_nondim,Mt_nondim] = distributed_plate_nondim(a,b,F_heave/4,nu,rho); - w_nondim = -w_nondim; - Mr_nondim = -Mr_nondim; - Mt_nondim = -Mt_nondim; % different sign convention on all three outputs -end + D_d = 2*a; + D_s = 2*b; + [P_hydrostatic,A_dt,theta_dt,h_d,A_c] = deal(0); % not used, so set to 0 + t_d = vb.t_d_nom * 1e-3; + L_dt = (D_d - D_s) / (2*cos(theta_d_tu)); % formula copy-pasted from geometry.m + h_stiff = p.h_over_h1_stiff_d * vb.h1_stiff_d; + width_stiff = p.w_over_h1_stiff_d * vb.h1_stiff_d; + + [sigma_vm,delta_max] = damping_plate_structures(F_heave, D_d, D_s,P_hydrostatic,t_d,A_dt,... + theta_dt, L_dt, h_d, A_c, p.E(1), nu, h_stiff, width_stiff,... + p.D_d_tu, p.t_d_tu, p.num_terms_plate, p.radial_mesh_plate, ... + p.num_stiff_d, true); + + + [Mr_over_F_heave, ... + Mt_over_F_heave, ... + Q_nondim, ... + delta_over_a2Fheave_2piDeq, ... + F_tube_over_F_heave] = combined_plate_nondim(rho, b, a, nu, N, plate_tube_K_ratio); -function [w_nondim,Mr_nondim,abcd] = calc_plate(lam,nu,theta,rho,N) - [w_nondim,Mr_nondim,abcd] = concentrated_plate_nondim(lam,nu,theta,rho,N); - w_nondim = -w_nondim; % different sign convention on w but same convention on M end diff --git a/mdocean/simulation/modules/structures/distributed_plate_nondim.m b/mdocean/simulation/modules/structures/distributed_plate_nondim.m index 9209d355a..014e013a9 100644 --- a/mdocean/simulation/modules/structures/distributed_plate_nondim.m +++ b/mdocean/simulation/modules/structures/distributed_plate_nondim.m @@ -1,51 +1,82 @@ -function [w_nondim,Mr_nondim,Mt_nondim] = distributed_plate_nondim(a,b,F_heave,nu,rho) - A = pi*a^2; - q = F_heave/A; - P = F_heave; +function [w_nondim,Mr_nondim,Mt_nondim,Q_nondim] = distributed_plate_nondim(a,b,q,nu,rho) +% Roark's table 11.2, case 2L (page 467) +% a: outer radius +% b: inner radius +% r0: radius at which the distributed loading begins + +% A = pi*a^2; +% q = F_heave/A; +% P = F_heave; v = nu; r = rho*a; r0 = b; - C2 = 1/4*(1-(b/a)^2*(1+2*log(a/b))); + C2 = 1/4 * ( 1 - (b/a)^2 * (1 + 2*log(a/b)) ); C3 = b/4/a*(((b/a)^2+1)*log(a/b)+(b/a)^2-1); - C8 = 1/2*(1+v+(1-v)*(b/a)^2); - C9 = b/a*((1+v)/2*log(a/b)+(1-v)/4*(1-(b/a)^2)); + C8 = 1/2 * (1+v+(1-v)*(b/a)^2); + C9 = b/a * ( (1+v)/2*log(a/b) + (1-v)/4*(1-(b/a)^2) ); %L3 = r0/4/a*(((r0/a)^2+1)*log(a/r0)+(r0/a)^2-1); % for case 1L %L9 = r0/a*((1+v)/2*log(a/r0)+(1-v)/4*(1-(r0/a)^2)); % for case 1L - L11 = 1/64*(1+4*(r0/a)^2-5*(r0/a)^4-4*(r0/a)^2*(2+(r0/a)^2)*log(a/r0)); % for end deflection only - L17 = 1/4*(1-(1-v)/4*(1-(r0/a)^4)-(r0/a)^2*(1+(1+v)*log(a/r0))); + L11 = 1/64 * (1+4*(r0/a)^2-5*(r0/a)^4-4*(r0/a)^2*(2+(r0/a)^2)*log(a/r0)); % for end deflection only + L17 = 1/4 * (1-(1-v)/4*(1-(r0/a)^4)-(r0/a)^2*(1+(1+v)*log(a/r0))); F2 = 1/4 * (1 - (b./r).^2.*(1+2*log(r/b))); F3 = b./(4*r) .* ( ( (b./r).^2 + 1 ).*log(r/b) + (b./r).^2 - 1); + F5 = 1/2 * (1 - (b./r).^2); + F6 = b./(4*r) .* ( (b./r).^2 - 1 + 2*log(r/b) ); F8 = 1/2 * (1+v+(1-v)*(b./r).^2); F9 = b./r .* (1/2*(1+v)*(log(r/b)) + 1/4*(1-v)*(1-(b./r).^2)); + bracket = zeros(size(r)); bracket(r > r0) = 1; ratio = r0./r; G11 = 1/64 * (1 + 4*ratio.^2 - 5*ratio.^4 - 4*ratio.^2.*(2+ratio.^2).*log(1./ratio) ) .* bracket; + G14 = 1/16 * (1 - (r0./r).^4 - 4*(r0./r).^2 .* log(r/r0)) .* bracket; G17 = 1/4 * (1 - ((1-v)/4)*(1-ratio.^4) - (ratio).^2.*(1+(1+v)*log(1./ratio))) .* bracket; + % Mrb: reaction moment per unit length on inner edge Mrb = -q*a^2/C8 * (C9*(a^2-r0^2)/(2*a*b) - L17); + % Qb: shear reaction force per unit length on inner edge Qb = q/2/b * (a^2 - r0^2); - E = 0;% fixme - only needed for Mt - D = 0; % fixme E*h^3/12/(1-v^2); - only needed for Mt - y_over_D = Mrb * r.^2 .* F2 ... + % case 2 general deflection equation p463, with first two terms + % ommitted since yb and theta_b are zero for case 2L + y_times_D = Mrb * r.^2 .* F2 ... + Qb * r.^3 .* F3 ... - q * r.^4 .* G11; - w_nondim = y_over_D * 2*pi/(P*a^2); - tilt_angle = 0; % fixme use equation for theta on p463 + + theta_times_D = Mrb * r .* F5 ... + + Qb * r.^2 .* F6 ... + - q * r.^3 .* G14; + + Mr = Mrb * F8 ... + + Qb * r .* F9 ... + - q * r.^2 .* G17; + Mt = theta_times_D * (1-nu^2) ./ r ... + + nu * Mr; + Q = Qb * b./r ... + - q./(2*r).*(r.^2 - r0.^2) .* bracket; - Mr = Mrb*F8 + Qb*r.*F9 - q*r.^2.*G17; - Mt = tilt_angle*D*(1-nu^2)./r + nu*Mr; + % Roark's p463 nondimensionalization + K_y = y_times_D / (q * a^4); + K_Mt = Mt / (q * a^2); + K_Mr = Mr / (q * a^2); + K_theta = theta_times_D / (q * a^3); + K_Q = Q / (q * a); + % nondimensionalize to be consistent with concentrated loading plate in Boedo and Prantil + % note: w_nondim = 2*Ky/(1-(b/a)^2) + P = q * pi * (a^2-b^2); % equivalent concentrated load + w_nondim = y_times_D * 2*pi/(P*a^2); Mr_nondim = Mr * 2*pi/P; Mt_nondim = Mt * 2*pi/P; + Q_nondim = Q * 2*pi/(P*a); % use sign convention that positive F_heave results in positive % deflection at outer end and positive moment at inner end w_nondim = -w_nondim; Mr_nondim = -Mr_nondim; Mt_nondim = -Mt_nondim; + Q_nondim = -Q_nondim; end \ No newline at end of file diff --git a/mdocean/simulation/modules/structures/distributed_plate_nondim_ortho.m b/mdocean/simulation/modules/structures/distributed_plate_nondim_ortho.m new file mode 100644 index 000000000..c74535a56 --- /dev/null +++ b/mdocean/simulation/modules/structures/distributed_plate_nondim_ortho.m @@ -0,0 +1,161 @@ +function [w_nondim,Mr_nondim,Mt_nondim,Q_nondim] = distributed_plate_nondim_ortho(a,b,q,nu,rho, ... + Dtheta_over_Dr, Drt_over_Dr) +% Roark's table 11.2, case 2L (page 467) +% a: outer radius +% b: inner radius +% r0: radius at which the distributed loading begins +% Now extended for orthotropic bending rigidity by allowing +% different radial/tangential equivalent thicknesses. +% + + % if orthotropic constants aren't given, assume isotropic (all axes the same) + if nargin < 6 + Dtheta_over_Dr = 1; + end + % approximation for coupling: D_rt = nu * sqrt(Dr * Dtheta), reduces to nu*D for isotropic + if nargin < 7 + Drt_over_Dr = nu*sqrt(Dtheta_over_Dr); + end + + v = nu; + r = rho*a; + r0 = b; + + C8 = C_function(8, a, b, v); + C9 = C_function(9, a, b, v); + + F2 = F_function(2, r, b, v); + F3 = F_function(3, r, b, v); + F5 = F_function(5, r, b, v); + F6 = F_function(6, r, b, v); + F8 = F_function(8, r, b, v); + F9 = F_function(9, r, b, v); + + L17 = L_function(17, r0, a, Dtheta_over_Dr); + + G11 = G_function(11, r0, r, Dtheta_over_Dr); + G14 = G_function(14, r0, r, Dtheta_over_Dr); + G17 = G_function(17, r0, r, Dtheta_over_Dr); + + % Mrb: reaction moment per unit length on inner edge + Mrb = -q*a^2/C8 * (C9*(a^2-r0^2)/(2*a*b) - L17); + % Qb: shear reaction force per unit length on inner edge + Qb = q/2/b * (a^2 - r0^2); + + % ================================================================= + % === ORTHOTROPIC MODIFY: nondimensional rigidity ratios ========== + % ================================================================= + Deff_over_Dr = Dtheta_over_Dr - Drt_over_Dr.^2; % effective rigidity ratio + + % ================================================================= + % === y*D and theta*D (only q-terms need D_eff) ================== + % ================================================================= + + % case 2 general deflection equation p463, with first two terms + % ommitted since yb and theta_b are zero for case 2L + y_times_Dr = Mrb * r.^2 .* F2 ... + + Qb * r.^3 .* F3 ... + - q * r.^4 .* G11; + + % theta_times_D: + % Mrb and Qb contributions use Dr + % q term uses D_eff + theta_times_Dr = Mrb * r .* F5 ... + + Qb * r.^2 .* F6 ... + - q * r.^3 .* (Deff_over_Dr * G14); + + % ================================================================= + % === Mr and Mt =================================================== + % ================================================================= + % Mr depends only on radial rigidity Dr, so unchanged: + Mr = Mrb * F8 ... + + Qb * r .* F9 ... + - q * r.^2 .* G17; + + % Orthotropic tangential moment: + % Mt = (Drt/Dr)*Mr - (D_eff)*(theta/r) + Mt = Drt_over_Dr .* Mr - Deff_over_Dr .* (theta_times_Dr ./ r); + + Q = Qb * b./r ... + - q./(2*r).*(r.^2 - r0.^2) .* bracket; + + % Roark's p463 nondimensionalization + K_y = y_times_D / (q * a^4); + K_Mt = Mt / (q * a^2); + K_Mr = Mr / (q * a^2); + K_theta = theta_times_D / (q * a^3); + K_Q = Q / (q * a); + + % nondimensionalize to be consistent with concentrated loading plate in Boedo and Prantil + % note: w_nondim = 2*Ky/(1-(b/a)^2) + P = q * pi * (a^2-b^2); % equivalent concentrated load + w_nondim = y_times_Dr * 2*pi/(P*a^2); + Mr_nondim = Mr * 2*pi/P; + Mt_nondim = Mt * 2*pi/P; + Q_nondim = Q * 2*pi/(P*a); + + % use sign convention that positive F_heave results in positive + % deflection at outer end and positive moment at inner end + w_nondim = -w_nondim; + Mr_nondim = -Mr_nondim; + Mt_nondim = -Mt_nondim; + Q_nondim = -Q_nondim; +end + +function L = L_function(number, r0, a, Dtheta_over_Dr) + ratio = r0/a; + L = LG_function(number, ratio, Dtheta_over_Dr); +end + +function G = G_function(number, r0, r, Dtheta_over_Dr) + + bracket = zeros(size(r)); + bracket(r > r0) = 1; + + ratio = r0./r; + + G = LG_function(number, ratio, Dtheta_over_Dr) .* bracket; + +end + +function F = F_function(number, r, b, v) + ratio = r/b; + F = CF_function(number, ratio, v); +end + +function C = C_function(number, a, b, v) + ratio = a/b; + C = CF_function(number, ratio, v); +end + +function CF = CF_function(number, ratio, v) + + if number==2 + CF = 1/4 * (1 - (1./ratio).^2.*(1+2*log(ratio))); + elseif number==3 + CF = 1./ratio/(4) .* ( ( (1./ratio).^2 + 1 ).*log(ratio) + (1./ratio).^2 - 1); + elseif number==5 + CF = 1/2 * (1 - (1./ratio).^2); + elseif number==6 + CF = 1./ratio/(4) .* ( (1./ratio).^2 - 1 + 2*log(ratio) ); + elseif number==8 + CF = 1/2 * (1+v+(1-v)*(1./ratio).^2); + elseif number==9 + CF = 1./ratio .* (1/2*(1+v)*(log(ratio)) + 1/4*(1-v)*(1-(1./ratio).^2)); + else + error(['number=' number ' is not implemented']) + end +end +function LG = LG_function(number, ratio, Dtheta_over_Dr) + + if number==11 + LG = 1/(8*(9-Dtheta_over_Dr)) * (1 + 4*ratio.^2 - 5*ratio.^4 - 4*ratio.^2.*(2+ratio.^2).*log(1./ratio) ); + elseif number==14 + LG = 1/(2*(9-Dtheta_over_Dr)) * (1 - (ratio).^4 - 4*(ratio).^2 .* log(1./ratio)); + elseif number==17 + LG = 2/(9-Dtheta_over_Dr) * (1 - ((1-v)/4)*(1-ratio.^4) - (ratio).^2.*(1+(1+v)*log(1./ratio))); + else + error(['number=' number ' is not implemented']) + end + +end diff --git a/mdocean/simulation/modules/structures/structures.m b/mdocean/simulation/modules/structures/structures.m index 5598584dc..4be59d4d3 100644 --- a/mdocean/simulation/modules/structures/structures.m +++ b/mdocean/simulation/modules/structures/structures.m @@ -49,7 +49,7 @@ %% Damping plate plot_on = false; - radial_stress_damping_plate = damping_plate_structures(F_heave, D_d, D_s,P_hydrostatic,t_d,A_dt,... + sigma_damping_plate = damping_plate_structures(F_heave, D_d, D_s,P_hydrostatic,t_d,A_dt,... theta_dt,L_dt,h_d,A_c,E,nu, h_stiff_d,w_stiff_d, ... D_d_tu, t_d_tu, num_terms_plate, radial_mesh_plate, ... num_stiff_d, plot_on); @@ -57,7 +57,7 @@ %% Factor of Safety (FOS) Calculations FOS1Y = sigma_max / sigma_float_bot; FOS2Y = FOS_spar; - FOS3Y = sigma_max / radial_stress_damping_plate; + FOS3Y = sigma_max / sigma_damping_plate; end @@ -117,25 +117,113 @@ a = D_d/2; b = D_s/2; - % nondimensional annular plate solutions rho = linspace(b/a, 1, radial_mesh_plate); % evaluate at these radial points + + % stiffness ratio and dimensional section properties + [plate_tube_K_ratio, h_eq_vec, y_max_vec, D_eq] = get_plate_tube_stiffness(... + h_stiff, num_stiffeners, ... + rho, a, t_d, width_stiff, E, nu, ... + D_d_tu, t_d_tu, ... + L_dt, D_d, D_s); + + % nondimensional moment, displacement, and tube force ratio + [Mr_over_F_heave, ... + Mt_over_F_heave, ... + Q_nondim, ... + delta_over_a2Fheave_2piDeq, ... + F_tube_over_F_heave] = combined_plate_nondim(rho, b, a, nu, N, ... + plate_tube_K_ratio); + + % dimensional forces/moments + F_tube = F_tube_over_F_heave * F_heave; + Mr = Mr_over_F_heave / (2*pi) * F_heave; + Mt = Mt_over_F_heave / (2*pi) * F_heave; + Q = Q_nondim * a / (2*pi) * F_heave; + + % dimensional displacement + delta_total = a^2 * F_heave / (2*pi * D_eq) * delta_over_a2Fheave_2piDeq; + delta_max = max(abs(delta_total)); + + % dimensional stress + sigma_rr_vec = get_plate_stress(Mr, y_max_vec, h_eq_vec); + sigma_tt_vec = get_plate_stress(Mt, t_d/2, t_d); % stiffeners don't affect section properties in theta + sigma_rz_vec = Q / t_d; + + sigma_zz = 0; sigma_rt = 0; sigma_tz = 0; + sigma_vm_vec = von_mises(sigma_rr_vec, sigma_tt_vec, sigma_zz, ... + sigma_rt, sigma_tz, sigma_rz_vec); + + % maximum stress at all radii + sigma_vm = max(abs(sigma_vm_vec)); + + if plot_on + plot_damping_plate(a, D_eq, F_heave, F_tube, r, ... + delta_plate_dis_nondim_vec, delta_plate_con_nondim_vec, ... + Mr_con_nondim, Mr_dis_nondim_vec, Mr, ... + y_max_vec, h_eq_vec, sigma_rr_vec) + end +end + +function [Mr_over_F_heave, ... + Mt_over_F_heave, ... + Q_nondim, ... + delta_over_a2Fheave_2piDeq, ... + F_tube_over_F_heave] = combined_plate_nondim(rho, b, a, nu, N, plate_tube_K_ratio) + + % nondimensional annular plate solutions theta = [0 pi/2 pi 3*pi/2]; % evaluate at these angles - simulating a % single concentrated load applied at theta=0 evaluated at these 4 angles % and summing is the same as simulating 4 concentrated loads applied at % these angles evaluated at theta=0, due to superposition. - [delta_plate_dis_nondim_vec, Mr_dis_nondim_vec, ... - Mt_dis_nondim_vec] = distributed_plate_nondim(a,b,F_heave,nu,rho); - [delta_plate_con_nondim_vec, Mr_con_nondim_vec] = concentrated_plate_nondim(b/a,nu,theta,rho,N); + % distrubuted load + q = 1; % actually q = F_heave / (pi * (a^2-b^2)) but q=1 makes it nondimensional + [delta_plate_dis_nondim_vec, ... % defined as delta_dis * 2pi * D / (F_heave * a^2) evaluated at all r's + Mr_dis_nondim_vec, ... % defined as Mr_dis * 2pi / F_heave evaluated at all r's + Mt_dis_nondim_vec, ... % defined as Mt_dis * 2pi / F_heave evaluated at all r's + Q_dis_nondim_vec] ... % defined as Q_dis * 2pi / (F_heave * a) evaluated at all r's + = distributed_plate_nondim(a,b,q,nu,rho); + + % concentrated load + [delta_plate_con_nondim_vec, ... % defined as delta_con * 2pi * D / (F_tube * a^2) for a single F_tube evaluated at 4 thetas and all r's + Mr_con_nondim_vec] ... % defined as Mr_con * 2pi / (F_tube) for a single F_tube evaluated at 4 thetas and all r's + = concentrated_plate_nondim(b/a,nu,theta,rho,N); + Mt_con_nondim = 0; % solutions aren't available for theta bending moment and shear for concentrated load, so assume zero + Q_con_nondim = 0; % use deflection at outer edge for compatibility delta_plate_dis_nondim = delta_plate_dis_nondim_vec(end); delta_plate_con_nondim = delta_plate_con_nondim_vec(end,:); % sum the four angles for concentrated solution - delta_plate_con_nondim = sum(delta_plate_con_nondim); - Mr_con_nondim = sum(Mr_con_nondim_vec,2); + delta_plate_con_nondim = sum(delta_plate_con_nondim); % total for four equally-spaced F_tubes + Mr_con_nondim = sum(Mr_con_nondim_vec,2); % total for four equally-spaced F_tubes + + % compatibility using tube stiffness (derived from equating plate and tube deflection) + % see notebook 7 p70 1/16/25 and notebook 9 p1 11/20/25 + F_tube_over_F_heave = delta_plate_dis_nondim / (plate_tube_K_ratio - delta_plate_con_nondim); + % nondim moment superposition + Mr_con_over_F_heave = Mr_con_nondim.' * F_tube_over_F_heave; % defined as Mr_con * 2 * pi / F_heave for four equally-spaced F_tubes + Mr_dis_over_F_heave = Mr_dis_nondim_vec; % defined as Mr_dis * 2 * pi / F_heave + Mr_over_F_heave = Mr_con_over_F_heave + Mr_dis_over_F_heave; % defined as Mr * 2 * pi / F_heave for combined loading (dist load F_heave and four equally-spaced point loads F_tube) + + Mt_con_over_F_heave = Mt_con_nondim.' * F_tube_over_F_heave; + Mt_dis_over_F_heave = Mt_dis_nondim_vec; + Mt_over_F_heave = Mt_con_over_F_heave + Mt_dis_over_F_heave; + + % nondim shear force superposition + Q_nondim = Q_con_nondim.' * F_tube_over_F_heave + Q_dis_nondim_vec; + + % nondim displacement superposition + delta_over_a2Fheave_2piDeq = delta_plate_dis_nondim_vec ... + + F_tube_over_F_heave * sum(delta_plate_con_nondim_vec.',1); % defined as delta * 2pi * D / (F_heave * a^2) for combined loading +end + +function [plate_tube_K_ratio, h_eq_vec, y_max_vec, D_eq] = get_plate_tube_stiffness(h_stiff,num_stiffeners, ... + rho, a, t_d, width_stiff, E, nu, ... + D_d_tu, t_d_tu, ... + L_dt, D_d, D_s) % stiffeners num_unique_stiffeners = length(h_stiff)/2; num_stiffener_repeats = num_stiffeners / num_unique_stiffeners; @@ -154,31 +242,8 @@ I_tube = pi/64 * ( D_d_tu^4 - (D_d_tu - 2*t_d_tu)^4 ); K_tube = 6*E*I_tube / (L_dt^2 * (D_d - D_s)); - % compatbility (derived from equating plate and tube deflection) - F_tube = F_heave * delta_plate_dis_nondim / (D_eq/(a^2*K_tube) - delta_plate_con_nondim); - - % moment superposition - Mr_con = Mr_con_nondim.' * F_tube; - Mr_dis = Mr_dis_nondim_vec * F_heave; - Mr = Mr_dis + Mr_con; - - % displacement superposition - delta_total = a^2/D_eq * (F_heave * delta_plate_dis_nondim_vec ... - + F_tube * sum(delta_plate_con_nondim_vec.',1)); - delta_max = max(abs(delta_total)); - - % stress - sigma_r_vec = get_plate_stress(Mr, y_max_vec, h_eq_vec); - sigma_r = max(abs(sigma_r_vec)); - - sigma_vm = sigma_r; % ignore Mt for now - - if plot_on - plot_damping_plate(a, D_eq, F_heave, F_tube, r, ... - delta_plate_dis_nondim_vec, delta_plate_con_nondim_vec, ... - Mr_con_nondim, Mr_dis_nondim_vec, Mr, ... - y_max_vec, h_eq_vec, sigma_r_vec) - end + % notebook p1 11/20/25 + plate_tube_K_ratio = D_eq*2*pi/(a^2*K_tube); end function [] = plot_damping_plate(r, delta_total, ... diff --git a/mdocean/simulation/simulation.m b/mdocean/simulation/simulation.m index a33f936fb..2e8661c3c 100644 --- a/mdocean/simulation/simulation.m +++ b/mdocean/simulation/simulation.m @@ -78,7 +78,7 @@ 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_f, ... % bulk dimensions - in.D_f_tu, in.D_d, L_dt, in.theta_d_tu,in.D_d_tu,... % more bulk dimensions + in.D_f_tu, in.D_d, L_dt, in.theta_d_tu, in.D_d_tu,... % more bulk dimensions in.t_s_r, I, A_c, A_lat_sub, in.t_f_b, in.t_f_t, in.t_d, in.t_d_tu, in.h_d, A_dt, ... % structural dimensions in.h_stiff_f, in.w_stiff_f, in.h_stiff_d, in.w_stiff_d,... % stiffener thicknesses in.M, in.rho_w, in.g, in.sigma_y, in.sigma_e, in.E, in.nu, ... % constants