Question

Can you provide a sample programming for 2D plane frame in MATLAB ?

Answer 1

clear all

clc

format short

E=210e6;

I=9*10^-5;

A=1*10^-2;

L1=[sqrt(13),5,sqrt(13)];

theta=[atand(3/2) 0 360-atand(3/2)];

ndcon=[1 2;2 3;3 4];

rsdof=[1 2 3 10 11 12];

nelem = size(ndcon,1);

nodes = 4;

ndofn = 3;

nnode = 2;

tdofs = 12;

K=zeros(tdofs);

gstif=zeros(12,12);

for ielem = 1:nelem

    gbdof=[];

    a=A*E/L1(ielem);

    b=12*E*I/(L1(ielem))^3;

    c=b*L1(ielem)/2;

    d=(b*L1(ielem)^2)/3;

    e=d/2;

    C=cosd(theta(ielem));

    S=sind(theta(ielem));

    estif=[a 0 0 -a 0 0;

           0 b c 0 -b c;

           0 c d 0 -c e;

          -a 0 0 a 0 0;

           0 -b -c 0 b -c;

           0 c e 0 -c d];

       T=[C S 0 0 0 0;

           -S C 0 0 0 0;

           0 0 1 0 0 0;

           0 0 0 C S 0;

           0 0 0 -S C 0;

           0 0 0 0 0 1];

    elcon=ndcon(ielem,:);

for inode = 1:nnode

for idofn = 1:ndofn

            eldof = (elcon(inode)-1)*ndofn+idofn;

            gbdof= [gbdof, eldof];

end

end

    gstif(gbdof,gbdof) = gstif(gbdof,gbdof)+T'*estif*T;

end

gstif

w=5;

Ft=[20 (-w*5/2) -(w*25/12) 0 (-w*5/2) (w*25/12)];

force=[0 0 0 20 (-w*5/2) -(w*25/12) 0 (-w*5/2) (w*25/12) 0 0 0]';

activ = setdiff(1:tdofs,rsdof);

u=gstif(activ,activ)\force(activ);

U=[0;0;0;u;0;0;0]

F=gstif*U-[0 0 0 Ft 0 0 0]';

f=zeros(6,3);

for ielem = 1:nelem

    gbdof=[]; %global DOF

    a=A*E/L1(ielem);

    b=12*E*I/(L1(ielem))^3;

    c=b*L1(ielem)/2;

    d=(b*L1(ielem)^2)/3;

    e=d/2;

    C=cosd(theta(ielem));

    S=sind(theta(ielem));

    elcon=ndcon(ielem,:);

    estif=[a 0 0 -a 0 0;

           0 b c 0 -b c;

           0 c d 0 -c e;

          -a 0 0 a 0 0;

           0 -b -c 0 b -c;

           0 c e 0 -c d];

       T=[C S 0 0 0 0;

           -S C 0 0 0 0;

           0 0 1 0 0 0;

           0 0 0 C S 0;

           0 0 0 -S C 0;

           0 0 0 0 0 1];

for inode = 1:nnode

for idofn = 1:ndofn

            eldof = (elcon(inode)-1)*ndofn+idofn;

            gbdof= [gbdof, eldof];

end

end

       f(:,ielem)=estif*T*U(gbdof);

end

for ielem=1:nelem

x=[0,L1(ielem)];

y1=[-f(3,ielem),f(6,ielem)];

y2=[f(2,ielem),-f(5,ielem)];

y3=[-f(1,ielem),f(4,ielem)];

figure

plot(x,y1,x,y2,x,y3)

legend('BM','SF','AF')

title('member')

end

%OUTPUT

% gstif =

%

%    1.0e+05 *

%

%     1.8256    2.6658   -0.0726   -1.8256   -2.6658   -0.0726         0         0         0         0         0         0

%     2.6658    4.0471    0.0484   -2.6658   -4.0471    0.0484         0         0         0         0         0         0

%    -0.0726    0.0484    0.2097    0.0726   -0.0484    0.1048         0         0         0         0         0         0

%    -1.8256   -2.6658    0.0726    6.0256    2.6658    0.0726   -4.2000         0         0         0         0         0

%    -2.6658   -4.0471   -0.0484    2.6658    4.0653   -0.0030         0   -0.0181    0.0454         0         0         0

%    -0.0726    0.0484    0.1048    0.0726   -0.0030    0.3609         0   -0.0454    0.0756         0         0         0

%          0         0         0   -4.2000         0         0    6.0256   -2.6658    0.0726   -1.8256    2.6658   0.0726

%          0         0         0         0   -0.0181   -0.0454   -2.6658    4.0653    0.0030    2.6658   -4.0471    0.0484

%          0         0         0         0    0.0454    0.0756    0.0726    0.0030    0.3609   -0.0726   -0.0484    0.1048

%          0         0         0         0         0         0   -1.8256    2.6658   -0.0726    1.8256   -2.6658   -0.0726

%          0         0         0         0         0         0    2.6658   -4.0471   -0.0484   -2.6658    4.0471   -0.0484

%          0         0         0         0         0         0    0.0726    0.0484    0.1048   -0.0726   -0.0484    0.2097

%

% U =

%          0

%          0

%          0

%     0.0013

%    -0.0009

%    -0.0005

%     0.0012

%     0.0008

%     0.0003

%          0

%          0

%          0