FEM_A. J. M. _11.17

目录

文档中包含的MATLAB程序如下：

离散系统相关程序

  1. problem1.m：解决弹簧问题，计算系统刚度矩阵、位移和反力。
% MATLAB codes for Finite Element Analysis
% problem1.m % A.J.M. Ferreira, N. Fantuzzi 2019
%%
% clear memory
clear
elementNodes = [1 2;2 3;2 4]; % elementNodes: connections at elements
% numberElements: number of Elements
numberElements = size(elementNodes,1);
% numberNodes: number of nodes
numberNodes = 4;
% for structure:
% displacements: displacement vector
% stiffness: stiffness matrix
displacements = zeros(numberNodes,1);
force = zeros(numberNodes,1);
stiffness = zeros(numberNodes);
% applied load at node 2
force(2) = 10.0;
% computation of the system stiffness matrix
for e = 1:numberElements
    % elementDof: element degrees of freedom (Dof)
    elementDof = elementNodes(e,:);
    stiffness(elementDof,elementDof) =...
        stiffness(elementDof,elementDof) + [1 -1;-1 1];
end
% boundary conditions and solution
% prescribed dofs
prescribedDof = [1;3;4];
% free Dof: activeDof
activeDof = setdiff((1:numberNodes)',prescribedDof);
% solution
displacements(activeDof) =...
    stiffness(activeDof,activeDof)\force(activeDof);
% output displacements/reactions
outputDisplacementsReactions(displacements,stiffness,...
    numberNodes,prescribedDof)

杆或桁架相关程序

  1. problem2.m：解决等参杆在均布载荷下的问题，计算位移、应力等。
% MATLAB codes for Finite Element Analysis
% A.J.M. Ferreira, N. Fantuzzi 2019 % problem2.m
%% % clear memory
clear; close all
% E; modulus of elasticity
% A: area of cross section % L: length of bar
E = 30e6; A = 1; EA = E*A; L = 90; p = 50;
% generation of coordinates and connectivities
% numberElements: number of elements
numberElements = 3; % generation equal spaced coordinates
nodeCoordinates = linspace(0,L,numberElements+1);
xx = nodeCoordinates; % numberNodes: number of nodes
numberNodes = size(nodeCoordinates,2);
% elementNodes: connections at elements
ii = 1:numberElements;
elementNodes(:,1) = ii;
elementNodes(:,2) = ii+1;
% for structure:
% displacements: displacement vector
% force : force vector
% stiffness: stiffness matrix
displacements = zeros(numberNodes,1);
force = zeros(numberNodes,1);
stiffness = zeros(numberNodes,numberNodes);
% computation of the system stiffness matrix and force vector
for e = 1:numberElements % elementDof: element degrees of freedom (Dof)
    elementDof = elementNodes(e,:);
    nn = length(elementDof);
    length_element = nodeCoordinates(elementDof(2))...
        -nodeCoordinates(elementDof(1));
    detJacobian = length_element/2;
    invJacobian = 1/detJacobian;
    % central Gauss point (xi=0, weight W=2)
    [shape,naturalDerivatives] = shapeFunctionL2(0.0);
    Xderivatives = naturalDerivatives*invJacobian;
    % B matrix
    B = zeros(1,nn);
    B(1:nn) = Xderivatives(:);
    stiffness(elementDof,elementDof) = stiffness(elementDof,elementDof) + B'*B*2*detJacobian*EA;
    force(elementDof,1) =...
        force(elementDof,1) + 2*shape*p*detJacobian;
end
% prescribed dofs
prescribedDof = find(xx==min(nodeCoordinates(:))...
    | xx==max(nodeCoordinates(:)))';
% free Dof: activeDof
activeDof = setdiff((1:numberNodes)',prescribedDof);
% solution
GDof = numberNodes;
displacements = solution(GDof,prescribedDof,stiffness,force);
% output displacements/reactions
outputDisplacementsReactions(displacements,stiffness,...
    numberNodes,prescribedDof)
% stresses at elements
sigma = zeros(numberElements,1);
for e = 1:numberElements
    elementDof = elementNodes(e,:); % elementDof: element degrees of freedom (Dof)
    nn = length(elementDof);
    sigma(e) = E/length_element*([-1 1]*displacements(elementDof));
    length_element = nodeCoordinates(elementDof(2))...
        -nodeCoordinates(elementDof(1));
end
% drawing nodal displacements
figure; axes1 = axes;
hold on; box on; % displacements
figure
plot(axes1,nodeCoordinates,displacements,...
    'ok','markersize',8,'linewidth',1.5)
% graphical representation with interpolation for each element
interpNodes = 10;
for e = 1:numberElements
    nodeA = elementNodes(e,1);
    nodeB = elementNodes(e,2);
    XX = linspace(nodeCoordinates(nodeA),nodeCoordinates(nodeB),...
        interpNodes);
    ll = nodeCoordinates(nodeB)-nodeCoordinates(nodeA);
    % dimensionless coordinate
    xi = (XX - nodeCoordinates(nodeA))*2/ll - 1;
    phi1 = 0.5*(1 - xi);
    phi2 = 0.5*(1 + xi); % linear shape function
    u = phi1*displacements(nodeA) + phi2*displacements(nodeB); % displacement at the element
    plot(axes1,XX,u,'-k','linewidth',1.5)
    plot(axes1,XX,p*L*XX/2/EA.*(1 - XX/L), '--b','linewidth',1.5)
    % stress at the element
    plot(axes1,XX,sigma,'-k','linewidth',1.5)
    sigma = E/ll * ones(1,interpNodes) *...
        (displacements(nodeB) - displacements(nodeA));
    plot(axes1,XX,p*L/A*(0.5 - XX/L), '--b','linewidth',1.5)
end
set(axes1,'fontsize',18);
  1. problem3.m：解决固定杆与弹簧支撑问题（直接刚度法）。
% MATLAB codes for Finite Element Analysis
% problem3.m
% direct stiffness method % A.J.M. Ferreira, N. Fantuzzi 2019 % ref: D. Logan, A first course in the finite element method,
% third Edition, page 121, exercise P3-10
%% % clear memory
clear
% E; modulus of elasticity % A: area of cross section
% L: length of bar % k: spring stiffness
E = 70000; A = 200; k = 2000;
% generation of coordinates and connectivities
% numberElements: number of elements
numberElements = 3;
numberNodes = 4;
elementNodes = [1 2; 2 3; 3 4];
nodeCoordinates = [0 2000 4000 4000];
xx = nodeCoordinates;
% for structure:
% displacements: displacement vector
% force : force vector
% stiffness: stiffness matrix
displacements = zeros(numberNodes,1);
stiffness = zeros(numberNodes,numberNodes);
force = zeros(numberNodes,1);
% applied load at node 2
force(2) = 8000;
% computation of the system stiffness matrix
ea = zeros(1,numberElements);
for e = 1:numberElements
    L = nodeCoordinates(elementDof(2))...
        % elementDof: element degrees of freedom (Dof)
        elementDof = elementNodes(e,:);
        - nodeCoordinates(elementDof(1));
    if e < 3
        ea(e) = E*A/L;
    else
        ea(e) = k;
    end
    stiffness(elementDof,elementDof) =...
        stiffness(elementDof,elementDof) + ea(e)*[1 -1;-1 1];
end
% boundary conditions and solution
prescribedDof = [1;4];
% free Dof: activeDof
activeDof = setdiff((1:numberNodes)',prescribedDof);
% solution
displacements = solution(numberNodes,prescribedDof,stiffness,force);
% output displacements/reactions
outputDisplacementsReactions(displacements,stiffness,...
    numberNodes,prescribedDof)
  1. problem3a.m：解决固定杆与弹簧支撑问题（等参公式）。
% MATLAB codes for Finite Element Analysis % problem3a.m % ref: D. Logan, A first course in the finite element method,
% with isoparametric formulation % third Edition, page 121, exercise P3-10
% A.J.M. Ferreira, N. Fantuzzi 2019
%%
% clear memory
clear
% E; modulus of elasticity
% A: area of cross section
% L: length of bar
E = 70000; A = 200; EA = E*A; k = 2000;
% generation of coordinates and connectivities
numberElements = 3;
numberNodes = 4;
elementNodes = [1 2; 2 3; 3 4];
nodeCoordinates = [0 2000 4000 4000];
xx = nodeCoordinates;
% for structure: % displacements: displacement vector
% force : force vector
% stiffness: stiffness matrix
displacements = zeros(numberNodes,1);
stiffness = zeros(numberNodes,numberNodes);
force = zeros(numberNodes,1);
% applied load at node 2
force(2) = 8000.0;
% computation of the system stiffness matrix
ea = zeros(1,numberElements);
for e = 1:numberElements % elementDof: element degrees of freedom (Dof)
    elementDof = elementNodes(e,:);
    if e < 3 % bar elements
        nn = length(elementDof);
        length_element = nodeCoordinates(elementDof(2))...
            -nodeCoordinates(elementDof(1));
        detJacobian = length_element/2;
        invJacobian = 1/detJacobian;
        % central Gauss point (xi=0, weight W=2)
        [shape,naturalDerivatives] = shapeFunctionL2(0.0);
        Xderivatives = naturalDerivatives*invJacobian;
        % B matrix
        B = zeros(1,nn);
        B(1:nn) = Xderivatives(:);
        ea(e) = E*A;
        stiffness(elementDof,elementDof) =...
            stiffness(elementDof,elementDof)...
            + B'*B*2*detJacobian*ea(e);
    else % spring element
        stiffness(elementDof,elementDof) =...
            stiffness(elementDof,elementDof)...
            + k*[1 -1;-1 1];
    end
end
% boundary conditions and solution
prescribedDof = [1;4];
% solution
displacements = solution(numberNodes,prescribedDof,stiffness,force);
% output displacements/reactions
outputDisplacementsReactions(displacements,stiffness,...
    numberNodes,prescribedDof)
  1. problem3vib.m：解决杆的自由振动问题，考虑不同的质量矩阵计算方法。
% MATLAB codes for Finite Element Analysis
% problem3vib.m % ref: J.N. Reddy, An introduction to the Finite Element Method,
% third Edition, page 86, example 2.5.4
% A.J.M. Ferreira, N. Fantuzzi 2019
%% % clear memory
clear
% E; modulus of elasticity % A: area of cross section % L: length of bar % rho: density
E = 70000; A = 200; EA = E*A; k = EA/4000; rho = 1000;
% generation of coordinates and connectivities
numberElements = 3;
numberNodes = 4;
elementNodes = [1 2; 2 3; 3 4];
nodeCoordinates = [0 2000 4000 4000];
xx = nodeCoordinates;
% for structure:
% displacements: displacement vector
% force : force vector
% stiffness: stiffness

### 二维桁架相关程序
1. **problem4.m**：解决二维桁架的静力分析问题，计算位移、应力等。
```matlab
% MATLAB codes for Finite Element Analysis
% problem4.m % A.J.M. Ferreira, N. Fantuzzi 2019
%%
% clear memory
clear
% E; modulus of elasticity
% A: area of cross section
E = 30e6; A = 2; EA = E*A;
% generation of coordinates and connectivities
numberElements = 3;
numberNodes = 4;
elementNodes = [1 2;1 3;1 4];
nodeCoordinates = [0 0;0 120;120 120;120 0];
xx = nodeCoordinates(:,1);
yy = nodeCoordinates(:,2);
% for structure:
% displacements: displacement vector
% force : force vector
% GDof: total number of degrees of freedom
% stiffness: stiffness matrix
displacements = zeros(GDof,1);
force = zeros(GDof,1);
GDof = 2*numberNodes;
% applied load at node 2
force(2) = -10000.0;
% computation of the system stiffness matrix
[stiffness] =...
    formStiffness2Dtruss(GDof,numberElements,...
    elementNodes,numberNodes,nodeCoordinates,xx,yy,EA);
% boundary conditions and solution
prescribedDof = [3:8]';
% solution
displacements = solution(GDof,prescribedDof,stiffness,force);
% drawing displacements
us = 1:2:2*numberNodes-1;
vs = 2:2:2*numberNodes;
figure
L = xx(2)-xx(1);
XX = displacements(us);
YY = displacements(vs);
dispNorm = max(sqrt(XX.^2+YY.^2));
scaleFact = 15000*dispNorm;
hold on
drawingMesh(nodeCoordinates+scaleFact*[XX YY],elementNodes,...
    'L2','k.--');
drawingMesh(nodeCoordinates,elementNodes,'L2','k.-');
axis equal
set(gca,'fontsize',18)
% stresses at elements
stresses2Dtruss(numberElements,elementNodes,...
    xx,yy,displacements,E)
% output displacements/reactions
outputDisplacementsReactions(displacements,stiffness,...
    GDof,prescribedDof)
  1. problem5.m：解决另一个二维桁架的静力分析问题。
% MATLAB codes for Finite Element Analysis % problem5.m % A.J.M. Ferreira, N. Fantuzzi 2019
%%
% clear memory
clear
% E; modulus of elasticity
% A: area of cross section
E = 70000;
A = 300;
EA = E*A;
% generation of coordinates and connectivities
elementNodes = [1 2;1 3;2 3;2 4;1 4;3 4;3 6;4 5;4 6;3 5;5 6];
nodeCoordinates = [0 0;0 3000;3000 0;3000 3000;6000 0;6000 3000];
numberElements = size(elementNodes,1);
numberNodes = size(nodeCoordinates,1);
xx = nodeCoordinates(:,1);
yy = nodeCoordinates(:,2);
% for structure:
% displacements: displacement vector
% GDof = 2*numberNodes;
% stiffness: stiffness matrix
% force : force vector
U = zeros(GDof,1);
force = zeros(GDof,1);
% applied load at node 2
force(4) = -50000;
force(8) = -100000;
force(12) = -50000;
% computation of the system stiffness matrix
[stiffness] =...
    formStiffness2Dtruss(GDof,numberElements,...
    elementNodes,numberNodes,nodeCoordinates,xx,yy,EA);
% boundary conditions and solution
prescribedDof = [1 2 10]';
% solution
displacements = solution(GDof,prescribedDof,stiffness,force);
us = 1:2:2*numberNodes-1;
vs = 2:2:2*numberNodes;
% drawing displacements
figure
L = xx(2)-xx(1);
XX = displacements(us);
YY = displacements(vs);
dispNorm = max(sqrt(XX.^2+YY.^2));
scaleFact = 2*dispNorm;
hold on
drawingMesh(nodeCoordinates+scaleFact*[XX YY],...
    elementNodes,'L2','k.--');
drawingMesh(nodeCoordinates,elementNodes,'L2','k.-');
axis equal
set(gca,'fontsize',18)
% output displacements/reactions
outputDisplacementsReactions(displacements,stiffness,...
    GDof,prescribedDof)
% stresses at elements
stresses2Dtruss(numberElements,elementNodes,...
    xx,yy,displacements,E)
  1. problem6.m：解决包含弹簧的二维桁架问题。
% MATLAB codes for Finite Element Analysis
% problem6.m
% ref: D. Logan, A first course in the finite element method,
% third Edition, mixing trusses with springs % A.J.M. Ferreira, N. Fantuzzi 2019
%% % clear memory
clear
% E; modulus of elasticity
% A: area of cross section
E = 210000;
A = 500;
EA = E*A;
% generation of coordinates and connectivities
nodeCoordinates = [0 0;-5000*cos(pi/4) 5000*sin(pi/4); -10000 0];
elementNodes = [1 2;1 3;1 4];
numberElements = size(elementNodes,1);
numberNodes = size(nodeCoordinates,1)+1; % spring added
xx = nodeCoordinates(:,1);
yy = nodeCoordinates(:,2);
% for structure:
% displacements: displacement vector
% stiffness: stiffness matrix
% force : force vector
GDof = 2*numberNodes;
U = zeros(GDof,1);
stiffness = zeros(GDof);
force = zeros(GDof,1);
% applied load at node 2
force(2) = -25000;
% computation of the system stiffness matrix
[stiffness] =...
    formStiffness2Dtruss(GDof,numberElements-1,...
    elementNodes,numberNodes,nodeCoordinates,xx,yy,EA);
% spring stiffness in global Dof
stiffness([2 7],[2 7]) = stiffness([2 7],[2 7]) + 2000*[1 -1;-1 1];
% boundary conditions and solution
prescribedDof = (3:8)';
% solution
displacements = solution(GDof,prescribedDof,stiffness,force);
% output displacements/reactions
outputDisplacementsReactions(displacements,stiffness,...
    GDof,prescribedDof)
% stresses at elements
stresses2Dtruss(numberElements-1,elementNodes,...
    xx,yy,displacements,E)
  1. problem5vib.m：解决二维桁架的自由振动问题。
% MATLAB codes for Finite Element Analysis % problem5vib.m % A.J.M. Ferreira, N. Fantuzzi 2019
%%
% clear memory
clear
% E; modulus of elasticity
% A: area of cross section
E = 70000;
A = 300;
EA = E*A;
rho = 1000;
rhoA = rho*A;
% generation of coordinates and connectivities
elementNodes = [1 2;1 3;2 3;2 4;1 4;3 4;3 6;4 5;4 6;3 5;5 6];
nodeCoordinates = [0 0;0 3000;3000 0;3000 3000;6000 0;6000 3000];
numberElements = size(elementNodes,1);
numberNodes = size(nodeCoordinates,1);
xx = nodeCoordinates(:,1);
yy = nodeCoordinates(:,2);
% for structure:
% displacements: displacement vector
% GDof = 2*numberNodes;
% stiffness: stiffness matrix
% force : force vector
U = zeros(GDof,1);
GDof = 2*numberNodes;
% computation of the system stiffness matrix
[stiffness] =...
    formStiffness2Dtruss(GDof,numberElements,...
    elementNodes,numberNodes,nodeCoordinates,xx,yy,EA);
% computation of the system stiffness matrix
[mass] =...
    formMass2Dtruss(GDof,numberElements,...
    elementNodes,numberNodes,nodeCoordinates,xx,yy,rhoA);
% boundary conditions and solution
prescribedDof = [1 2 10]';
% free vibration problem
[modes,eigenvalues] = eigenvalue(GDof,prescribedDof,...
    stiffness,mass,0);
us = 1:2:2*numberNodes-1;
vs = 2:2:2*numberNodes;
modeNumber = 1;
% drawing displacements
figure
L = xx(2)-xx(1);
XX = modes(us,modeNumber);
YY = modes(vs,modeNumber);
dispNorm = max(sqrt(XX.^2+YY.^2));
scaleFact = 1e12*dispNorm;
hold on
drawingMesh(nodeCoordinates+scaleFact*[XX YY],...
    elementNodes,'L2','k.--');
drawingMesh(nodeCoordinates,elementNodes,'L2','k.-');
axis equal
set(gca,'fontsize',18)
omega = sqrt(eigenvalues)

三维桁架相关程序

  1. problem7.m：解决三维桁架的静力分析问题。
% MATLAB codes for Finite Element Analysis % problem7.m % ref: D. Logan, A first course in the finite element method,
% A.J.M. Ferreira, N. Fantuzzi 2019 % third Edition, A 3D truss example
8
%%
% clear memory
clear
% E; modulus of elasticity
% A: area of cross section
E = 1.2e6;
A = [0.302;0.729;0.187]; % area for various sections
% generation of coordinates and connectivities
nodeCoordinates = [72 0 0; 0 36 0; 0 36 72; 0 0 -48];
numberElements = size(elementNodes,1);
numberNodes = size(nodeCoordinates,1);
xx = nodeCoordinates(:,1);
yy = nodeCoordinates(:,2);
elementNodes = [1 2;1 3;1 4];
% for structure:
% displacements: displacement vector
% force : force vector
% stiffness: stiffness matrix
% GDof: global number of degrees of freedom
GDof = 3*numberNodes;
U = zeros(GDof,1);
force = zeros(GDof,1);
% applied load at node 2
force(3) = -1000;
% stiffness matrix
[stiffness] =...
    formStiffness3Dtruss(GDof,numberElements,...
    elementNodes,numberNodes,nodeCoordinates,E,A);
% boundary conditions and solution
prescribedDof = [2 4:12]';
% solution
displacements = solution(GDof,prescribedDof,stiffness,force);
% output displacements/reactions
outputDisplacementsReactions(displacements,stiffness,...
    GDof,prescribedDof)
% stresses at elements
stresses3Dtruss(numberElements,elementNodes,nodeCoordinates,...
    displacements,E)
  1. problem8.m：解决另一个三维桁架的静力分析问题。
% MATLAB codes for Finite Element Analysis
% problem8.m
% ref: D. Logan, A first course in the finite element method,
% third Edition, A second 3D truss example
% A.J.M. Ferreira, N. Fantuzzi 2019
%%
% clear memory
clear
% E; modulus of elasticity
A = [100 100 100 100]; % area for various sections % A: area of cross section
E = 210000;
% generation of coordinates and connectivities
nodeCoordinates = [4000 4000 3000; 0 4000

### 三维桁架相关程序（续）
1. **problem7vib.m**：解决三维桁架的自由振动问题。
```matlab
% MATLAB codes for Finite Element Analysis
% A 3D truss example in free vibrations % problem7vib.m
% A.J.M. Ferreira, N. Fantuzzi 2019
%%
% clear memory
clear
% E; modulus of elasticity % A: area of cross section
E = 1.2e6;
A = [0.302;0.729;0.187]; % area for various sections
rho = [1;1;1]; % density for various sections
% generation of coordinates and connectivities
nodeCoordinates = [72 0 0; 0 36 0; 0 36 72; 0 0 -48];
elementNodes = [1 2;1 3;1 4];
numberElements = size(elementNodes,1);
numberNodes = size(nodeCoordinates,1);
xx = nodeCoordinates(:,1);
yy = nodeCoordinates(:,2);
% for structure:
% displacements: displacement vector
% force : force vector
% stiffness: stiffness matrix
% GDof: global number of degrees of freedom
GDof = 3*numberNodes;
U = zeros(GDof,1);
force = zeros(GDof,1);
% stiffness matrix
[stiffness] =...
    formStiffness3Dtruss(GDof,numberElements,...
    elementNodes,numberNodes,nodeCoordinates,E,A);
% mass matrix
[mass] =...
    formMass3Dtruss(GDof,numberElements,...
    elementNodes,numberNodes,nodeCoordinates,rho,A);
% boundary conditions and solution
prescribedDof = [2 4:12]';
% free vibration problem
[modes,eigenvalues] = eigenvalue(GDof,prescribedDof,...
    stiffness,mass,0);
us = 1:3:3*numberNodes-2;
vs = 2:3:3*numberNodes-1;
ws = 3:3:3*numberNodes;
modeNumber = 1;
% drawing displacements
figure
L = xx(2)-xx(1);
XX = modes(us,modeNumber);
YY = modes(vs,modeNumber);
ZZ = modes(ws,modeNumber);
dispNorm = max(sqrt(XX.^2+YY.^2+ZZ.^2));
scaleFact = 1e3*dispNorm;
hold on
drawingMesh(nodeCoordinates+scaleFact*[XX YY ZZ],...
    elementNodes,'L3','k.--');
drawingMesh(nodeCoordinates,elementNodes,'L3','k.-');
axis equal
set(gca,'fontsize',18)
view(45,45)
omega = sqrt(eigenvalues)

伯努利梁相关程序

  1. problem9.m：解决伯努利梁在不同边界条件下的静力分析问题。
% MATLAB codes for Finite Element Analysis % problem9.m % A.J.M. Ferreira, N. Fantuzzi 2019
%%
% clear memory
clear
% E; modulus of elasticity
% I: second moment of area % L: length of bar
E = 1;
I = 1;
EI = E*I;
% generation of coordinates and connectivities
numberElements = 2;
nodeCoordinates = linspace(0,1,numberElements+1)';
L = max(nodeCoordinates);
numberNodes = size(nodeCoordinates,1);
xx = nodeCoordinates(:,1);
elementNodes = zeros(numberElements,2);
for i = 1:numberElements
    elementNodes(i,1) = i;
    elementNodes(i,2) = i + 1;
end
% distributed load
P = -1;
% for structure:
% displacements: displacement vector
% force : force vector
% GDof = 2*numberNodes;
% stiffness: stiffness matrix
% GDof: global number of degrees of freedom
[stiffness,force] =...
    formStiffnessBernoulliBeam(GDof,numberElements,...
    elementNodes,numberNodes,xx,EI,P);
% boundary conditions and solution
% clamped-clamped
% fixedNodeU =[1 2*numberNodes+1]';
% fixedNodeV =[2 2*numberNodes+2]';
% simply supported-simply supported
fixedNodeU =[1 2*numberNodes+1]';
fixedNodeV = []';
%fixedNodeU =[1]';
% fixedNodeV =[2]'; % clamped at x=0
prescribedDof = [fixedNodeU;fixedNodeV];
% solution
displacements = solution(GDof,prescribedDof,stiffness,force);
% output displacements/reactions
outputDisplacementsReactions(displacements,stiffness,...
    GDof,prescribedDof)
% reordering displacements and rotations
W = displacements(1:2:2*numberNodes);
R = displacements(2:2:2*numberNodes);
% drawing nodal displacements
figure
plot(nodeCoordinates,W,'ok','markersize',8,'linewidth',1.5)
set(gca,'fontsize',18)
% graphical representation with interpolation for each element
drawInterpolatedBeam
  1. problem9a.m：解决一端固定一端弹簧支撑的伯努利梁在均布载荷下的问题。
% MATLAB codes for Finite Element Analysis
% problem9a.m % A.J.M. Ferreira, N. Fantuzzi 2019
%%
% clear memory
clear
% E;

### 伯努利梁相关程序（续）
1. **problem9a.m（续）**
```matlab
% E; modulus of elasticity
% I: second moment of area % L: length of bar
E = 1e6;
L = 10;
t = L / 1000;
I = 1 * t ^ 3 / 12;
EI = E * I;
% generation of coordinates and connectivities
numberElements = 3;
nodeCoordinates = linspace(0,L,numberElements + 1)';
L = max(nodeCoordinates);
numberNodes = size(nodeCoordinates,1);
xx = nodeCoordinates(:,1);
elementNodes = zeros(numberElements,2);
for i = 1:numberElements
    elementNodes(i,1) = i;
    elementNodes(i,2) = i + 1;
end
% distributed force
P = -1000;
% for structure:
% displacements: displacement vector
% stiffness: stiffness matrix
% force : force vector
% GDof = 2*numberNodes;
% GDof: global number of degrees of freedom
stiffnessSpring = zeros(GDof + 1);
forceSpring = zeros(GDof + 1,1);
% stiffess matrix and force vector
[stiffness,force] =...
    formStiffnessBernoulliBeam(GDof,numberElements,...
    elementNodes,numberNodes,xx,EI,P);
% spring added
stiffnessSpring(1:GDof,1:GDof) = stiffness;
forceSpring(1:GDof) = force;
k = 10;
stiffnessSpring([GDof - 1 GDof + 1],[GDof - 1 GDof + 1]) =...
    stiffnessSpring([GDof - 1 GDof + 1],[GDof - 1 GDof + 1]) + [k -k;-k k];
% boundary conditions and solution
fixedNodeU = [1]';
fixedNodeV = [2]';
prescribedDof = [fixedNodeU;fixedNodeV;GDof + 1];
% solution
displacements = solution(GDof + 1,prescribedDof,...
    stiffnessSpring,forceSpring);
% displacements
disp('Displacements')
jj = 1:GDof + 1;
format
[jj' displacements]
% reordering displacements and rotations
W = displacements(1:2:2*numberNodes);
R = displacements(2:2:2*numberNodes);
% drawing nodal displacements
figure
plot(nodeCoordinates,W,'ok','markersize',8,'linewidth',1.5)
set(gca,'fontsize',18)
% graphical representation with interpolation for each element
drawInterpolatedBeam
% exact solution by Bathe (Solutions Manual of Procedures...)
load = [L * P / 3;L * P / 3;L * P / 6];
K = E * I / L ^ 3 * [189 -108 27;-108 135 -54;27 -54 27 + k * L ^ 3 / E / I];
X = K \ load
plot([0; 3.3333; 6.6667; 10.0000],[0; X],'-xb',...
    'markersize',8,'linewidth',1.5)
  1. problem9vib.m：解决伯努利梁的自由振动问题。
% MATLAB codes for Finite Element Analysis
% A.J.M. Ferreira, N. Fantuzzi 2019 % problem9vib.m
%% % clear memory
clear
% E; modulus of elasticity
% L: length of bar
E = 1;
I = 1;
EI = E * I;
rho = 1;
A = 2.3;
rhoA = rho * A;
% generation of coordinates and connectivities
numberElements = 64;
nodeCoordinates = linspace(0,1,numberElements + 1)';
L = max(nodeCoordinates);
numberNodes = size(nodeCoordinates,1);
xx = nodeCoordinates(:,1);
elementNodes = zeros(numberElements,2);
for i = 1:numberElements
    elementNodes(i,1) = i;
    elementNodes(i,2) = i + 1;


    ### 伯努利梁相关程序（续）
1. **problem9vib.m（续）**
```matlab
end
% for structure:
% displacements: displacement vector
% stiffness: stiffness matrix
% mass: mass matrix
% GDof: global number of degrees of freedom
GDof = 2*numberNodes;
% stiffess matrix
[stiffness,~] =...
    formStiffnessBernoulliBeam(GDof,numberElements,...
    elementNodes,numberNodes,xx,EI,1);
% stiffess matrix
[mass] =...
    formMassBernoulliBeam(GDof,numberElements,...
    elementNodes,numberNodes,xx,rhoA);
% boundary conditions and solution
% clamped-clamped
% fixedNodeU =[1 2*numberNodes+1]';
% simply supported-simply supported
% fixedNodeV =[2 2*numberNodes+2]';
% clamped at x=0
%fixedNodeU =[1]';
% fixedNodeV =[2]';
fixedNodeU =[1 2*numberNodes+1]';
fixedNodeV = []';
prescribedDof = [fixedNodeU;fixedNodeV];
% free vibration problem
[modes,eigenvalues] = eigenvalue(GDof,prescribedDof,...
    stiffness,mass,0);
% natural frequencies
omega = sqrt(eigenvalues);
% exact frequencies
omega_exact(1,1) = pi^2*sqrt(EI/rhoA/L^4);
omega_exact(2,1) = 4*pi^2*sqrt(EI/rhoA/L^4);
omega_exact(3,1) = 9*pi^2*sqrt(EI/rhoA/L^4);
  1. problem9buk.m：解决伯努利梁的稳定性分析问题。
% MATLAB codes for Finite Element Analysis % problem9buk.m % A.J.M. Ferreira, N. Fantuzzi 2019
%%
% clear memory
clear
% E; modulus of elasticity
% I: second moment of area % L: length of bar
E = 1;
I = 1;
EI = E * I;
% generation of coordinates and connectivities
numberElements = 64;
nodeCoordinates = linspace(0,1,numberElements + 1)';
L = max(nodeCoordinates);
numberNodes = size(nodeCoordinates,1);
xx = nodeCoordinates(:,1);
elementNodes = zeros(numberElements,2);
for i = 1:numberElements
    elementNodes(i,1) = i;
    elementNodes(i,2) = i + 1;
end
% for structure:
% displacements: displacement vector
% stiffness: stiffness matrix
% GDof = 2*numberNodes;
% stability: geometric matrix
% GDof: global number of degrees of freedom
[stiffness,~] =...
    formStiffnessBernoulliBeam(GDof,numberElements,...
    elementNodes,numberNodes,xx,EI,1);
% stability matrix
[stability] =...
    formStabilityBernoulliBeam(GDof,numberElements,...
    elementNodes,numberNodes,xx);
% boundary conditions and solution
% clamped-clamped
% fixedNodeU =[1 2*numberNodes+1]';
% fixedNodeV =[2 2*numberNodes+2]';
fixedNodeU =[1 2*numberNodes+1]';
fixedNodeV = []'; % simply supported-simply supported
% clamped at x=0
%fixedNodeU =[1]';
% fixedNodeV =[2]';
prescribedDof = [fixedNodeU;fixedNodeV];
% free vibration problem
[modes,eigenvalues] = eigenvalue(GDof,prescribedDof,...
    stiffness,stability,0);
% natural frequencies
N0 = eigenvalues;
% exact frequencies simply-supported beam
N0_exact(1,1) = pi^2*EI/L^2;
N0_exact(2,1) = 4*pi^2*EI/L^2;
N0_exact(3,1) = 9*pi^2*EI/L^2;

二维框架相关程序

  1. problem10.m：解决二维框架的静力分析问题。
% MATLAB codes for Finite Element Analysis
% A.J.M. Ferreira
posted @ 2024-11-17 03:08  redufa  阅读(56)  评论(0)    收藏  举报