Produce plots depending on the conditions that are listed within the code

%% Input parameters
A = 264000; %Rectangular cross-sectional area [mm^2]
d_max = 1250; %maximum value for depth of cross section of beam [mm]
b_min = 200; %minimum value for width of cross section of beam [mm]
d_min = sqrt(A);
d = b_min:d_max; %Array of values for d [mm]
b = b_min:d_max;
F = 300e3; %Transverse downwards force [N]
T = 20e7; % Horizontal twisting moment (Nmm)
sigma_adm = 21; %maximum admissible normal stress of beam [MPa]
tau_adm = 3.5; %maximum admissible shear stress of beam [MPa]
AB = 5000; %Length of horizontal beam AB [mm]
BC = 1650; %Length of horizontal beam BC [mm]
AC = AB + BC; %Length of full horizontal beam AC [mm]

%% Support reactions in statically determined beam (Checks for hand calculations)
Vb = ((AB+BC)*F)/AB; %Vertical reaction at support B [N]
Va = F - Vb; %Vertical reaction at support A [N]

%% Alpha values and interpolate in required range
aspect_ratios = [1.0, 1.2, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 10.0];
alpha_tabled = [0.208, 0.219, 0.231, 0.246, 0.258, 0.267, 0.282, 0.291, 0.299, 0.312];
n_pts = 100;                          % number of points in the plots
Delta_d = (d_max-d_min) / (n_pts-1);  % increment in the section's depth
list_d = (d_min : Delta_d : d_max);   % list of cross sections's depths
list_alpha = interp1(aspect_ratios, alpha_tabled, d / d_max);
figure
hold on
plot (aspect_ratios,alpha_tabled,'k-o','LineWidth',2)               % black line in Figure 3
plot ([1,10],[1/3,1/3],'r--','LineWidth',1)                         % red dashed line in Figure 3
plot (d / d_max, list_alpha, 'r*')                                 % plot red * where values are interpolated
title ('Torsion coefficients for rectangular cross section')
xlabel('d/b');
ylabel('\alpha');
hold off
legend('Interpolation','Asymptotic value','Location','southeast')
grid on

%% loop to determine normal and shear stress as functions of beams depth (A must remain constant)
NL = length(d);
sigma_max = zeros(1,NL); %Initialise array for max sigma
tau_max = zeros(1,NL); %Initialise array for max tau
for i = 1:NL
    if ((d(i)>=b(i)) || (d(i)<=d_max)) && ((b(i)>=b_min) || (b(i)*d(i)==A))
        My(i) = F*1650;
        Iyy(i) = (b(i)*d(i)^3)/12;
        Qy(i) = (d(i)*b(i)^2)/8;
        alpha = list_alpha(i);
        J(i) = alpha*d(i)*b(i);
        sigma_max(i) = (My(i)/Iyy(i))*(d(i)/2);
        tauF(i) = (F/Iyy(i)*b(i))*Qy(i);
        tauT(i) = (T/J(i))*b(i);
        tau_max(i) = tauF(i) + tauT(i);
    end
end

%% Plot for normal stress
figure
plot(d,sigma_max,'r')
xlim([min(d) max(d)])
title('Maximum Normal Stress as functions of the beams depth d')
xlabel('Depth (mm)')
ylabel('Maximum Normal Stress (MPa)')

%% Plot for shear stress
figure
plot(d,tau_max,'b')
xlim([min(d) max(d)])
title('Maximum Shear Stress as functions of the beams depth d')
xlabel('Depth (mm)')
ylabel('Maximum Shear Stress (MPa)')

%% Task B
%% Plot for Parametric relationship between sigma and tau
figure
plot(sigma_max,tau_max)
for i = 1:NL
    hold on
    if(tau_max(i) <= tau_adm && sigma_max(i) <= sigma_adm)
        plot(sigma_max(i),tau_max(i),"r*");
    end
    hold off
end
title('Parametric relationship between sigma max (d) and tau max (d)')
xlabel('Sigma max (d)')
ylabel('Tau max (d)')
legend('Parametric relationship','Admissible design solutions')

Warning: Ignoring extra legend entries.

[Luca Ferro]

what is the problem?

[Zara]

When I press run my graph for shear stress won't plot and I am not sure why this is

Answer 1

[Cris LaPierre]

Inspect your variable values to see if they are what you are expecting. The reason your last 2 plots are not appearing is because the values are all NaNs. This suggests that one of your calculations is not doing what you expect, or your data is not what you expected it to be.

You might find the debugging tools in MATLAB helpful for this.

Here, the root issue is that you are using interp1 to extrapolate values when creating list_alpha. By default, it won't do that, so you get NaNs. You must explicitely tell interp1 you want to extrapolate.

See this example:

x = [1 2 3 4 5];
v = [12 16 31 10 6];

Specify the query points, xq, that extend beyond the domain of x.

xq = [0 0.5 1.5 5.5 6];

Evaluate v at xq using the default method. Note that everything outside the range of values in x are NaN

vq1 = interp1(x,v,xq)

vq1 = 1×5

NaN NaN 14 NaN NaN

Now, use the 'linear' method with the 'extrap' option.

vq3 = interp1(x,v,xq,'linear','extrap')

vq3 = 1×5

8 10 14 4 2

You can also use a method that supports extrapolation, like 'pchip'

vq1 = interp1(x,v,xq,'pchip')

vq1 = 1×5

19.3684 13.6316 13.2105 7.4800 12.5600