有限元法中，求解2D梁问题

如图所示的平面框架，已知杨氏模量E=200GPa，梁横截面积A=1.0e-2㎡,横截面转动惯量I=1E-4m⁴确定节点位移和约束反力。

 

 

需要使用到主程序fem.m

及两个子程序：

Beam2D2Node_Stiffness.m

Beam2D2Node_Assemble.m

----------------------------------------------------------------------------------------------------------------

主程序：

fem.m

% %--------------------------------------------

E=200e9;
A=1e-2;
I=1e-4;
l=2;% length of beam
% elemental stiffness
k1=Beam2D2Node_Stiffness(E,I,A,l);
theta=-60/180*pi;
lambda=[cos(theta), sin(theta), 0 0 0 0;-sin(theta), cos(theta), 0 0 0 0; 0 0 1 0 0 0;0 0 0 cos(theta), sin(theta) 0; 0 0 0 -sin(theta), cos(theta) 0; 0 0 0 0 0 1];
k2=Beam2D2Node_Stiffness(E,I,A,l);
k2=lambda'*k2*lambda;

% Assemble global stiffness
K=zeros(9,9);
K=Beam2D2Node_Assemble(K,k1,1,2);
K=Beam2D2Node_Assemble(K,k2,2,3);

% converting forces to nodes
u=zeros(9,1);
b=zeros(9,1);

x=1;
xl=x/l;
N=[1-xl 0 0 xl 0 0;
0 1-3*xl^2+2*xl^3 x-2*x^2/l+x^3/l^2 0 3*xl^2-2*xl^3 -x^2/l+x^3/l^2;
0 -6*x/l^2+6*x^2/l^3 1-4*x/l+3*xl^2 0 6*x/l^2-6*x^2/l^3 -2*xl+3*xl^2];
N'*[0;-20e3;0];
b([2;3;5;6])=[-10e3;-5e3;-10e3;-5e3]

% boundry condition
Freenodes=[3,4,5,6,9];

% solving
u(Freenodes)=inv(K(Freenodes,Freenodes))*b(Freenodes)

----------------------------------------------------------------------------------------------------------------

子程序Beam2D2Node_Stiffness.m：

function k = Beam2D2Node_ElementStiffness(E,I,A,L)
%determines the local stiffness matrix
%Young's module E,Area A, moment of inertia I, length L of beam element
%Output: k(6×6)
k = [E*A/L,0,0,-E*A/L,0,0;

0,12*E*I/(L^3),6*E*I/(L^2),0,-12*E*I/(L^3),6*E*I/(L^2);

0,6*E*I/(L^2),4*E*I/L,0,-6*E*I/(L^2),2*E*I/L;

-E*A/L,0,0,E*A/L,0,0;

0,-12*E*I/(L^3),-6*E*I/(L^2),0,12*E*I/(L^3),-6*E*I/(L^2);

0,6*E*I/(L^2),2*E*I/L,0,-6*E*I/(L^2),4*E*I/L];

----------------------------------------------------------------------------------------------------------------

子程序Beam2D2Node_Assemble.m：

function y =Beam2D2Node_Assemble(KK,k,i,j)
%assemble stiffness matrix
%Input: global stiffness matrix KK, elemental stiffness matrix k, element nodes indexing in the global i,j
%Output: global stiffness matrix y
%-----------------------------------------
DOF(1)=3*i-2;
DOF(2)=3*i-1;
DOF(3)=3*i;
DOF(4)=3*j-2;
DOF(5)=3*j-1;
DOF(6)=3*j;
for n1=1:6
    for n2=1:6
        KK(DOF(n1),DOF(n2))=KK(DOF(n1),DOF(n2))+k(n1,n2);
    end
end
y = KK;

----------------------------------------------------------------------------------------------------------------