Civil engineering coding problem ( I dont know where the problem in my code is and how to solve it, thanks in advance )

Question

Mahmoud Josefsson 2024-4-16，10:39

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2107431-civil-engineering-coding-problem-i-dont-know-where-the-problem-in-my-code-is-and-how-to-solve-it

移动：Torsten 2024-4-16，11:46

clear, clc, close all
% Given parameters
a    = 4;                               % width (m)
b    = 3;                               % height (m)
E    = 210E9;                           % E-modulus (Pa)
I    = 0.5E-3;                          % Area moment of inertia (m^4)
R    = 1000E3;                          % Point load (N)
% Element connectivity
edof = [1 2 3 4; 3 4 5 6 ];             % Corrected [EL1; EL2]
% Node coordinates
coord = [0 0; a 0; a b];                % Corrected coordinates
% Boundary conditions
fdof = [1 2 3 4];                       % Free dofs
fixed_dofs = [5 6];                     % Fixed dofs
% Initialize global stiffness matrix
K = zeros(6,6);
% Global load vector
F = zeros(6,1);
F(4) = -R;
% Assemble global stiffness matrix
for i = 1:size(edof,1)
    n1 = edof(i,1)
    n2 = edof(i,2)
    x1 = coord(n1,1)
    y1 = coord(n1,2)
    x2 = coord(n2,1)
    y2 = coord(n2,2)
    L = sqrt((x2-x1)^2 + (y2-y1)^2);
    c = (x2-x1)/L;
    s = (y2-y1)/L;
    ke = beam2d_stiff(E,I,L,c,s);
    dofs = [2*n1-1 2*n1 2*n2-1 2*n2];
    K(dofs,dofs) = K(dofs,dofs) + ke;
end
n1 = 1
n2 = 2
x1 = 0
y1 = 0
x2 = 4
y2 = 0
n1 = 3
n2 = 4
x1 = 4
y1 = 3
Index in position 1 exceeds array bounds. Index must not exceed 3.
% Apply boundary conditions
Kf = K(fdof,fdof);
Ff = F(fdof) - K(fdof,fixed_dofs)*zeros(length(fixed_dofs),1);
% Solve for displacements
Df = Kf\Ff;
% Compute reaction forces
Rf = K(fdof,fixed_dofs)*zeros(length(fixed_dofs),1) - F(fdof);
% Compute axial forces
N = zeros(2,1);
for i = 1:size(edof,1)
    n1 = edof(i,1);
    n2 = edof(i,2);
    x1 = coord(n1,1);
    y1 = coord(n1,2);
    x2 = coord(n2,1);
    y2 = coord(n2,2);
    L = sqrt((x2-x1)^2 + (y2-y1)^2);
    c = (x2-x1)/L;
    s = (y2-y1)/L;
    ue = [Df(2*n1-1); Df(2*n1); Df(2*n2-1); Df(2*n2)]; 
    [N(i),~,~] = beam2d_forces(E,I,L,c,s,ue);
end
% Display results
disp('          Node     dx (mm)      dy (mm)')
disp('         ---------------------------------')
disp([(1:3)' Df(1:2:end)*1E3 Df(2:2:end)*1E3])
disp(' ');
disp('          Node     Rx (kN)      Ry (kN)')
disp('         ---------------------------------')
disp([(1:3)' Rf(1:2:end)/1E3 Rf(2:2:end)/1E3])
disp(' ');
disp('           EL        N (kN)');
disp('         -------------------');
disp([(1:2)' N/1E3]);
% Plot geometry and displacements
fac = 500;                              % scale factor, deformed shape
figure(1), clf
x  = [0 a a+b];   y = [0 b 0];          % coordinates, undeformed geometry
dx = Df(1:2:end); dy = Df(2:2:end);     % nodal displacements
plot(x,y,'ko-','LineWidth',1.5), hold all, axis equal
plot(x+fac*dx', y+fac*dy','Color',[0.4 0.4 0.4])
title('\rm 2D bar model')
legend('undeformed','deformed')
% Beam stiffness matrix function
function ke = beam2d_stiff(E,I,L,c,s)
ke = E*I/L^3 * [ c^2*L^2   c*s*L^2  -c^2*L^2  -c*s*L^2;
    c*s*L^2   s^2*L^2  -c*s*L^2  -s^2*L^2;
    -c^2*L^2  -c*s*L^2   c^2*L^2   c*s*L^2;
    -c*s*L^2  -s^2*L^2   c*s*L^2   s^2*L^2 ];
end
% Beam forces calculation function
function [N,V,M] = beam2d_forces(E,I,L,c,s,ue)
N = E*I/L^2 * [-c -s  c  s] * ue;
V = E*I/L^2 * [-s  c -s  c] * ue;
M = E*I/L   * [ s -c -s  c] * ue;
end