Matlab trợ giúp tính toán Phân Tử Hữu Hạn

Matlab trợ giúp tính toán Phân Tử Hữu Hạn

Trước khi giới thiệu với các bác các hàm Matlab được viết để hỗ trợ FEA, tui xin giới thiệu lại phương pháp Phần Tử Hữu Hạn một chút.
FEA bao gồm 6 bước cơ bản.
1. Chia phần tử (Discretizing the domain).
2. Lập ma trận độ cứng và viết phương trình cân bằng cho mỗi phần tử (Writing the element stiffness matrix). MATLAB hỗ trợ bước này.
3. Lắp ghép ma trận toàn thể - ghép nối phần tử (Assembling the global stiffness matrix). MATLAB hỗ trợ bước này.
4. Áp dụng điều kiện đầu, điều kiện biên để giản ước hệ phương trình của các nút toàn thể. (Applying the boundary con***ions).
5. Giải phương trình các Nút. (Solving the equations) MATLAB hỗ trợ bước này.
6. Tính toán các thông số cần tìm. (Post-processing) MATLAB hỗ trợ bước này.

Phần Tử Lò Xo.
Hàm Matlab lập ma trận độ cứng (bước 2):
--------------------------------------------------------
function y = SpringElementStiffness(k)
%SpringElementStiffness This function returns the element stiffness
% matrix for a spring with stiffness k.
% The size of the element stiffness matrix
% is 2 x 2.
y = [k -k ; -k k];
--------------------------------------------------------

Hàm Matlab lắp ghép ma trận độ cứng toàn thể. (bước 3) Với hàm này khi sử dụng bạn cần khai báo trước ma trận K=zeros(n,n) với n là kích thước của ma trận toàn thể.
--------------------------------------------------------
function y = SpringAssemble(K,k,i,j)
%SpringAssemble This function assembles the element stiffness
% matrix k of the spring with nodes i and j into the
% global stiffness matrix K.
% This function returns the global stiffness matrix K
% after the element stiffness matrix k is assembled.
K(i,i) = K(i,i) + k(1,1);
K(i,j) = K(i,j) + k(1,2);
K(j,i) = K(j,i) + k(2,1);
K(j,j) = K(j,j) + k(2,2);
y = K;
----------------------------------------------------------

Bước 5 - Giải phương trình ma trận, trong Matlab đã có sẵn các công cụ để giải.
Bước 6 - Hàm Matlab tính lực tác dụng lên lò xo tại các nút.
----------------------------------------------------------
function y = SpringElementForces(k,u)
%SpringElementForces This function returns the element nodal force
% vector given the element stiffness matrix k
% and the element nodal displacement vector u.
y = k * u;
----------------------------------------------------------

Được Ben_Franklin sửa chữa / chuyển vào 23:24 ngày 06/07/2006

Được Ben_Franklin sửa chữa / chuyển vào 23:40 ngày 06/07/2006

Phần tử thanh tuyến tính - The Linear Bar Element
Hàm lập ma trận độ cứng -bước 2.
---------------------------------------------------
function y = LinearBarElementStiffness(E,A,L)
%LinearBarElementStiffness This function returns the element
% stiffness matrix for a linear bar with
% modulus of elasticity E, cross-sectional
% area A, and length L. The size of the
% element stiffness matrix is 2 x 2.
y = [E*A/L -E*A/L ; -E*A/L E*A/L];
----------------------------------------------------
Hàm ghép nối phần tử - bước 3
----------------------------------------------------
function y = LinearBarAssemble(K,k,i,j)
%LinearBarAssemble This function assembles the element stiffness
% matrix k of the linear bar with nodes i and j
% into the global stiffness matrix K.
% This function returns the global stiffness
% matrix K after the element stiffness matrix
% k is assembled.
K(i,i) = K(i,i) + k(1,1);
K(i,j) = K(i,j) + k(1,2);
K(j,i) = K(j,i) + k(2,1);
K(j,j) = K(j,j) + k(2,2);
y = K;
----------------------------------------------------
Tính lực tác dụng lên mỗi thanh tại mỗi nút - bước 6
----------------------------------------------------
function y = LinearBarElementForces(k,u)
%LinearBarElementForces This function returns the element nodal
% force vector given the element stiffness
% matrix k and the element nodal displacement
% vector u.
y = k * u;
-----------------------------------------------------
Tính ứng suất - bước 6
-----------------------------------------------------
function y = LinearBarElementStresses(k, u, A)
%LinearBarElementStresses This function returns the element nodal
% stress vector given the element stiffness
% matrix k, the element nodal displacement
% vector u, and the cross-sectional area A.
y = k * u/A;
-----------------------------------------------------
Được Ben_Franklin sửa chữa / chuyển vào 00:01 ngày 07/07/2006

Phần tử thanh toàn phương (3 nút) - The Quadratic Bar Element
Các bác sử dụng tương tự phần trước nhé.
---------------------------------------------------------------------------------------
function y = QuadraticBarElementStiffness(E,A,L)
%QuadraticBarElementStiffness This function returns the element
% stiffness matrix for a quadratic bar
% with modulus of elasticity E,
% cross-sectional area A, and length L.
% The size of the element stiffness
% matrix is 3 x 3.
y = E*A/(3*L)*[7 1 -8 ; 1 7 -8 ; -8 -8 16];
---------------------------------------------------------------------------------------
---------------------------------------------------------------------------------------
function y = QuadraticBarAssemble(K,k,i,j,m)
%QuadraticBarAssemble This function assembles the element stiffness
% matrix k of the quadratic bar with nodes i, j
% and m into the global stiffness matrix K.
% This function returns the global stiffness
% matrix K after the element stiffness matrix
% k is assembled.
K(i,i) = K(i,i) + k(1,1);
K(i,j) = K(i,j) + k(1,2);
K(i,m) = K(i,m) + k(1,3);
K(j,i) = K(j,i) + k(2,1);
K(j,j) = K(j,j) + k(2,2);
K(j,m) = K(j,m) + k(2,3);
K(m,i) = K(m,i) + k(3,1);
K(m,j) = K(m,j) + k(3,2);
K(m,m) = K(m,m) + k(3,3);
y = K;
---------------------------------------------------------------------------------------
---------------------------------------------------------------------------------------
function y = QuadraticBarElementForces(k,u)
%QuadraticBarElementForces This function returns the element nodal
% force vector given the element stiffness
% matrix k and the element nodal
% displacement vector u.
y = k * u;
---------------------------------------------------------------------------------------
---------------------------------------------------------------------------------------
function y = QuadraticBarElementStresses(k, u, A)
%QuadraticBarElementStresses This function returns the element
% nodal stress vector given the element
% stiffness matrix k, the element nodal
% displacement vector u, and the
% cross-sectional area A.
y = k * u/A;
---------------------------------------------------------------------------------------
Được Ben_Franklin sửa chữa / chuyển vào 00:06 ngày 07/07/2006

Dạng giàn khung - The Plane Truss Element.
Chúng ta phải quan tâm đến chiều dài L và góc theta của mỗi phần tử.
---------------------------------------------------------------------------------------
function y = PlaneTrussElementLength(x1,y1,x2,y2)
%PlaneTrussElementLength This function returns the length of the
% plane truss element whose first node has
% coordinates (x1,y1) and second node has
% coordinates (x2,y2).
y = sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1));
---------------------------------------------------------------------------------------
---------------------------------------------------------------------------------------
function y = PlaneTrussElementStiffness(E,A,L, theta)
%PlaneTrussElementStiffness This function returns the element
% stiffness matrix for a plane truss
% element with modulus of elasticity E,
% cross-sectional area A, length L, and
% angle theta (in degrees).
% The size of the element stiffness
% matrix is 4 x 4.
x = theta*pi/180;
C = cos(x);
S = sin(x);
y = E*A/L*[C*C C*S -C*C -C*S ; C*S S*S -C*S -S*S ;
-C*C -C*S C*C C*S ; -C*S -S*S C*S S*S];
---------------------------------------------------------------------------------------
---------------------------------------------------------------------------------------
function y = PlaneTrussAssemble(K,k,i,j)
%PlaneTrussAssemble This function assembles the element stiffness
% matrix k of the plane truss element with nodes
% i and j into the global stiffness matrix K.
% This function returns the global stiffness
% matrix K after the element stiffness matrix
% k is assembled.
K(2*i-1,2*i-1) = K(2*i-1,2*i-1) + k(1,1);
K(2*i-1,2*i) = K(2*i-1,2*i) + k(1,2);
K(2*i-1,2*j-1) = K(2*i-1,2*j-1) + k(1,3);
K(2*i-1,2*j) = K(2*i-1,2*j) + k(1,4);
K(2*i,2*i-1) = K(2*i,2*i-1) + k(2,1);
K(2*i,2*i) = K(2*i,2*i) + k(2,2);
K(2*i,2*j-1) = K(2*i,2*j-1) + k(2,3);
K(2*i,2*j) = K(2*i,2*j) + k(2,4);
K(2*j-1,2*i-1) = K(2*j-1,2*i-1) + k(3,1);
K(2*j-1,2*i) = K(2*j-1,2*i) + k(3,2);
K(2*j-1,2*j-1) = K(2*j-1,2*j-1) + k(3,3);
K(2*j-1,2*j) = K(2*j-1,2*j) + k(3,4);
K(2*j,2*i-1) = K(2*j,2*i-1) + k(4,1);
K(2*j,2*i) = K(2*j,2*i) + k(4,2);
K(2*j,2*j-1) = K(2*j,2*j-1) + k(4,3);
K(2*j,2*j) = K(2*j,2*j) + k(4,4);
y = K;
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
function y = PlaneTrussElementForce(E,A,L,theta,u)
%PlaneTrussElementForce This function returns the element force
% given the modulus of elasticity E, the
% cross-sectional area A, the length L,
% the angle theta (in degrees), and the
% element nodal displacement vector u.
x = theta * pi/180;
C = cos(x);
S = sin(x);
y = E*A/L*[-C -S C S]* u;
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
function y = PlaneTrussElementStress(E,L,theta,u)
%PlaneTrussElementStress This function returns the element stress
% given the modulus of elasticity E, the
% the length L, the angle theta (in
% degrees), and the element nodal
% displacement vector u.
x = theta * pi/180;
C = cos(x);
S = sin(x);
y = E/L*[-C -S C S]* u;
---------------------------------------------------------------------------------------
---------------------------------------------------------------------------------------
function y = PlaneTrussInclinedSupport(T,i,alpha)
%PlaneTrussInclinedSupport This function calculates the
% tranformation matrix T of the inclined
% support at node i with angle of
% inclination alpha (in degrees).
x = alpha*pi/180;
T(2*i-1,2*i-1) = cos(x);
T(2*i-1,2*i) = sin(x);
T(2*i,2*i-1) = -sin(x) ;
T(2*i,2*i) = cos(x);
y = T;
---------------------------------------------------------------------------------------

Giàn ko gian - The Space Truss Element
Sử dụng các hàm tương tự như trên.
---------------------------------------------------------------------------------------
function y = SpaceTrussElementLength(x1,y1,z1,x2,y2,z2)
%SpaceTrussElementLength This function returns the length of the
% space truss element whose first node has
% coordinates (x1,y1,z1) and second node has
% coordinates (x2,y2,z2).
y = sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) + (z2-z1)*(z2-z1));
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
function y = SpaceTrussElementStiffness(E,A,L,thetax,thetay,thetaz)
%SpaceTrussElementStiffness This function returns the element
% stiffness matrix for a space truss
% element with modulus of elasticity E,
% cross-sectional area A, length L, and
% angles thetax, thetay, thetaz
% (in degrees). The size of the element
% stiffness matrix is 6 x 6.
x = thetax*pi/180;
u = thetay*pi/180;
v = thetaz*pi/180;
Cx = cos(x);
Cy = cos(u);
Cz = cos(v);
w = [Cx*Cx Cx*Cy Cx*Cz ; Cy*Cx Cy*Cy Cy*Cz ; Cz*Cx Cz*Cy Cz*Cz];
y = E*A/L*[w -w ; -w w];
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
function y = SpaceTrussAssemble(K,k,i,j)
%SpaceTrussAssemble This function assembles the element stiffness
% matrix k of the space truss element with nodes
% i and j into the global stiffness matrix K.
% This function returns the global stiffness
% matrix K after the element stiffness matrix
% k is assembled.
K(3*i-2,3*i-2) = K(3*i-2,3*i-2) + k(1,1);
K(3*i-2,3*i-1) = K(3*i-2,3*i-1) + k(1,2);
K(3*i-2,3*i) = K(3*i-2,3*i) + k(1,3);
K(3*i-2,3*j-2) = K(3*i-2,3*j-2) + k(1,4);
K(3*i-2,3*j-1) = K(3*i-2,3*j-1) + k(1,5);
K(3*i-2,3*j) = K(3*i-2,3*j) + k(1,6);
K(3*i-1,3*i-2) = K(3*i-1,3*i-2) + k(2,1);
K(3*i-1,3*i-1) = K(3*i-1,3*i-1) + k(2,2);
K(3*i-1,3*i) = K(3*i-1,3*i) + k(2,3);
K(3*i-1,3*j-2) = K(3*i-1,3*j-2) + k(2,4);
K(3*i-1,3*j-1) = K(3*i-1,3*j-1) + k(2,5);
K(3*i-1,3*j) = K(3*i-1,3*j) + k(2,6);
K(3*i,3*i-2) = K(3*i,3*i-2) + k(3,1);
K(3*i,3*i-1) = K(3*i,3*i-1) + k(3,2);
K(3*i,3*i) = K(3*i,3*i) + k(3,3);
K(3*i,3*j-2) = K(3*i,3*j-2) + k(3,4);
K(3*i,3*j-1) = K(3*i,3*j-1) + k(3,5);
K(3*i,3*j) = K(3*i,3*j) + k(3,6);
K(3*j-2,3*i-2) = K(3*j-2,3*i-2) + k(4,1);
K(3*j-2,3*i-1) = K(3*j-2,3*i-1) + k(4,2);
K(3*j-2,3*i) = K(3*j-2,3*i) + k(4,3);
K(3*j-2,3*j-2) = K(3*j-2,3*j-2) + k(4,4);
K(3*j-2,3*j-1) = K(3*j-2,3*j-1) + k(4,5);
K(3*j-2,3*j) = K(3*j-2,3*j) + k(4,6);
K(3*j-1,3*i-2) = K(3*j-1,3*i-2) + k(5,1);
K(3*j-1,3*i-1) = K(3*j-1,3*i-1) + k(5,2);
K(3*j-1,3*i) = K(3*j-1,3*i) + k(5,3);
K(3*j-1,3*j-2) = K(3*j-1,3*j-2) + k(5,4);
K(3*j-1,3*j-1) = K(3*j-1,3*j-1) + k(5,5);
K(3*j-1,3*j) = K(3*j-1,3*j) + k(5,6);
K(3*j,3*i-2) = K(3*j,3*i-2) + k(6,1);
K(3*j,3*i-1) = K(3*j,3*i-1) + k(6,2);
K(3*j,3*i) = K(3*j,3*i) + k(6,3);
K(3*j,3*j-2) = K(3*j,3*j-2) + k(6,4);
K(3*j,3*j-1) = K(3*j,3*j-1) + k(6,5);
K(3*j,3*j) = K(3*j,3*j) + k(6,6);
y = K;
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
function y = SpaceTrussElementForce(E,A,L,thetax,thetay,thetaz,u)
%SpaceTrussElementForce This function returns the element force
% given the modulus of elasticity E, the
% cross-sectional area A, the length L,
% the angles thetax, thetay, thetaz
% (in degrees), and the element nodal
% displacement vector u.
x = thetax * pi/180;
w = thetay * pi/180;
v = thetaz * pi/180;
Cx = cos(x);
Cy = cos(w);
Cz = cos(v);
y = E*A/L*[-Cx -Cy -Cz Cx Cy Cz]*u;
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
function y = SpaceTrussElementStress(E,L,thetax,thetay,thetaz,u)
%SpaceTrussElementStress This function returns the element stress
% given the modulus of elasticity E, the
% length L, the angles thetax, thetay,
% thetaz (in degrees), and the element
% nodal displacement vector u.
x = thetax * pi/180;
w = thetay * pi/180;
v = thetaz * pi/180;
Cx = cos(x);
Cy = cos(w);
Cz = cos(v);
y = E/L*[-Cx -Cy -Cz Cx Cy Cz]*u;
---------------------------------------------------------------------------------------

Hay quá! Ben_Franklin đúng là một chuyên gia về món này đấy! Sau này mình sẽ học, và khi đấy sẽ nhờ cậu hướng dẫn cho nhé!
WJT.

Dầm - The Beam Element
---------------------------------------------------------------------------------------
function y = BeamElementStiffness(E,I,L)
%BeamElementStiffness This function returns the element
% stiffness matrix for a beam
% element with modulus of elasticity E,
% moment of inertia I, and length L.
% The size of the element stiffness
% matrix is 4 x 4.
y = E*I/(L*L*L)*[12 6*L -12 6*L ; 6*L 4*L*L -6*L 2*L*L ;
-12 -6*L 12 -6*L ; 6*L 2*L*L -6*L 4*L*L];
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
function y = BeamAssemble(K,k,i,j)
%BeamAssemble This function assembles the element stiffness
% matrix k of the beam element with nodes
% i and j into the global stiffness matrix K.
% This function returns the global stiffness
% matrix K after the element stiffness matrix
% k is assembled.
K(2*i-1,2*i-1) = K(2*i-1,2*i-1) + k(1,1);
K(2*i-1,2*i) = K(2*i-1,2*i) + k(1,2);
K(2*i-1,2*j-1) = K(2*i-1,2*j-1) + k(1,3);
K(2*i-1,2*j) = K(2*i-1,2*j) + k(1,4);
K(2*i,2*i-1) = K(2*i,2*i-1) + k(2,1);
K(2*i,2*i) = K(2*i,2*i) + k(2,2);
K(2*i,2*j-1) = K(2*i,2*j-1) + k(2,3);
K(2*i,2*j) = K(2*i,2*j) + k(2,4);
K(2*j-1,2*i-1) = K(2*j-1,2*i-1) + k(3,1);
K(2*j-1,2*i) = K(2*j-1,2*i) + k(3,2);
K(2*j-1,2*j-1) = K(2*j-1,2*j-1) + k(3,3);
K(2*j-1,2*j) = K(2*j-1,2*j) + k(3,4);
K(2*j,2*i-1) = K(2*j,2*i-1) + k(4,1);
K(2*j,2*i) = K(2*j,2*i) + k(4,2);
K(2*j,2*j-1) = K(2*j,2*j-1) + k(4,3);
K(2*j,2*j) = K(2*j,2*j) + k(4,4);
y = K;
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
function y = BeamElementForces(k,u)
%BeamElementForces This function returns the element nodal force
% vector given the element stiffness matrix k
% and the element nodal displacement vector u.
y = k * u;
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
function y = BeamElementShearDiagram(f, L)
%BeamElementShearDiagram This function plots the shear force
% diagram for the beam element with nodal
% force vector f and length L.
x = [0 ; L];
z = [f(1) ; -f(3)];
hold on;
title(''Shear Force Diagram'');
plot(x,z);
y1 = [0 ; 0];
plot(x,y1,''k'')
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
function y = BeamElementMomentDiagram(f, L)
%BeamElementMomentDiagram This function plots the bending moment
% diagram for the beam element with nodal
% force vector f and length L.
x = [0 ; L];
z = [-f(2) ; f(4)];
hold on;
title(''Bending Moment Diagram'');
plot(x,z);
y1 = [0 ; 0];
plot(x,y1,''k'')
---------------------------------------------------------------------------------------

Khung phẳng - The Plane Frame Element
---------------------------------------------------------------------------------------
function y = PlaneFrameElementLength(x1,y1,x2,y2)
%PlaneFrameElementLength This function returns the length of the
% plane frame element whose first node has
% coordinates (x1,y1) and second node has
% coordinates (x2,y2).
y = sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1));
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
function y = PlaneFrameElementStiffness(E,A,I,L,theta)
%PlaneFrameElementStiffness This function returns the element
% stiffness matrix for a plane frame
% element with modulus of elasticity E,
% cross-sectional area A, moment of
% inertia I, length L, and angle
% theta (in degrees).
% The size of the element stiffness
% matrix is 6 x 6.
x = theta*pi/180;
C = cos(x);
S = sin(x);
w1 = A*C*C + 12*I*S*S/(L*L);
w2 = A*S*S + 12*I*C*C/(L*L);
w3 = (A-12*I/(L*L))*C*S;
w4 = 6*I*S/L;
w5 = 6*I*C/L;
y = E/L*[w1 w3 -w4 -w1 -w3 -w4 ; w3 w2 w5 -w3 -w2 w5 ;
-w4 w5 4*I w4 -w5 2*I ; -w1 -w3 w4 w1 w3 w4 ;
-w3 -w2 -w5 w3 w2 -w5 ; -w4 w5 2*I w4 -w5 4*I];
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
function y = PlaneFrameAssemble(K,k,i,j)
%PlaneFrameAssemble This function assembles the element stiffness
% matrix k of the plane frame element with nodes
% i and j into the global stiffness matrix K.
% This function returns the global stiffness
% matrix K after the element stiffness matrix
% k is assembled.
K(3*i-2,3*i-2) = K(3*i-2,3*i-2) + k(1,1);
K(3*i-2,3*i-1) = K(3*i-2,3*i-1) + k(1,2);
K(3*i-2,3*i) = K(3*i-2,3*i) + k(1,3);
K(3*i-2,3*j-2) = K(3*i-2,3*j-2) + k(1,4);
K(3*i-2,3*j-1) = K(3*i-2,3*j-1) + k(1,5);
K(3*i-2,3*j) = K(3*i-2,3*j) + k(1,6);
K(3*i-1,3*i-2) = K(3*i-1,3*i-2) + k(2,1);
K(3*i-1,3*i-1) = K(3*i-1,3*i-1) + k(2,2);
K(3*i-1,3*i) = K(3*i-1,3*i) + k(2,3);
K(3*i-1,3*j-2) = K(3*i-1,3*j-2) + k(2,4);
K(3*i-1,3*j-1) = K(3*i-1,3*j-1) + k(2,5);
K(3*i-1,3*j) = K(3*i-1,3*j) + k(2,6);
K(3*i,3*i-2) = K(3*i,3*i-2) + k(3,1);
K(3*i,3*i-1) = K(3*i,3*i-1) + k(3,2);
K(3*i,3*i) = K(3*i,3*i) + k(3,3);
K(3*i,3*j-2) = K(3*i,3*j-2) + k(3,4);
K(3*i,3*j-1) = K(3*i,3*j-1) + k(3,5);
K(3*i,3*j) = K(3*i,3*j) + k(3,6);
K(3*j-2,3*i-2) = K(3*j-2,3*i-2) + k(4,1);
K(3*j-2,3*i-1) = K(3*j-2,3*i-1) + k(4,2);
K(3*j-2,3*i) = K(3*j-2,3*i) + k(4,3);
K(3*j-2,3*j-2) = K(3*j-2,3*j-2) + k(4,4);
K(3*j-2,3*j-1) = K(3*j-2,3*j-1) + k(4,5);
K(3*j-2,3*j) = K(3*j-2,3*j) + k(4,6);
K(3*j-1,3*i-2) = K(3*j-1,3*i-2) + k(5,1);
K(3*j-1,3*i-1) = K(3*j-1,3*i-1) + k(5,2);
K(3*j-1,3*i) = K(3*j-1,3*i) + k(5,3);
K(3*j-1,3*j-2) = K(3*j-1,3*j-2) + k(5,4);
K(3*j-1,3*j-1) = K(3*j-1,3*j-1) + k(5,5);
K(3*j-1,3*j) = K(3*j-1,3*j) + k(5,6);
K(3*j,3*i-2) = K(3*j,3*i-2) + k(6,1);
K(3*j,3*i-1) = K(3*j,3*i-1) + k(6,2);
K(3*j,3*i) = K(3*j,3*i) + k(6,3);
K(3*j,3*j-2) = K(3*j,3*j-2) + k(6,4);
K(3*j,3*j-1) = K(3*j,3*j-1) + k(6,5);
K(3*j,3*j) = K(3*j,3*j) + k(6,6);
y = K;
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
function y = PlaneFrameElementForces(E,A,I,L,theta,u)
%PlaneFrameElementForces This function returns the element force
% vector given the modulus of elasticity E,
% the cross-sectional area A, the moment of
% inertia I, the length L, the angle theta
% (in degrees), and the element nodal
% displacement vector u.
x = theta * pi/180;
C = cos(x);
S = sin(x);
w1 = E*A/L;
w2 = 12*E*I/(L*L*L);
w3 = 6*E*I/(L*L);
w4 = 4*E*I/L;
w5 = 2*E*I/L;
kprime = [w1 0 0 -w1 0 0 ; 0 w2 w3 0 -w2 w3 ;
0 w3 w4 0 -w3 w5 ; -w1 0 0 w1 0 0 ;
0 -w2 -w3 0 w2 -w3 ; 0 w3 w5 0 -w3 w4];
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];
y = kprime*T* u;
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
function y = PlaneFrameElementAxialDiagram(f, L)
%PlaneFrameElementAxialDiagram This function plots the axial force
% diagram for the plane frame element
% with nodal force vector f and length
% L.
x = [0 ; L];
z = [-f(1) ; f(4)];
hold on;
title(''Axial Force Diagram'');
plot(x,z);
y1 = [0 ; 0];
plot(x,y1,''k'')
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
function y = PlaneFrameElementShearDiagram(f, L)
%PlaneFrameElementShearDiagram This function plots the shear force
% diagram for the plane frame element
% with nodal force vector f and length
% L.
x = [0 ; L];
z = [f(2) ; -f(5)];
hold on;
title(''Shear Force Diagram'');
plot(x,z);
y1 = [0 ; 0];
plot(x,y1,''k'')
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
function y = PlaneFrameElementMomentDiagram(f, L)
%PlaneFrameElementMomentDiagram This function plots the bending
% moment diagram for the plane frame
% element with nodal force vector f
% and length L.
x = [0 ; L];
z = [-f(3) ; f(6)];
hold on;
title(''Bending Moment Diagram'');
plot(x,z);
y1 = [0 ; 0];
plot(x,y1,''k'')
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
function y = PlaneFrameInclinedSupport(T,i,alpha)
%PlaneFrameInclinedSupport This function calculates the
% tranformation matrix T of the inclined
% support at node i with angle of
% inclination alpha (in degrees).
x = alpha*pi/180;
T(3*i-2,3*i-2) = cos(x);
T(3*i-2,3*i-1) = sin(x);
T(3*i-2,3*i) = 0;
T(3*i-1,3*i-2) = -sin(x);
T(3*i-1,3*i-1) = cos(x);
T(3*i-1,3*i) = 0;
T(3*i,3*i-2) = 0;
T(3*i,3*i-1) = 0;
T(3*i,3*i) = 1;
y = T;

The Grid Element
---------------------------------------------------------------------------------------
function y = GridElementLength(x1,y1,x2,y2)
%GridElementLength This function returns the length of the
% grid element whose first node has
% coordinates (x1,y1) and second node has
% coordinates (x2,y2).
y = sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1));
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
function y = GridElementStiffness(E,G,I,J,L,theta)
%GridElementStiffness This function returns the element
% stiffness matrix for a grid
% element with modulus of elasticity E,
% shear modulus of elasticity G, moment of
% inertia I, torsional constant J, length L,
% and angle theta (in degrees).
% The size of the element stiffness
% matrix is 6 x 6.
x = theta*pi/180;
C = cos(x);
S = sin(x);
w1 = 12*E*I/(L*L*L);
w2 = 6*E*I/(L*L);
w3 = G*J/L;
w4 = 4*E*I/L;
w5 = 2*E*I/L;
kprime = [w1 0 w2 -w1 0 w2 ; 0 w3 0 0 -w3 0 ;
w2 0 w4 -w2 0 w5 ; -w1 0 -w2 w1 0 -w2 ;
0 -w3 0 0 w3 0 ; w2 0 w5 -w2 0 w4];
R = [1 0 0 0 0 0 ; 0 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];
y = R''*kprime*R;
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
function y = GridAssemble(K,k,i,j)
%GridAssemble This function assembles the element stiffness
% matrix k of the grid element with nodes
% i and j into the global stiffness matrix K.
% This function returns the global stiffness
% matrix K after the element stiffness matrix
% k is assembled.
K(3*i-2,3*i-2) = K(3*i-2,3*i-2) + k(1,1);
K(3*i-2,3*i-1) = K(3*i-2,3*i-1) + k(1,2);
K(3*i-2,3*i) = K(3*i-2,3*i) + k(1,3);
K(3*i-2,3*j-2) = K(3*i-2,3*j-2) + k(1,4);
K(3*i-2,3*j-1) = K(3*i-2,3*j-1) + k(1,5);
K(3*i-2,3*j) = K(3*i-2,3*j) + k(1,6);
K(3*i-1,3*i-2) = K(3*i-1,3*i-2) + k(2,1);
K(3*i-1,3*i-1) = K(3*i-1,3*i-1) + k(2,2);
K(3*i-1,3*i) = K(3*i-1,3*i) + k(2,3);
K(3*i-1,3*j-2) = K(3*i-1,3*j-2) + k(2,4);
K(3*i-1,3*j-1) = K(3*i-1,3*j-1) + k(2,5);
K(3*i-1,3*j) = K(3*i-1,3*j) + k(2,6);
K(3*i,3*i-2) = K(3*i,3*i-2) + k(3,1);
K(3*i,3*i-1) = K(3*i,3*i-1) + k(3,2);
K(3*i,3*i) = K(3*i,3*i) + k(3,3);
K(3*i,3*j-2) = K(3*i,3*j-2) + k(3,4);
K(3*i,3*j-1) = K(3*i,3*j-1) + k(3,5);
K(3*i,3*j) = K(3*i,3*j) + k(3,6);
K(3*j-2,3*i-2) = K(3*j-2,3*i-2) + k(4,1);
K(3*j-2,3*i-1) = K(3*j-2,3*i-1) + k(4,2);
K(3*j-2,3*i) = K(3*j-2,3*i) + k(4,3);
K(3*j-2,3*j-2) = K(3*j-2,3*j-2) + k(4,4);
K(3*j-2,3*j-1) = K(3*j-2,3*j-1) + k(4,5);
K(3*j-2,3*j) = K(3*j-2,3*j) + k(4,6);
K(3*j-1,3*i-2) = K(3*j-1,3*i-2) + k(5,1);
K(3*j-1,3*i-1) = K(3*j-1,3*i-1) + k(5,2);
K(3*j-1,3*i) = K(3*j-1,3*i) + k(5,3);
K(3*j-1,3*j-2) = K(3*j-1,3*j-2) + k(5,4);
K(3*j-1,3*j-1) = K(3*j-1,3*j-1) + k(5,5);
K(3*j-1,3*j) = K(3*j-1,3*j) + k(5,6);
K(3*j,3*i-2) = K(3*j,3*i-2) + k(6,1);
K(3*j,3*i-1) = K(3*j,3*i-1) + k(6,2);
K(3*j,3*i) = K(3*j,3*i) + k(6,3);
K(3*j,3*j-2) = K(3*j,3*j-2) + k(6,4);
K(3*j,3*j-1) = K(3*j,3*j-1) + k(6,5);
K(3*j,3*j) = K(3*j,3*j) + k(6,6);
y = K;
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
function y = GridElementForces(E,G,I,J,L,theta,u)
%GridElementForces This function returns the element force
% vector given the modulus of elasticity E,
% the shear modulus of elasticity G, the
% moment of inertia I, the torsional constant J,
% the length L, the angle theta (in degrees),
% and the element nodal displacement vector u.
x = theta*pi/180;
C = cos(x);
S = sin(x);
w1 = 12*E*I/(L*L*L);
w2 = 6*E*I/(L*L);
w3 = G*J/L;
w4 = 4*E*I/L;
w5 = 2*E*I/L;
kprime = [w1 0 w2 -w1 0 w2 ; 0 w3 0 0 -w3 0 ;
w2 0 w4 -w2 0 w5 ; -w1 0 -w2 w1 0 -w2 ;
0 -w3 0 0 w3 0 ; w2 0 w5 -w2 0 w4];
R = [1 0 0 0 0 0 ; 0 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];
y = kprime*R* u;