Maximum Deflection for simply supported beam under arbitrary triangular loading

Question

I am trying to get the:

• maximum deflection
• location of maximum deflection
• deflection at any point

For a beam under an arbitrary triangular load (i.e. a, b, c, and w can vary). I tried to solve this via double integration by splitting it up into 3 piecewise functions but get stuck with the constants where I can't seem to get enough boundary conditions. I know that deflection at x = 0 is 0, deflection at x = L is 0, and the slope and deflection must be compatible at x = a and x = (a+b) at the cut points for the piecewise functions. (Where x = 0 at the left end of the beam, and positive going to the right).

Is there a smarter way to do this? I'm hoping to get solutions in terms of formulas, where I can just solve this once and plug and play in the future for a small coding project I'm doing.

Answer 1

Have a look at the last response in this link: Beam deformation under trapezoidal load

It uses singularity functions.

Answer 2

If you are doing this for a coding product, a good possibility is doing it by numerically integrating the applied load $w$. The advantage of this is that it can be easily extended to general forms of $w$ for which closed expressions may be hard to get, the disadvantage is that it's more computationally expensive than evaluating a known formula.

The relevant equations are as follows, for the simply supported case shown in the image.

Reactions $$ \begin{aligned} P_2&=\frac{1}{L}\int_0^Lw(x)x\,dx \\ P_1&=\int_0^Lw(x)\,dx -P_2 \end{aligned} $$

Internal shear and bending moment $$ \begin{aligned} V(x)&=\int_0^xw(x)\,dx-P_1 \\ M(x)&=-\int_0^xV(x)\,dx \end{aligned} $$

And, assuming Euler-Bernoulli beam theory, $$ \begin{aligned} \phi(x)&=\frac{M(x)}{EI}\\ \theta(x)&=\int\phi(x)\,dx +C \\ \delta(x)&=\int\theta(x)\,dx +Cx+D \end{aligned} $$ For the simply supported case $D=0$ and $C=-\left[\int\theta(x)\,dx\right](L)$

A MATLAB implementation of this is shown below

clear variables
close all
clc

%Inputs
%--Beam properties--
L=10;       %Length, meters
E=200;      %Modulus of elasticity, GPa
I=1;        %Moment of inertia about axis considered, m^4

%--Load properties--
a=3;        %Start of triangular load, m
b=4;        %Length of triangular load,m
w=1;        %Maximum magnitude, kN/m

%Boundary conditions (Only isostatic cases considered in this code)
fixed_slope=[0 0];          %0 indicates free dof, 1 means fixed dof.
fixed_deflection=[1 1];


%Simulation parameters
delta=0.01; %Max element length, m

%Calculations

% ADD A BIT HERE CHECKING THAT fixed_slope+fixed_deflection = 2 AND
% fixed_deflection has at leas one element or the problem is ill-posed
% ADD A BIT HERE CHECKING THAT a+b<=L and other validations
%

n=ceil(L/delta);
x=linspace(0,L,n)';
delta=x(2)-x(1);        %Recalculated because original delta may have changed

q=(x>a).*(x<=a+b).*(w/b*(x-a)); %Conditional definition of applied load

% Equivalently using for loops and conditionals
% q=zeros(length(x),1);
% for i=1:length(x)
%    if x(i)>=a && x(i)<=a+b
%        q(i)=w/b*(x(i)-a);
%    end
% end

figure
plot(x,q); xlabel('x [m]'); ylabel('Applied load [kN/m]')

%Determine reactions
qV=delta*trapz(q);      %Total force applied by w
qM=delta*trapz(q.*x);   %Moment generated by w around the leftmost end

P1=0;   %Leftmost reaction, kN
P2=0;   %Rightmost reaction, kN
M1=0;   %Leftmost moment, kN-m
M2=0;   %Rightmost moment, kN-m

% Five possibilities for reactions
if fixed_slope(1)==1 && fixed_deflection(1)==1
    P1=qV;  M1=qM;
end
if fixed_slope(1)==1 && fixed_deflection(2)==1
    P2=qV;  M1=qM-qV*L;
end
if fixed_slope(2)==1 && fixed_deflection(1)==1
    P1=qV;  M2=qM;
end
if fixed_slope(2)==1 && fixed_deflection(2)==1
    P2=qV;  M2=qM-qV/L;
end
if fixed_deflection(1)==1 && fixed_deflection(2)==1
    P2=qM/L;    P1=qV-P2;
end

%Calculate internal shear and bending moment

V=zeros(length(x),1);   %Shear, kN
M=zeros(length(x),1);   %Bending moment, kN-m

for i=1:length(x)
    V(i)=delta*trapz(q(1:i))-P1;
end

figure
plot(x,V); xlabel('x [m]'); ylabel('Shear [kN]')
for i=1:length(x)
    M(i)=-delta*trapz(V(1:i))-M1;
end

figure
plot(x,-M)      %I plot moments on tensioned fibre side, just a convention
xlabel('x [m]'); ylabel('Bending moment [kN-m]')

%Some unit transformation, due to the way inputs were stated.
E=E*1e6;        %kPa, or kN/m2

%Euler-Bernoulli beam theory calculation
phi=M./(E.*I);

theta=zeros(length(phi),1);
def=zeros(length(phi),1);

%Integrate radius of curvature to get slope
for i=1:length(theta)
    theta(i)=delta*trapz(phi(1:i));
end

%Integrate slope to get deflection
for i=1:length(def)
    def(i)=delta*trapz(theta(1:i));
end

% Five cases for boundary conditions
if fixed_slope(1)==1 && fixed_deflection(1)==1
    C=-theta(1); D=-def(1);
end
if fixed_slope(1)==1 && fixed_deflection(2)==1
    C=-theta(1); D=-def(end)-C*L;
end
if fixed_slope(2)==1 && fixed_deflection(1)==1
    C=-theta(end); D=-def(1);
end
if fixed_slope(2)==1 && fixed_deflection(2)==1
    C=-theta(end);  D=-def(end)-C*L;
end
if fixed_deflection(1)==1 && fixed_deflection(2)==1
    D=-def(1);  C=-def(end)/L;
end

%Recall that we integrated twice, so we have a constant of integration and 
%a constant of integration multiplied by x.
def=def+C*x+D;

figure
plot(x,def); xlabel('x [m]'); ylabel('deflection [m]')
%Plot a black line showing the original shape
hold on
plot([0 L],[0 0],'k','linewidth',2)

This produces the following results: Applied load

Internal shear

Bending moment

Deflection

If you have any questions, just leave a comment and I'll get to it.