Damping plate stress

Roark_func.m

@@ -0,0 +1,101 @@
+function [y,l1,l2] = Roark_func(r0_point,r0_dist,a,w,t,E,q,b,v)
+%Roark_func
+%    [Y,L1,L2] = Roark_func(R0_POINT,R0_DIST,A,W,T,E,Q,B,V)
+
+%    This function was generated by the Symbolic Math Toolbox version 24.1.
+%    04-Dec-2024 15:29:33
+
+t2 = a.^2;
+t3 = a.^3;
+t5 = b.^2;
+t6 = r0_dist.^2;
+t8 = r0_point.^2;
+t9 = v+1.0;
+t10 = v.^2;
+t11 = 1.0./E;
+t12 = 1.0./a;
+t15 = 1.0./b;
+t16 = 1.0./r0_dist;
+t17 = 1.0./r0_point;
+t18 = 1.0./t.^3;
+t19 = v-1.0;
+t20 = v./2.0;
+t21 = v./4.0;
+t4 = t2.^2;
+t7 = t6.^2;
+t13 = 1.0./t2;
+t22 = -t6;
+t23 = t10-1.0;
+t24 = a.*t15;
+t25 = a.*t16;
+t26 = a.*t17;
+t34 = t20+1.0./2.0;
+t36 = t21-1.0./4.0;
+t14 = 1.0./t4;
+t27 = log(t24);
+t28 = log(t25);
+t29 = log(t26);
+t30 = t5.*t13;
+t31 = t6.*t13;
+t33 = t8.*t13;
+t38 = t2+t22;
+t32 = t7.*t14;
+t35 = t27.*2.0;
+t37 = t30+1.0;
+t39 = t31+2.0;
+t40 = t33+1.0;
+t41 = t9.*t28;
+t42 = t30-1.0;
+t44 = t33-1.0;
+t46 = t31./1.6e+1;
+t50 = (t19.*t30)./2.0;
+t52 = t27.*t34;
+t53 = t29.*t34;
+t43 = t32-1.0;
+t45 = t35+1.0;
+t47 = t41+1.0;
+t48 = t32.*(5.0./6.4e+1);
+t51 = -t50;
+t54 = t27.*t37;
+t55 = t29.*t40;
+t57 = t36.*t42;
+t58 = t36.*t44;
+t67 = t28.*t31.*t39.*(-1.0./1.6e+1);
+t49 = -t48;
+t56 = (t30.*t45)./4.0;
+t59 = (t31.*t47)./4.0;
+t60 = (t36.*t43)./4.0;
+t61 = t34+t51;
+t65 = t42+t54;
+t66 = t44+t55;
+t71 = t52+t57;
+t72 = t53+t58;
+t62 = t56-1.0./4.0;
+t63 = 1.0./t61;
+t68 = (r0_point.*t12.*t65)./4.0;
+t69 = (r0_point.*t12.*t66)./4.0;
+t73 = r0_point.*t12.*t71;
+t74 = r0_point.*t12.*t72;
+t76 = (t13.*t38.*t71)./2.0;
+t77 = q.*t2.*t11.*t18.*t23.*t38.*t65.*(3.0./2.0);
+t79 = t46+t49+t67+1.0./6.4e+1;
+t70 = -t69;
+t75 = -t74;
+t78 = -t77;
+t80 = q.*t4.*t11.*t18.*t23.*t79.*1.2e+1;
+t82 = t59+t60+t76-1.0./4.0;
+t81 = t73+t75;
+t84 = q.*t4.*t11.*t18.*t23.*t62.*t63.*t82.*1.2e+1;
+t83 = t62.*t63.*t81;
+t85 = -t84;
+t86 = t68+t70+t83;
+t87 = t3.*t11.*t18.*t23.*t86.*w.*1.2e+1;
+t88 = -t87;
+y = t78+t80+t85+t88;
+if nargout > 1
+    l1 = t88;
+end
+if nargout > 2
+    l2 = t78+t80+t85;
+end
+end

damping_plate_stress.m

@@ -0,0 +1,85 @@
+    syms D_d A_dt theta_dt L_dt h_d D_s t_d
+    syms A_c [1,3]
+    syms P_hydrostatic [1,3]
+    syms E [1,3]
+
+    %What codey thing dictates which material we're using so we don't have
+    %to hardcode this?
+    E = E(1);
+
+    %calculate the deflection of the damping plate
+    r = D_d/2;
+    r_bending = r-(D_s/2);
+    F_water = P_hydrostatic(2) * A_c(3)/2;
+    syms F_supportsym xsym
+    Msym = -(F_water+F_supportsym*sin(theta_dt))*xsym + (F_water*xsym^2/(2*r_bending)) + ...
+        (F_water*r_bending/2) + (F_supportsym*r_bending*sin(theta_dt));
+    Isym = (1/6) * sqrt(r_bending^2 - xsym^2) * (h_d^3 - t_d^3);
+    y_double_prime = Msym/(E*Isym);
+    y_prime = int(y_double_prime,xsym);
+    y = int(y_prime,xsym);
+    y_tip = subs(y,xsym,r_bending-0.5);
+    eqn = F_supportsym == y_tip * E * A_dt / L_dt;
+    F_support = solve(eqn,F_supportsym);
+
+    I_d = (1/12) * r_bending * (h_d^3 - t_d^3);
+    %F_support = 3*F_water*r_bending^3*A_dt/(8*sin(theta_dt)*(3*L_dt*I_d - (r_bending^3*A_dt)));
+
+    % x1 = linspace(0,r-t_d,1500);
+    % A = 2*sqrt(r^2-x1.^2) * h_d - (2*sqrt(r^2-x1.^2) * t_d);
+    % sigma_axial = F_support * cos(theta_dt) ./ A;
+    % v = abs(-F_water - (F_support*(sin(theta_dt))) + (F_water*x1/r));
+    % shear = max(v./A);
+
+    x = linspace(0,r_bending-0.01,1500);
+    x1 = linspace(0,r-t_d,1500);
+
+    A = 2*sqrt(r_bending^2-x.^2) * h_d - (2*sqrt(r_bending^2-x.^2) * t_d);
+    sigma_axial = F_support * cos(theta_dt) ./ A;
+    v = -F_water - (F_support*(sin(theta_dt))) + (F_water*x1/r_bending);
+    shear = max(v./A);
+
+    M_x = -(F_water+F_support*sin(theta_dt))*x + (F_water.*x.^2/(2*r_bending)) + ...
+        (F_water*r_bending/2) + (F_support*r_bending*sin(theta_dt));
+    I_x = (1/12) * 2*sqrt(r_bending^2-x.^2) * h_d^3 - (((1/12) * 2*sqrt(r_bending^2-x.^2) * t_d^3));
+    sigma_bending = abs(M_x .* h_d./(2.*I_x));
+    sigma_xx = max(sigma_bending+sigma_axial);
+
+    %debugging
+    % figure
+    % hold on
+    % title("Bending and Axial Stress vs. Radial Position")
+    % plot(x,sigma_bending)
+    % plot(x,sigma_axial)
+    % legend("bending","axial")
+    % xlabel("x")
+    % ylabel("stress")
+    % hold off
+    % 
+    % figure
+    % hold on
+    % title("Shear and Bending Moment Diagram")
+    % plot(x,v)
+    % plot(x,M_x)
+    % legend("shear","bending moment")
+    % xlabel("x")
+    % hold off
+    % 
+    % figure
+    % hold on
+    % title("Moment of Inertia vs. Radial Position")
+    % plot(x,I_x)
+    % legend("moment of inertia")
+    % xlabel("x")
+    % ylabel("I")
+    % hold off
+    % 
+    % figure
+    % hold on
+    % plot(x,sigma_axial)
+    % title("Axial Stress vs. Radial Position")
+    % xlabel("x")
+    % ylabel("Axial Stress")
+    % hold off
+
+    ht = matlabFunction(sigma_xx, shear, 'vars', [E, D_d, A_dt, theta_dt, L_dt, h_d, D_s, t_d, A_c, P_hydrostatic])

mdocean/simulation/modules/damping_plate_func.m

@@ -0,0 +1,147 @@
+function [sigma_xx,shear,sigma_bending,sigma_axial,M_x,I_x,F_water,F_support,y] = damping_plate_func(E1,D_d,A_dt,theta_dt,L_dt,h_d,D_s,t_d,A_c1,A_c2,A_c3,P_hydrostatic1,P_hydrostatic2,P_hydrostatic3,x,x1,F_heave)
+%DAMPING_PLATE_FUNC
+%    [SIGMA_XX,SHEAR,SIGMA_BENDING,SIGMA_AXIAL,M_x,I_x,F_water,F_support,Y] = DAMPING_PLATE_FUNC(E1,D_d,A_dt,THETA_DT,L_dt,H_D,D_s,T_D,A_c1,A_c2,A_c3,P_hydrostatic1,P_hydrostatic2,P_hydrostatic3,X,X1,F_heave)
+
+%    This function was generated by the Symbolic Math Toolbox version 24.1.
+%    25-Nov-2024 15:16:54
+
+t2 = cos(theta_dt);
+t3 = sin(theta_dt);
+t4 = D_d.^2;
+t5 = D_s.^2;
+t6 = h_d.^3;
+t7 = t_d.^3;
+t8 = x.^2;
+t9 = -D_s;
+t10 = 1.0./L_dt;
+t11 = D_d./2.0;
+t12 = D_s./2.0;
+t13 = F_heave./2.0;
+t14 = D_d.*D_s.*F_heave.*1.8e+1;
+t26 = D_d.*D_s.*(-1.0./2.0);
+t15 = -t8;
+t16 = D_d+t9;
+t17 = -t11;
+t18 = -t12;
+t20 = t4./4.0;
+t21 = t5./4.0;
+t22 = -t14;
+t23 = F_heave.*t4.*9.0;
+t24 = F_heave.*t5.*9.0;
+t25 = D_d.*D_s.*t3.*4.8e+1;
+t27 = D_d.*E1.*t6.*8.0;
+t28 = D_s.*E1.*t6.*8.0;
+t29 = D_d.*E1.*t7.*8.0;
+t30 = D_s.*E1.*t7.*8.0;
+t32 = t3.*t4.*2.4e+1;
+t33 = t3.*t5.*2.4e+1;
+t31 = -t25;
+t34 = 1.0./t16;
+t35 = -t28;
+t36 = -t29;
+t37 = t11+t18;
+t38 = t12+t17+1.0;
+t48 = t22+t23+t24;
+t49 = t20+t21+t26;
+t39 = t37.^2;
+t40 = t38.^2;
+t41 = t8.*t13.*t34;
+t42 = (F_heave.*t37)./4.0;
+t50 = t31+t32+t33;
+t52 = sqrt(t49);
+t63 = t27+t30+t35+t36;
+t43 = -t40;
+t44 = sqrt(t39);
+t47 = t15+t39;
+t67 = 1.0./t63;
+t70 = F_heave.*t16.*t38.*t52.*2.4e+1;
+t71 = t3.*t16.*t38.*t52.*4.8e+1;
+t45 = t44.^3;
+t51 = sqrt(t47);
+t62 = t39+t43;
+t69 = t44.*t48;
+t89 = t44.*t50.*t67;
+t46 = F_heave.*t45.*4.0;
+t53 = 1.0./t51;
+t54 = h_d.*t51.*2.0;
+t55 = t51.*t_d.*2.0;
+t59 = (t6.*t51)./3.0;
+t60 = (t7.*t51)./3.0;
+t64 = sqrt(t62);
+t56 = t53.*x;
+t58 = -t55;
+t61 = -t60;
+t65 = t64.^3;
+t66 = 1.0./t64;
+t75 = t48.*t64;
+t76 = F_heave.*t16.*t38.*t64.*1.2e+1;
+t78 = t3.*t16.*t38.*t64.*2.4e+1;
+t81 = t46+t69;
+t82 = t50.*t64;
+t57 = atan(t56);
+t68 = F_heave.*t65.*4.0;
+t72 = t38.*t66;
+t74 = t54+t58;
+t79 = -t76;
+t80 = -t78;
+t83 = t59+t61;
+t97 = t67.*t81;
+t73 = atan(t72);
+t77 = 1.0./t74;
+t84 = 1.0./t83;
+t85 = F_heave.*t16.*t40.*t73.*1.2e+1;
+t86 = t3.*t16.*t40.*t73.*2.4e+1;
+t91 = t38.*t48.*t73;
+t92 = t38.*t50.*t73;
+t93 = F_heave.*t16.*t62.*t73.*1.2e+1;
+t95 = t3.*t16.*t62.*t73.*2.4e+1;
+t87 = -t85;
+t88 = -t86;
+t94 = -t93;
+t96 = -t95;
+t98 = t71+t80+t82+t88+t92+t96;
+t99 = t68+t70+t75+t79+t87+t91+t94;
+t100 = t67.*t98;
+t101 = t67.*t99;
+t102 = -t101;
+t106 = -1.0./(A_dt.*E1.*t10.*(t89-t100)+1.0);
+t107 = t97+t102;
+t108 = A_dt.*E1.*t3.*t10.*t106.*t107;
+t109 = (A_dt.*E1.*t3.*t10.*t107.*-2.0)./(A_dt.*E1.*t10.*(t89-t100)+1.0);
+t112 = (A_dt.*D_d.*D_s.*E1.*t3.*t10.*t107.*4.8e+1)./(A_dt.*E1.*t10.*(t89-t100)+1.0);
+t119 = A_dt.*E1.*t2.*t10.*t77.*t106.*t107;
+t110 = F_heave+t109;
+t115 = t13+t108;
+t118 = t37.*t108;
+t116 = t115.*x;
+t117 = -t116;
+t120 = t41+t42+t117+t118;
+t121 = h_d.*t84.*t120;
+sigma_xx = t119+t121;
+if nargout > 1
+    shear = [t77.*(-t13+(t13.*x1)./t37+(A_dt.*E1.*t3.*t10.*t107)./(A_dt.*E1.*t10.*(t89-t100)+1.0))];
+end
+if nargout > 2
+    sigma_bending = t121;
+end
+if nargout > 3
+    sigma_axial = t119;
+end
+if nargout > 4
+    M_x = t120;
+end
+if nargout > 5
+    I_x = (t6.*t51)./6.0-(t7.*t51)./6.0;
+end
+if nargout > 6
+    F_water = t13;
+end
+if nargout > 7
+    F_support = A_dt.*E1.*t10.*t106.*t107;
+end
+if nargout > 8
+    et1 = -t67.*(t51.*(t14-t23-t24-t112+(A_dt.*E1.*t10.*t32.*t107)./(A_dt.*E1.*t10.*(t89-t100)+1.0)+(A_dt.*E1.*t10.*t33.*t107)./(A_dt.*E1.*t10.*(t89-t100)+1.0))-F_heave.*1.0./t53.^3.*4.0+t16.*x.*(F_heave.*t52+t52.*t109).*2.4e+1+t57.*x.*(t14-t23-t24-t112+(A_dt.*E1.*t10.*t32.*t107)./(A_dt.*E1.*t10.*(t89-t100)+1.0)+(A_dt.*E1.*t10.*t33.*t107)./(A_dt.*E1.*t10.*(t89-t100)+1.0))-t8.*t16.*t57.*t110.*1.2e+1-t16.*t51.*t110.*x.*1.2e+1+t16.*t57.*t110.*(t8-t39).*1.2e+1);
+    et2 = -t67.*(t46-t44.*(t14-t23-t24-t112+(A_dt.*E1.*t10.*t32.*t107)./(A_dt.*E1.*t10.*(t89-t100)+1.0)+(A_dt.*E1.*t10.*t33.*t107)./(A_dt.*E1.*t10.*(t89-t100)+1.0)));
+    y = et1+et2;
+end
+end

mdocean/simulation/modules/damping_plate_stress.m

@@ -0,0 +1,56 @@
+syms D_d A_dt theta_dt L_dt h_d D_s t_d F_heave
+syms A_c [1,3]
+syms P_hydrostatic [1,3]
+syms E [1,3]
+syms x x1 c1 c2
+
+%What codey thing dictates which material we're using so we don't have
+%to hardcode this?
+E = E(1);
+
+%set up the dimensions and pressures
+r = D_d/2;
+r_bending = r-(D_s/2);
+%F_water = P_hydrostatic(2) * A_c(3)/2;
+P_hydrodynamic = F_heave / A_c(3);
+F_water = P_hydrodynamic * A_c(3) / 2;
+
+%solve y" = M/EI to find the deflection of the damping plate
+syms F_supportsym xsym
+Msym = -(F_water+F_supportsym*sin(theta_dt))*xsym + (F_water*xsym^2/(2*r_bending)) + ...
+    (F_water*r_bending/2) + (F_supportsym*r_bending*sin(theta_dt));
+Isym = (1/6) * sqrt(r_bending^2 - xsym^2) * (h_d^3 - t_d^3);
+y_double_prime = Msym/(E*Isym);
+y_prime = int(y_double_prime,xsym) + c1;
+find_c1 = 0 == subs(y_prime,xsym,0);
+c1_solved = solve(find_c1,c1);
+y_prime = subs(y_prime, c1, c1_solved);
+y = int(y_prime,xsym) + c2;
+find_c2 = 0 == subs(y, xsym, 0);
+c2_solved = solve(find_c2,c2);
+y = subs(y,c2,c2_solved);
+
+%use the tip deflection (which is equal to the deflection of the tubular
+%support) to find the force of the tubular support on the damping plate
+y_tip = subs(y,xsym,r_bending-1);
+eqn = F_supportsym == y_tip * E * A_dt / L_dt;
+F_support = solve(eqn,F_supportsym);
+
+%solve for deflection as a function of x for debugging
+y = subs(y, F_supportsym, F_support);
+y = subs(y,xsym,x);
+
+%solve for the shear and axial stresses
+A = 2*sqrt(r_bending^2-x^2) * h_d - (2*sqrt(r_bending^2-x^2) * t_d);
+sigma_axial = F_support * cos(theta_dt) ./ A;
+v = -F_water - (F_support*(sin(theta_dt))) + (F_water*x1/r_bending);
+shear = max(v./A);
+
+%solve for the bending stress
+M_x = -(F_water+F_support*sin(theta_dt))*x + (F_water.*x^2/(2*r_bending)) + ...
+    (F_water*r_bending/2) + (F_support*r_bending*sin(theta_dt));
+I_x = (1/12) * 2*sqrt(r_bending^2-x^2) * h_d^3 - (((1/12) * 2*sqrt(r_bending^2-x^2) * t_d^3));
+sigma_bending = M_x .* h_d./(2.*I_x);
+sigma_xx = sigma_bending+sigma_axial;
+
+matlabFunction(sigma_xx, shear, sigma_bending, sigma_axial, M_x,I_x, F_water, F_support, y, 'File','damping_plate_func','Vars', [E, D_d, A_dt, theta_dt, L_dt, h_d, D_s, t_d, A_c, P_hydrostatic,x,x1,F_heave]);