FEM一维问题MATLAB

求解任意数量elements的位移，应力和支反力

范例

代码

下面展示 代码片。

clc
clear
close all
for n = 2:100
% element length
L = 1/n*ones(1,n);

% element area
bottol = 25e-3;
top = 100e-3;
ele = (top-bottol)/n;
A = ones(1,n);
for i = 1:n
A(n+1-i) = ((bottol + (i-1)*ele)^2 +  (bottol + i*ele)^2 + (bottol + (i-1)*ele)*(bottol + i*ele))/3;
end

% elastic modulus
E = 200e9*ones(1,n);

F_glo = zeros(1,n+1);

% elemental stiffness matrix
for ii = 1 : n
    f = L(ii)*A(ii)*77e3/2*[1 1];
    k(ii,1:2,1:2) = E(ii)*A(ii)/L(ii)*[1 -1;-1 1];
    F_glo(ii) = F_glo(ii) + f(1);
    F_glo(ii+1) = F_glo(ii+1) + f(2);
end

% global stiffness matrix setup
K_glo = zeros(n+1,n+1);
% elemental stiffness matrix setup
k_loc = zeros(2,2);

% stiffness matrix assembly
for ii =1 : n
    index = [(ii-1)+1 (ii-1)+2];
    k_loc(:,:) = k(ii,:,:);
    K_glo(index,index) = K_glo(index,index) + k_loc;
end

% reduced stiffness matrix
K_r = K_glo(2:n+1,2:n+1);
% reduced force vector
F_r = F_glo(2:n+1)';
% reduced node displacement
Q_r = K_r\F_r;
% node displacement
Q = [0; Q_r];

% stress
for ii = 1 : 2
str(ii,:) = -[-1 1]*[Q(ii+1);Q(ii)]/L(ii)*E(ii);
end

% Reaction force
R = K_glo*Q - F_glo;

end