solving beam deflection curve

11 次查看（过去 30 天）
mo
mo 2023-12-5
回答： VBBV 2023-12-6
Ran in:
I had a question about plotting Beam deflection curve
The drawing was almost correct, but i had a percentage of error when measuring the deflection at 1500 and 3000 compared to the results I got from skyciv.
I do not understand why this break occurred in the third period, and despite my application of the boundary condition equation at 2000 that the equation value equal zero, but in the loop and on the curve the value is not zero at the support third support
clc,clear
E = 250e3; % Young's modulus
I = 2e7; % Moment of inertia
L = 1000; % Length of the beam
F = 5000; % Force at A
W1 = 200; % Distributed load 1
W2 = 20; % Distributed load 2
W3 = 0.1; % Distributed load 3
m = W3*(1000^2)/2
m = 50000
syms R1 R2 R3 c
eqf = R1 +R2+R3 - F -W1*L/3-W3*L-W2*L ==0
eqf = 
eqm = R2*L+2*R3*L-F*L/3-W1*L/3*0.5*L-W3*L*1.5*L+W3*L^2/2-2.5*W2*L^2==0
eqm = 
cdef = (R1*L^3)/6-(4/81)*F*L^3-W1*L^4/144+c*L == 0
cdef = 
def2 = R1*4*L^3/3+R2*L^3/6-(125/162)*F*L^3-3/16*W1*L^4-W3*L^4/48-W3*L^4/16+2*L*c==0
def2 = 
sol = solve(eqf,eqm,cdef,def2)
sol = struct with fields:
R1: 2183306111051851699891/67108864000000000 R2: 1119568478814815206359/33554432000000000 R3: 3471827368770369759157/134217728000000000 c: -1524645038459259098803/402653184000
R1 = double(sol.R1)
R1 = 3.2534e+04
R2 = double(sol.R2)
R2 = 3.3366e+04
R3 = double (sol.R3)
R3 = 2.5867e+04
c = double(sol.c)
c = -3.7865e+09
x = linspace(0, 3000, 300001);
y = zeros(size(x));
% Compute y values
for n = 1:length(x)
if x(n) < (1000/3)
y(n) = (1/(E*I)) * ((R1*x(n)^3)/6 + c*x(n));
elseif x(n) >= (1000/3) && x(n) < (2000/3)
y(n) = (1/(E*I)) * ((R1*x(n)^3)/6 - F*(x(n)-(L/3))^3/6 - W1*(x(n)-(L/3))^4/24 + c*x(n));
elseif x(n) >= (2000/3) && x(n) < 1000
y(n) = (1/(E*I)) * ((R1*x(n)^3)/6 - F*(x(n)-(L/3))^3/6 - W1*L/3*(x(n)-(L/2))^3/6 + c*x(n));
elseif x(n) >= 1000 && x(n) < 1500
y(n) = (1/(E*I)) * ((R1*x(n)^3)/6 - F*(x(n)-(L/3))^3/6 - W1*L/3*(x(n)-(L/2))^3/6 + R2*(x(n)-L)^3/6 + c*x(n));
elseif x(n) >= 1500 && x(n) < 2000
y(n) = (1/(E*I)) * ((R1*x(n)^3)/6 - F*(x(n)-(L/3))^3/6 - W1*L/3*(x(n)-(L/2))^3/6 + R2*(x(n)-L)^3/6 - (W3*((x(n)-3/2*L)^3)/6 + m*((x(n)-(3/2)*L)^2)/2) + c*x(n));
elseif x(n) >= 2000 && x(n) <= 3000
y(n) = (1/(E*I)) * ((R1*x(n)^3)/6 - F*(x(n)-(L/3))^3/6 - W1*L/3*(x(n)-(L/2))^3/6 + R2*(x(n)-L)^3/6 - (W3*((x(n)-3/2*L)^3)/6 + m*((x(n)-(3/2)*L)^2)/2) + R3*((x(n)-2*L)^3)/6 - W2*((x(n)-2*L)^4)/24 + c*x(n));
end
end
plot(x,y)
grid on
YC = (1/(E*I)) * ((R1*1500^3)/6 - F*(1500-(L/3))^3/6 - W1*L/3*(1500-(L/2))^3/6 + R2*(1500-L)^3/6 + c*1500)
YC = 0.1762
YE = (1/(E*I)) * ((R1*3000^3)/6 - F*(3000-(L/3))^3/6 - W1*L/3*(3000-(L/2))^3/6 + R2*(3000-L)^3/6 - (W3*((3000-3/2*L)^3)/6 + m*((3000-(3/2)*L)^2)/2) + R3*((3000-2*L)^3)/6 - W2*((3000-2*L)^4)/24 + c*3000)
YE = -1.2924
Y_C = (1/(E*I))*((R1*3*L^2)/2 -F*((3*L-(L/3))^2)/2 - W1*L/3*(3*L-(L/2)^2)/2+(R2*(3*L-L)^2)/2+c)
Y_C = 0.0204
Y_E = (1/(E*I))*((R1*3*L^2)/2 -F*((3*L-(L/3))^2)/2 - W1*L/3*(3*L-(L/2)^2)/2+(R2*(3*L-L)^2)/2-(W3*((3*L-3/2*L)^2)/2+ m*(3*L-(3/2)*L))+(R3*((3*L-2*L)^2)/2-W2*(((3*L-2*L)^3)/6))+c)
Y_E = 0.0223

请先登录，再回答此问题。

回答（1 个）

VBBV
VBBV 2023-12-6
Ran in:
if you use linspace function, it divides the length between initial and final points into non-equispaced values, however if you specify vector of values between 0 and 3000, then it does consider each value which includes points 1500 and 3000 of the beam.
x = 0:1:3000;
Ix = x(x==1500) % find the index of location 1500
Ix = 1500
x = linspace(0,3000,30000)
x = 1×30000
0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000 1.1000 1.2000 1.3000 1.4000 1.5001 1.6001 1.7001 1.8001 1.9001 2.0001 2.1001 2.2001 2.3001 2.4001 2.5001 2.6001 2.7001 2.8001 2.9001
Ix = x(x==1500) % no such point exists in the vector when you use linspace function
Ix = 1×0 empty double row vector

类别

在 Help Center 和 File Exchange 中查找有关 Mathematics 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by