In Materials and Mechanics (EGR 211), it is crucial to learn and discuss how steel beams bend and deflect when a load is acting on it. It is necessary for engineers to understand how materials react under forces in order to design the proper support for the presented scenario and ensure public safety. In this analysis, we explore two methods of calculating the deflection of a 14m long beam with the following parameters: elastic modulus, moment of inertia, uniform load, and the length of the beam. The deflection of a beam can be calculated analytically as seen in the Alternative Solution Steps section; however, another method is the application of a numerical approach or Euler's method. Since the deflection of a beam changes depending on its length, it is important to implement an analytical method providing an accurate point of comparison when compared to the numerical method. After establishing both methods to calculate the deflection of the beam, it is necessary to perform many iterations in order to obtain an accurate result. For this reason, we take advantage of algorithms created in MATLAB to acquire a precise outcome. Furthermore, after calculating the deflection with each method, these values are compared to the limited safe deflection level presented by L/360 (where L is the length of the beam). It was noted that the uniform load placed on the beam exceeded the allowed amount of weight the beam could support. In effect to the issue, calculations were performed to reveal the maximum amount of stress a uniform load could weigh according to the characteristics of the beam.
There are two initial assumptions we can make about the beam: when the distance (x) = 0m, the deflection (y) = 0m, and when the distance (x) = 7m, the slope of the deflection (dy/dx) = 0. We consider this assumption to be an initial condition and we use this information to discover relationships between variables. We are given the following properties of the beam:
Table 1: Parameter values for the beam
An important first step in solving this sort of problem is to create a free body diagram. Below, the beam is drawn with all forces we consider in our solution:
Figure 1: Free body diagram of the beam1
Now that we have determined our beam parameters and drawn the free body diagram, we can begin our numerical solution steps.
As denoted above, an alternate way to solve the ordinary differential equation (ODE) presented by the calculation of deflection of a beam depending on its length and uniform load, is Euler’s Method. Euler’s method is a first-order method, which means that the local error is proportional to the square of the step size and the global error is proportional to the step size. In this specific case, Euler’s method looks at the change in deflection depending on its length (dy/dx) as seen in equation 1.
(1)
The change of deflection in accordance to length can be associated with a relationship of the modulus of elasticity (E), the moment of inertia (I), the presented uniform load (w), the length of the beam (L), and the specific length at which the deflection is being calculated (a range between x = 0 to x = L) as seen to the right in equation 2:
(2)
Combining equation one and equation two a direct numerical relationship can be obtained with a step size (h) as seen in equation 3:
(3)
Equation 3 can be rearranged to obtained a direct relationship solving for the deflection of the beam in accordance to its length and parameters presented in Table 1. This relationship can be seen in equation 4:
(4)
Note that Euler’s method utilizes a step size (h) which is a chosen arbitrary value that can increase accuracy. In other words, in this case, the lower the step size the more accurate the calculations will be since there will be more values ranging from 0 to L. In this case, two step sizes were chosen: h = 1.0m & h = 0.25m. An example calculation is shown here:
The safety and efficiency of beams can be determined by dividing the length of the beam by a factor of safety of 360 (L/360). In this case, the length of the beam is 14 meters long, which when divided by 360 equals 0.0388. When looking at Figure 2, it is clear that any deflection after 2 meters will be unsustainable by the provided beam. This means that the beam can be subject to serious damage if the uniform load is kept on it. When looking at the equation for maximum deflection allowed, derived in the Analytical Method section, the maximum uniform load allowed on this beam is 5,270 N/m as shown on the right:
The numerical beam analysis method is very similar to the analytical method, with the exception that the uniform load is a numerical approximation depending on a range of values represented by 'n' and an accuracy factor, which in this case is 10. As n increases, w(n) decreases, approaching the limit of the theoretical max change in deflection. In other words, when looking at the example on the right, n is in a chosen range of 0 - 1,000. When n is 523, w(n) is 5270 N/m, which yields a value of approximately 0.03876 of max deflection, which is less than the theoretical allowable max deflection.
Now that we have completed the numerical solution steps, let's move on to the analytical solution steps before exploring the pseudocode for both methods.
An alternative solution to the problem is an analytical solution that also utilizes the change of deflection over distance and their relationship to the physical characteristics of the beam (parameters found on in table 1). This approach begins by looking at the relationship between the radius of curvature (P) and the curvature of the beam (K). This relationship with combination of the calculation of moments along a beam and the physical characteristics of a beam leads to a analytical ODE as seen to the right and below.
This method attempts to derive an equation for deflection as a function of ‘x.’ Note that for small deflections and angles 𝑑𝑠 ≈ 𝑑𝑥
Since a beam is supposed to be stationary it can be assumed that the sum of the moments is equal to zero, therefore the following relationship is defined:
Since 𝐸𝐼𝑦′′ = 𝑀, the previous equation can be arranged into the following relationship:
This relationship provides us with a direct way of solving for deflection after we divide by EI and integrate twice with respect to x:
The first integration gives us the basis equation to use and solve the problem numerically as seen in the above section; however, we want to keep on solving this problem analytically so we must solve for C1 and integrate again. Since x=L/2, ∆′= 0 because ∆′= 𝑑∆/dx:
Integrating again with the consideration that C2 is equal to zero as it is the derivative of another constant:
This is the final analytical solution to the problem when solving for the constants. In addition, from this relationship we can obtain the maximum level of deflection supported according to the characteristics of the beam as seen to the right:
As denoted above, the safety and efficiency of beams can be determined by dividing the length of the beam by a factor of safety of 360 (L/360). With a maximum length of 14m the deflection limit is 0.0388. From Figure 2, it can be easily seen on the analytical curve that this limit is also exceeded right after a length of two meters. This means the beam is not safe to support the presented load and can be subject to serious damage if the uniform load is kept on it. When looking at the equation for maximum deflection allowed, derived in the above section, the maximum uniform load allowed on this beam is 5,270 N/m as shown on the right:
To find an accurate result using the methods outlined above, it is necessary to iterate through these equations many times with increasing values of x. Using an algorithm makes this process more convenient and much faster than calculating by hand. The pseudo-code for these algorithms are as follows:
E = 2 ×10^11
I = 3.4 ×10^−4
w = 10,500
L = 14
step = 0.25
x1 = a vector starting from 0, increasing by the value of step, and ending when x1 = (14 + step) y1 = a vector of zeros, the same size as x1
index = 1
for each value of x,
use Euler’s equation and iterate through indices of y as x increases
index = index + 1
end
plot x2 and y1 in blue initialize x2 and y1 labels declare x2 and y1 limits
E = 2 ×10^11
I = 3.4 ×10^−4
w = 10,500
L = 14
x2 = a vector starting from 0, increasing by the value of step, and ending when x2 = (14 + step) y2 = analytical equation performed on each index of x2
hold on
plot x2 and y2 in red
create a legend for analytical and numerical results
Figure 2: Distance vs. deflection of a 14m beam using numerical and analytical methods
Based on Figure 2 above, we see that the numerical solution with a smaller step size (0.25m) is much more accurate compared to the numerical approach with step size = 1m. The numerical solution with a larger step size of 1m is less accurate, as it deviates further from the analytical solution at a quicker rate. This trend suggests that the greater the step size, the less accurate the results. This is logical because the accuracy and precision of each calculation will increase by performing more approximations to the actual solution. Both numerical approaches deviate further from the analytical approach as x approaches 14m. The maximum absolute value of the deflection is 0.0794m using numerical analysis when the step size is 0.25m and 0.0772m using analytical analysis. In each instance, the maximum absolute value of the deflection occurred when the length is 7m, or at half the length of the beam. Furthermore, it is interesting to see how must of the error begins to occur after the length exceeds 3 meters. This suggests that the beam may not be able to support the uniform load presented since the three curves deviate from each other as length increases.
Figure 3 below shows the relative error between numerical (Euler's) and analytical solutions for two step sizes (h = 0.25 m & h = 1.0 m). It is evident from this plot that when the step size is larger, the deflection value deviates from that of the analytical approach much quicker. This plot also suggests that the error of the deflection increases as the length of the beam increases, with the most error occurring towards the end of the beam and the least in the beginning. As mentioned in the section above there is a spot after 3 meters where the analytical and numerical curves diverge from each other due to the exceedance presented by the uniform load. This analysis can be clearly seen in Figure 3, where the error exponentially increases after the length of the beam is long enough to provide no support to the uniform load. This makes sense because after a length of 2 meters the uniform load is technically too heavy for the characteristics of the beam explaining why the error increases exponentially.
Figure 3: Relative error between numerical and analytical solutions for two step sizes
The relative error presented above was calculated by subtracting the numerical deflection from the analytical deflection and diving it by the analytical deflection as shown to the right:
Figure 4 below shows the true error between numerical and analytical solutions for the two step sizes, 0.25m and 1m. Again, we see a trend that error increases as the length increases and the error grows at a greater rate when the step size is larger.
Figure 4: The error between numerical and analytical solutions for two step sizes
The true error presented above was calculated by subtracting the numerical deflection from the analytical deflection as seen on the right in equation 10:
As seen above both methods provide an accurate version of calculations relative to a time step. In this case, both methods demonstrated the deflection of the beam was extremely big due to the oversized uniform load. In the future the beam can be designed with stronger physical characteristics or the uniform load has to be under the allowed stress in order to be safe and efficient. In addition, in order to increase accuracy and precision in further trials, the the numerical method would be the best approach since the time step can be as small as desired. In other words, the smaller the time step is the more accurate the deflection will be at that given length, reducing error. However, there will be a greater amount of data since the calculations must be performed until the end of the length. In further trials the most efficient way to reduce error is to use a uniform load that can be safely held by the presented beam. Once the physical components are compatible, error can be reduced by picking a small step size like 0.001m, in order to have extremely detailed data on how the beam will behave as its deflection is calculated along its length. Click the box below to access the MATLAB code for these solution steps. This code calculates the deflection using numerical and analytical approaches. The program also calculates error. Lastly, all the results including the error are plotted. Comments are provided on the code to clarify steps.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Katherine Reiss
% 02/04/21
% This program uses analytical and numerical methods with Euler's equation
% to find the deflection in a beam.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Declare Constants
E = 2*(10^11); % modulus of elasticity (Pa)
I = 3.4*(10^-4); % moment of inertia (m^4)
w = 10500; % uniform load (N/m)
L = 14; % length of beam (m)
%% Declare variables
step1 = 0.25; % choose one step size (m)
step2 = 1; % choose a different step size (m)
x1 = (0:step1:14+step1)'; % creates a vector, x1, storing values from 0 to 100+step with step size of value 'step1'
x2 = (0:step1:14+step1)'; % creates another vector, x2, identical to x (in order to differentiate between solution types)
x3 = (0:step2:14+step2)'; % creates a vector, x3, storing values from 0 to 100+step with step size of value 'step2'
x4 = (0:step2:14+step2)'; % creates a vector, x4, storing values from 0 to 100+step with step size of value 'step2'
size1 = length(x1); % size is a single value of the length of vector x1
size3 = length(x3); % size is a single value of the length of vector x3
y1 = zeros(size1,1); % initializes a zero matrix for the deflection values
y3 = zeros(size3,1); % initializes a zero matrix for the deflection values
index = 1; % tracks the index of the x1 vector
index2 = 1; % tracks the index of the x2 vector
%% Implement Numerical Solution
for x1 = 0:step1:14 % for loop iterates through x1 vector from 0 to 14 with a step size of value 'step1'
y1(index+1,1) = y1(index,1) + (((w.*L.*(x1.^2))./(4.*E.*I)) - ((w.*(x1.^3))./(6.*E.*I)) - ((w*(L^3))/(24*E*I)))*step1; % Euler's equation
index = index + 1; % increase the index by one for each iteration
end % for loop ends when x > 14
for x3 = 0:step2:14 % for loop iterates through x1 vector from 0 to 14 with a step size of value 'step2'
y3(index2+1,1) = y3(index2,1) + (((w.*L.*(x3.^2))./(4.*E.*I)) - ((w.*(x3.^3))./(6.*E.*I)) - ((w*(L^3))/(24*E*I)))*step2; % Euler's equation
index2 = index2 + 1; % increase the index by one for each iteration
end % for loop ends when x > 14
(1)
%% Implement Analytical Solution
y2 = ((w.*L.*(x2.^3))./(12.*E.*I)) - ((w.*(x2.^4))./(24.*E.*I)) - ((x2).*((w.*(L.^3))./(24.*E.*I))); % analytical solution
y4 = ((w.*L.*(x4.^3))./(12.*E.*I)) - ((w.*(x4.^4))./(24.*E.*I)) - ((x4).*((w.*(L.^3))./(24.*E.*I))); % analytical solution
(2)
%% Numerical Analysis Plot
plot(x2, y1, '-xb', 'linewidth', 1) % plot numerical solution in blue
xlabel('Distance (m)') % x label declaration
ylabel('Deflection (m)') % y label declaration
ylim([-0.09,0]) % y limits declaration
xlim([0, 14]) % x limits declaration
grid on % grid on plot is turned on
box on % box around plot is turned on
hold on % graphing solutions on the same plot: off
plot(x4, y3, '-xk', 'linewidth', 1) % plot analytical solution in black with step size 1
set (gca, 'fontsize', 20) % sets font size
%% Analytical Analysis Plot
hold on % graphing solutions on the same plot
plot(x2, y2, '-xr', 'linewidth', 1) % plot analytical solution in red
set (gca, 'fontsize', 20) % sets font size
legend('Numerical with step = 0.25m','Numerical with step = 1m', 'Analytical') % legend declaration
title('Distance vs. Deflection of a 14m Beam using Numerical and Analytical Methods'); % title declaration
%% Calculating Error
trueErrorPt25 = y2 - y1; % true error calculation for step size 0.25
relativeErrorPt25 = 100.*abs((y1-y2)./y2); % relative error calculation for step size 0.25
trueError1 = y4 - y3; % true error calculation for step size 1
relativeError1 = 100.*abs((y4-y3)./y4); % relative error calculation for step size 1
(3)
hold off % next plot will be in a different window
plot(x2, relativeErrorPt25, '-xr', 'linewidth', 1); % plot relative error for a step size of 0.25
hold on % plot new one on same plot
plot(x4, relativeError1, '-or', 'linewidth', 1); % plot relative error for a step size of 1
title('Relative Error between Numerical and Analytical Solutions for Two Step Sizes', 'fontsize', 12); % declare title
legend('Step = 0.25m', 'Step = 1m') % declare legend
xlabel('Distance (m)'); % declare x label
ylabel('Error (%)'); % declare y label
xlim([0,14]); % declare x limit
ylim([0,110]); % declare y limit
set (gca, 'fontsize', 20) % sets font size
box on % turn box on for plot
grid on % turn grid on for plot
hold off % next plot will be in a different window
plot(x2, trueErrorPt25, '-xr', 'linewidth', 1); % plot relative error for a step size of 0.25
hold on % plot new one on same plot
plot(x4, trueError1, '-or', 'linewidth', 1); % plot relative error for a step size of 1
title('True Error between Numerical and Analytical Solutions for Two Step Sizes', 'fontsize', 12); % declare title
legend('Step = 0.25m', 'Step = 1m') % declare legend
xlabel('Distance (m)'); % declare x label
ylabel('Error (m)'); % declare y label
xlim([0,14]); % declare x limit
ylim([0,0.02]); % declare y limit
set (gca, 'fontsize', 20) % sets font size
box on % turn box on for plot
grid on % turn grid on for plot
(4)
(1) The above code shows the setup of the program where the important beam parameters as seen in Table 1 are defined. In addition, the first and second numerical approaches are executed with a time steps (h) of 0.25m and 1.0m respectively. In this case, a vectors “size1” and "size3" are created to start from 0 increase by h (0.25 & 1.0) and stop at 14. The length of this vector is translated to the variables “y1” and "y3," which are then used to define the size of the deflection vectors. Variables "index" and "index2" are defined to act as a counter for the numerical method (counting from 0 – 14). A for loop is created to iterate 15 times starting at 0, increasing by h (0.25 & 1.0), and finishing at 14. The for loop will reiterate the numerical method 15 different times storing the deflection value according to its length. It is crucial to notice that the vector y1 will contain 60 deflection approximations while y3 will only contain 15. This is so because y1 contains a time step of 0.25 which means it will increase its calculating accuracy by a factor of 4, when compared to y3 which contains a step size of 1.
(2) In the second part of the code, the analytical solution is implemented for both step sizes as 'y2' calculates the analytical solution for the step size 0.25m and 'y4' calculates the analytical solution for the step size 1.0m. The deflection is then calculated utilizing vector multiplication, exponent raise, and division, while following the analytical equation defined above.
(3 & 4) The error between the analytical and the numerical methods are calculated as seen above. The true error is calculated for a 0.25m and a 1.0m step size and is assigned to the variables trueErrorPt25 and trueError1 respectively. As seen above true error is the difference between the analytical values and the numerical values. Furthermore, the relative error was calculated for both step sizes as seen in the variables relativeErrorPt25 and relativeError1, by taking the true error using vector division to divide by the analytical values and finally use vector multiplication to multiply by 100 and obtain a percentage.
In the last sections of the code the calculated numerical and analytical results are plotted. Plots are included for the deflection of both numerical and analytical methods, relative error, and true error. Each plot is added utilizing the plot(x, y, specs) function and specific graph characteristics are specified such as a title, axis labels, and legend. Some of the customization functions include: xlabel(), ylabel(), legend(), title(), xlim(), and ylim().
Figure 5. For a more in depth look at how beams deflect watch this video. Note this video explores several Analytical methods on how to calculate the deflection of a beam, including the analytical method seen above. YouTube
Throughout this process, the primary changes to the lab included adding more description to the written sections, providing more mathematical background, and increasing the size of plots and adding titles in order to ensure clarity. As outlined above, the comments from the peer review process and the TA comments were used to iterate the submission and provide a final lab that took all of these changes into consideration. From working together as a team, we were able to take bits and pieces of each other's initial submissions, showcasing both of our strengths while acknowledging and fixing the parts of the lab that were weaker. The work distribution as seen on the first page of this document was equal. Communication was present amongst our team, as we made ourselves available via facetime and phone calls. Overall, this experience was a positive one and our collaborative efforts allowed us to provide a great final product.
[1] EGR 211 Materials and Mechanics Class Notes by Professor Lutzweiler
[2] Set of Extensive Blogs, How To Structural Steel Beams Are Made And Are They Cost-Effective Or Not? https://setofextensiveblogs.wordpress.com/2016/05/25/how-to-structural-steel-beams-are-made-and-are-they-cost-effective-or-not/
[3] The Efficient Engineer, Understanding the Deflection of Beams, https://www.youtube.com/watch?v=MvBqCeZllpQ