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

The Space Frame Element
---------------------------------------------------------------------------------------
function y = SpaceFrameElementLength(x1,y1,z1,x2,y2,z2)
%SpaceFrameElementLength This function returns the length of the
% space frame 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 = SpaceFrameElementStiffness(E,G,A,Iy,Iz,J,x1,y1,z1,x2,y2,z2)
%SpaceFrameElementStiffness This function returns the element
% stiffness matrix for a space frame
% element with modulus of elasticity E,
% shear modulus of elasticity G, cross-
% sectional area A, moments of inertia
% Iy and Iz, torsional constant J,
% coordinates (x1,y1,z1) for the first
% node and coordinates (x2,y2,z2) for the
% second node.
% The size of the element stiffness
% matrix is 12 x 12.
L = sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) + (z2-z1)*(z2-z1));
w1 = E*A/L;
w2 = 12*E*Iz/(L*L*L);
w3 = 6*E*Iz/(L*L);
w4 = 4*E*Iz/L;
w5 = 2*E*Iz/L;
w6 = 12*E*Iy/(L*L*L);
w7 = 6*E*Iy/(L*L);
w8 = 4*E*Iy/L;
w9 = 2*E*Iy/L;
w10 = G*J/L;
kprime = [w1 0 0 0 0 0 -w1 0 0 0 0 0 ;
0 w2 0 0 0 w3 0 -w2 0 0 0 w3 ;
0 0 w6 0 -w7 0 0 0 -w6 0 -w7 0 ;
0 0 0 w10 0 0 0 0 0 -w10 0 0 ;
0 0 -w7 0 w8 0 0 0 w7 0 w9 0 ;
0 w3 0 0 0 w4 0 -w3 0 0 0 w5 ;
-w1 0 0 0 0 0 w1 0 0 0 0 0 ;
0 -w2 0 0 0 -w3 0 w2 0 0 0 -w3 ;
0 0 -w6 0 w7 0 0 0 w6 0 w7 0 ;
0 0 0 -w10 0 0 0 0 0 w10 0 0 ;
0 0 -w7 0 w9 0 0 0 w7 0 w8 0 ;
0 w3 0 0 0 w5 0 -w3 0 0 0 w4];
if x1 == x2 & y1 == y2
if z2 > z1
Lambda = [0 0 1 ; 0 1 0 ; -1 0 0];
else
Lambda = [0 0 -1 ; 0 1 0 ; 1 0 0];
end
else
CXx = (x2-x1)/L;
CYx = (y2-y1)/L;
CZx = (z2-z1)/L;
D = sqrt(CXx*CXx + CYx*CYx);
CXy = -CYx/D;
CYy = CXx/D;
CZy = 0;
CXz = -CXx*CZx/D;
CYz = -CYx*CZx/D;
CZz = D;
Lambda = [CXx CYx CZx ; CXy CYy CZy ; CXz CYz CZz];
end
R = [Lambda zeros(3) zeros(3) zeros(3) ;
zeros(3) Lambda zeros(3) zeros(3) ;
zeros(3) zeros(3) Lambda zeros(3) ;
zeros(3) zeros(3) zeros(3) Lambda];
y = R''*kprime*R;
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
function y = SpaceFrameAssemble(K,k,i,j)
%SpaceFrameAssemble This function assembles the element stiffness
% matrix k of the space 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(6*i-5,6*i-5) = K(6*i-5,6*i-5) + k(1,1);
K(6*i-5,6*i-4) = K(6*i-5,6*i-4) + k(1,2);
K(6*i-5,6*i-3) = K(6*i-5,6*i-3) + k(1,3);
K(6*i-5,6*i-2) = K(6*i-5,6*i-2) + k(1,4);
K(6*i-5,6*i-1) = K(6*i-5,6*i-1) + k(1,5);
K(6*i-5,6*i) = K(6*i-5,6*i) + k(1,6);
K(6*i-5,6*j-5) = K(6*i-5,6*j-5) + k(1,7);
K(6*i-5,6*j-4) = K(6*i-5,6*j-4) + k(1,8);
K(6*i-5,6*j-3) = K(6*i-5,6*j-3) + k(1,9);
K(6*i-5,6*j-2) = K(6*i-5,6*j-2) + k(1,10);
K(6*i-5,6*j-1) = K(6*i-5,6*j-1) + k(1,11);
K(6*i-5,6*j) = K(6*i-5,6*j) + k(1,12);
K(6*i-4,6*i-5) = K(6*i-4,6*i-5) + k(2,1);
K(6*i-4,6*i-4) = K(6*i-4,6*i-4) + k(2,2);
K(6*i-4,6*i-3) = K(6*i-4,6*i-3) + k(2,3);
K(6*i-4,6*i-2) = K(6*i-4,6*i-2) + k(2,4);
K(6*i-4,6*i-1) = K(6*i-4,6*i-1) + k(2,5);
K(6*i-4,6*i) = K(6*i-4,6*i) + k(2,6);
K(6*i-4,6*j-5) = K(6*i-4,6*j-5) + k(2,7);
K(6*i-4,6*j-4) = K(6*i-4,6*j-4) + k(2,8);
K(6*i-4,6*j-3) = K(6*i-4,6*j-3) + k(2,9);
K(6*i-4,6*j-2) = K(6*i-4,6*j-2) + k(2,10);
K(6*i-4,6*j-1) = K(6*i-4,6*j-1) + k(2,11);
K(6*i-4,6*j) = K(6*i-4,6*j) + k(2,12);
K(6*i-3,6*i-5) = K(6*i-3,6*i-5) + k(3,1);
K(6*i-3,6*i-4) = K(6*i-3,6*i-4) + k(3,2);
K(6*i-3,6*i-3) = K(6*i-3,6*i-3) + k(3,3);
K(6*i-3,6*i-2) = K(6*i-3,6*i-2) + k(3,4);
K(6*i-3,6*i-1) = K(6*i-3,6*i-1) + k(3,5);
K(6*i-3,6*i) = K(6*i-3,6*i) + k(3,6);
K(6*i-3,6*j-5) = K(6*i-3,6*j-5) + k(3,7);
K(6*i-3,6*j-4) = K(6*i-3,6*j-4) + k(3,8);
K(6*i-3,6*j-3) = K(6*i-3,6*j-3) + k(3,9);
K(6*i-3,6*j-2) = K(6*i-3,6*j-2) + k(3,10);
K(6*i-3,6*j-1) = K(6*i-3,6*j-1) + k(3,11);
K(6*i-3,6*j) = K(6*i-3,6*j) + k(3,12);
K(6*i-2,6*i-5) = K(6*i-2,6*i-5) + k(4,1);
K(6*i-2,6*i-4) = K(6*i-2,6*i-4) + k(4,2);
K(6*i-2,6*i-3) = K(6*i-2,6*i-3) + k(4,3);
K(6*i-2,6*i-2) = K(6*i-2,6*i-2) + k(4,4);
K(6*i-2,6*i-1) = K(6*i-2,6*i-1) + k(4,5);
K(6*i-2,6*i) = K(6*i-2,6*i) + k(4,6);
K(6*i-2,6*j-5) = K(6*i-2,6*j-5) + k(4,7);
K(6*i-2,6*j-4) = K(6*i-2,6*j-4) + k(4,8);
K(6*i-2,6*j-3) = K(6*i-2,6*j-3) + k(4,9);
K(6*i-2,6*j-2) = K(6*i-2,6*j-2) + k(4,10);
K(6*i-2,6*j-1) = K(6*i-2,6*j-1) + k(4,11);
K(6*i-2,6*j) = K(6*i-2,6*j) + k(4,12);
K(6*i-1,6*i-5) = K(6*i-1,6*i-5) + k(5,1);
K(6*i-1,6*i-4) = K(6*i-1,6*i-4) + k(5,2);
K(6*i-1,6*i-3) = K(6*i-1,6*i-3) + k(5,3);
K(6*i-1,6*i-2) = K(6*i-1,6*i-2) + k(5,4);
K(6*i-1,6*i-1) = K(6*i-1,6*i-1) + k(5,5);
K(6*i-1,6*i) = K(6*i-1,6*i) + k(5,6);
K(6*i-1,6*j-5) = K(6*i-1,6*j-5) + k(5,7);
K(6*i-1,6*j-4) = K(6*i-1,6*j-4) + k(5,8);
K(6*i-1,6*j-3) = K(6*i-1,6*j-3) + k(5,9);
K(6*i-1,6*j-2) = K(6*i-1,6*j-2) + k(5,10);
K(6*i-1,6*j-1) = K(6*i-1,6*j-1) + k(5,11);
K(6*i-1,6*j) = K(6*i-1,6*j) + k(5,12);
K(6*i,6*i-5) = K(6*i,6*i-5) + k(6,1);
K(6*i,6*i-4) = K(6*i,6*i-4) + k(6,2);
K(6*i,6*i-3) = K(6*i,6*i-3) + k(6,3);
K(6*i,6*i-2) = K(6*i,6*i-2) + k(6,4);
K(6*i,6*i-1) = K(6*i,6*i-1) + k(6,5);
K(6*i,6*i) = K(6*i,6*i) + k(6,6);
K(6*i,6*j-5) = K(6*i,6*j-5) + k(6,7);
K(6*i,6*j-4) = K(6*i,6*j-4) + k(6,8);
K(6*i,6*j-3) = K(6*i,6*j-3) + k(6,9);
K(6*i,6*j-2) = K(6*i,6*j-2) + k(6,10);
K(6*i,6*j-1) = K(6*i,6*j-1) + k(6,11);
K(6*i,6*j) = K(6*i,6*j) + k(6,12);
K(6*j-5,6*i-5) = K(6*j-5,6*i-5) + k(7,1);
K(6*j-5,6*i-4) = K(6*j-5,6*i-4) + k(7,2);
K(6*j-5,6*i-3) = K(6*j-5,6*i-3) + k(7,3);
K(6*j-5,6*i-2) = K(6*j-5,6*i-2) + k(7,4);
K(6*j-5,6*i-1) = K(6*j-5,6*i-1) + k(7,5);
K(6*j-5,6*i) = K(6*j-5,6*i) + k(7,6);
K(6*j-5,6*j-5) = K(6*j-5,6*j-5) + k(7,7);
K(6*j-5,6*j-4) = K(6*j-5,6*j-4) + k(7,8);
K(6*j-5,6*j-3) = K(6*j-5,6*j-3) + k(7,9);
K(6*j-5,6*j-2) = K(6*j-5,6*j-2) + k(7,10);
K(6*j-5,6*j-1) = K(6*j-5,6*j-1) + k(7,11);
K(6*j-5,6*j) = K(6*j-5,6*j) + k(7,12);
K(6*j-4,6*i-5) = K(6*j-4,6*i-5) + k(8,1);
K(6*j-4,6*i-4) = K(6*j-4,6*i-4) + k(8,2);
K(6*j-4,6*i-3) = K(6*j-4,6*i-3) + k(8,3);
K(6*j-4,6*i-2) = K(6*j-4,6*i-2) + k(8,4);
K(6*j-4,6*i-1) = K(6*j-4,6*i-1) + k(8,5);
K(6*j-4,6*i) = K(6*j-4,6*i) + k(8,6);
K(6*j-4,6*j-5) = K(6*j-4,6*j-5) + k(8,7);
K(6*j-4,6*j-4) = K(6*j-4,6*j-4) + k(8,8);
K(6*j-4,6*j-3) = K(6*j-4,6*j-3) + k(8,9);
K(6*j-4,6*j-2) = K(6*j-4,6*j-2) + k(8,10);
K(6*j-4,6*j-1) = K(6*j-4,6*j-1) + k(8,11);
K(6*j-4,6*j) = K(6*j-4,6*j) + k(8,12);
K(6*j-3,6*i-5) = K(6*j-3,6*i-5) + k(9,1);
K(6*j-3,6*i-4) = K(6*j-3,6*i-4) + k(9,2);
K(6*j-3,6*i-3) = K(6*j-3,6*i-3) + k(9,3);
K(6*j-3,6*i-2) = K(6*j-3,6*i-2) + k(9,4);
K(6*j-3,6*i-1) = K(6*j-3,6*i-1) + k(9,5);
K(6*j-3,6*i) = K(6*j-3,6*i) + k(9,6);
K(6*j-3,6*j-5) = K(6*j-3,6*j-5) + k(9,7);
K(6*j-3,6*j-4) = K(6*j-3,6*j-4) + k(9,8);
K(6*j-3,6*j-3) = K(6*j-3,6*j-3) + k(9,9);
K(6*j-3,6*j-2) = K(6*j-3,6*j-2) + k(9,10);
K(6*j-3,6*j-1) = K(6*j-3,6*j-1) + k(9,11);
K(6*j-3,6*j) = K(6*j-3,6*j) + k(9,12);
K(6*j-2,6*i-5) = K(6*j-2,6*i-5) + k(10,1);
K(6*j-2,6*i-4) = K(6*j-2,6*i-4) + k(10,2);
K(6*j-2,6*i-3) = K(6*j-2,6*i-3) + k(10,3);
K(6*j-2,6*i-2) = K(6*j-2,6*i-2) + k(10,4);
K(6*j-2,6*i-1) = K(6*j-2,6*i-1) + k(10,5);
K(6*j-2,6*i) = K(6*j-2,6*i) + k(10,6);
K(6*j-2,6*j-5) = K(6*j-2,6*j-5) + k(10,7);
K(6*j-2,6*j-4) = K(6*j-2,6*j-4) + k(10,8);
K(6*j-2,6*j-3) = K(6*j-2,6*j-3) + k(10,9);
K(6*j-2,6*j-2) = K(6*j-2,6*j-2) + k(10,10);
K(6*j-2,6*j-1) = K(6*j-2,6*j-1) + k(10,11);
K(6*j-2,6*j) = K(6*j-2,6*j) + k(10,12);
K(6*j-1,6*i-5) = K(6*j-1,6*i-5) + k(11,1);
K(6*j-1,6*i-4) = K(6*j-1,6*i-4) + k(11,2);
K(6*j-1,6*i-3) = K(6*j-1,6*i-3) + k(11,3);
K(6*j-1,6*i-2) = K(6*j-1,6*i-2) + k(11,4);
K(6*j-1,6*i-1) = K(6*j-1,6*i-1) + k(11,5);
K(6*j-1,6*i) = K(6*j-1,6*i) + k(11,6);
K(6*j-1,6*j-5) = K(6*j-1,6*j-5) + k(11,7);
K(6*j-1,6*j-4) = K(6*j-1,6*j-4) + k(11,8);
K(6*j-1,6*j-3) = K(6*j-1,6*j-3) + k(11,9);
K(6*j-1,6*j-2) = K(6*j-1,6*j-2) + k(11,10);
K(6*j-1,6*j-1) = K(6*j-1,6*j-1) + k(11,11);
K(6*j-1,6*j) = K(6*j-1,6*j) + k(11,12);
K(6*j,6*i-5) = K(6*j,6*i-5) + k(12,1);
K(6*j,6*i-4) = K(6*j,6*i-4) + k(12,2);
K(6*j,6*i-3) = K(6*j,6*i-3) + k(12,3);
K(6*j,6*i-2) = K(6*j,6*i-2) + k(12,4);
K(6*j,6*i-1) = K(6*j,6*i-1) + k(12,5);
K(6*j,6*i) = K(6*j,6*i) + k(12,6);
K(6*j,6*j-5) = K(6*j,6*j-5) + k(12,7);
K(6*j,6*j-4) = K(6*j,6*j-4) + k(12,8);
K(6*j,6*j-3) = K(6*j,6*j-3) + k(12,9);
K(6*j,6*j-2) = K(6*j,6*j-2) + k(12,10);
K(6*j,6*j-1) = K(6*j,6*j-1) + k(12,11);
K(6*j,6*j) = K(6*j,6*j) + k(12,12);
y = K;
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
function y = SpaceFrameElementForces(E,G,A,Iy,Iz,J,x1,y1,z1,x2,y2,z2,u)
%SpaceFrameElementForces This function returns the element force
% vector given the modulus of elasticity E,
% the shear modulus of elasticity G, the
% cross-sectional area A, moments of inertia
% Iy and Iz, the torsional constant J,
% the coordinates (x1,y1,z1) of the first
% node, the coordinates (x2,y2,z2) of the
% second node, and the element nodal
% displacement vector u.
L = sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) + (z2-z1)*(z2-z1));
w1 = E*A/L;
w2 = 12*E*Iz/(L*L*L);
w3 = 6*E*Iz/(L*L);
w4 = 4*E*Iz/L;
w5 = 2*E*Iz/L;
w6 = 12*E*Iy/(L*L*L);
w7 = 6*E*Iy/(L*L);
w8 = 4*E*Iy/L;
w9 = 2*E*Iy/L;
w10 = G*J/L;
kprime = [w1 0 0 0 0 0 -w1 0 0 0 0 0 ;
0 w2 0 0 0 w3 0 -w2 0 0 0 w3 ;
0 0 w6 0 -w7 0 0 0 -w6 0 -w7 0 ;
0 0 0 w10 0 0 0 0 0 -w10 0 0 ;
0 0 -w7 0 w8 0 0 0 w7 0 w9 0 ;
0 w3 0 0 0 w4 0 -w3 0 0 0 w5 ;
-w1 0 0 0 0 0 w1 0 0 0 0 0 ;
0 -w2 0 0 0 -w3 0 w2 0 0 0 -w3 ;
0 0 -w6 0 w7 0 0 0 w6 0 w7 0 ;
0 0 0 -w10 0 0 0 0 0 w10 0 0 ;
0 0 -w7 0 w9 0 0 0 w7 0 w8 0 ;
0 w3 0 0 0 w5 0 -w3 0 0 0 w4];
if x1 == x2 & y1 == y2
if z2 > z1
Lambda = [0 0 1 ; 0 1 0 ; -1 0 0];
else
Lambda = [0 0 -1 ; 0 1 0 ; 1 0 0];
end
else
CXx = (x2-x1)/L;
CYx = (y2-y1)/L;
CZx = (z2-z1)/L;
D = sqrt(CXx*CXx + CYx*CYx);
CXy = -CYx/D;
CYy = CXx/D;
CZy = 0;
CXz = -CXx*CZx/D;
CYz = -CYx*CZx/D;
CZz = D;
Lambda = [CXx CYx CZx ; CXy CYy CZy ; CXz CYz CZz];
end
R = [Lambda zeros(3) zeros(3) zeros(3) ;
zeros(3) Lambda zeros(3) zeros(3) ;
zeros(3) zeros(3) Lambda zeros(3) ;
zeros(3) zeros(3) zeros(3) Lambda];
y = kprime*R* u;
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
function y = SpaceFrameElementAxialDiagram(f, L)
%SpaceFrameElementAxialDiagram This function plots the axial force
% diagram for the space frame element
% with nodal force vector f and length
% L.
x = [0 ; L];
z = [-f(1) ; f(7)];
hold on;
title(''Axial Force Diagram'');
plot(x,z);
y1 = [0 ; 0];
plot(x,y1,''k'')
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
function y = SpaceFrameElementShearZDiagram(f, L)
%SpaceFrameElementShearZDiagram This function plots the shear force
% diagram for the space frame element
% with nodal force vector f and
% length L.
x = [0 ; L];
z = [f(3) ; -f(9)];
hold on;
title(''Shear Force Diagram in Z Direction'');
plot(x,z);
y1 = [0 ; 0];
plot(x,y1,''k'')
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
function y = SpaceFrameElementShearYDiagram(f, L)
%SpaceFrameElementShearYDiagram This function plots the shear force
% diagram for the space frame element
% with nodal force vector f and
% length L.
x = [0 ; L];
z = [f(2) ; -f(8)];
hold on;
title(''Shear Force Diagram in Y Direction'');
plot(x,z);
y1 = [0 ; 0];
plot(x,y1,''k'')
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
function y = SpaceFrameElementTorsionDiagram(f, L)
%SpaceFrameElementTorsionDiagram This function plots the torsion
% diagram for the space frame
% element with nodal force vector f
% and length L.
x = [0 ; L];
z = [f(4) ; -f(10)];
hold on;
title(''Torsion Diagram'');
plot(x,z);
y1 = [0 ; 0];
plot(x,y1,''k'')
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
function y = SpaceFrameElementMomentZDiagram(f, L)
%SpaceFrameElementMomentZDiagram This function plots the bending
% moment diagram for the space frame
% element with nodal force vector f
% and length L.
x = [0 ; L];
z = [f(6) ; -f(12)];
hold on;
title(''Bending Moment Diagram along Z Axis'');
plot(x,z);
y1 = [0 ; 0];
plot(x,y1,''k'')
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
function y = SpaceFrameElementMomentYDiagram(f, L)
%SpaceFrameElementMomentYDiagram This function plots the bending
% moment diagram for the space frame
% element with nodal force vector f
% and length L.
x = [0 ; L];
z = [f(5) ; -f(11)];
hold on;
title(''Bending Moment Diagram along Y Axis'');
plot(x,z);
y1 = [0 ; 0];
plot(x,y1,''k'')

The Linear Triangular Element
---------------------------------------------------------------------------------------
function y = LinearTriangleElementArea(xi,yi,xj,yj,xm,ym)
%LinearTriangleElementArea This function returns the area of the
% linear triangular element whose first
% node has coordinates (xi,yi), second
% node has coordinates (xj,yj), and
% third node has coordinates (xm,ym).
y = (xi*(yj-ym) + xj*(ym-yi) + xm*(yi-yj))/2;
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
function y = LinearTriangleElementStiffness(E,NU,t,xi,yi,xj,yj,xm,ym,p)
%LinearTriangleElementStiffness This function returns the element
% stiffness matrix for a linear
% triangular element with modulus of
% elasticity E, Poisson''s ratio NU,
% thickness t, coordinates of the
% first node (xi,yi), coordinates of
% the second node (xj,yj), and
% coordinates of the third node
% (xm,ym). Use p = 1 for cases of
% plane stress, and p = 2 for cases
% of plane strain.
% The size of the element stiffness
% matrix is 6 x 6.
A = (xi*(yj-ym) + xj*(ym-yi) + xm*(yi-yj))/2;
betai = yj-ym;
betaj = ym-yi;
betam = yi-yj;
gammai = xm-xj;
gammaj = xi-xm;
gammam = xj-xi;
B = [betai 0 betaj 0 betam 0 ;
0 gammai 0 gammaj 0 gammam ;
gammai betai gammaj betaj gammam betam]/(2*A);
if p == 1
D = (E/(1-NU*NU))*[1 NU 0 ; NU 1 0 ; 0 0 (1-NU)/2];
elseif p == 2
D = (E/(1+NU)/(1-2*NU))*[1-NU NU 0 ; NU 1-NU 0 ; 0 0 (1-2*NU)/2];
end
y = t*A*B''*D*B;
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
unction y = LinearTriangleAssemble(K,k,i,j,m)
%LinearTriangleAssemble This function assembles the element
% stiffness matrix k of the linear
% triangular element 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(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-1,2*m-1) = K(2*i-1,2*m-1) + k(1,5);
K(2*i-1,2*m) = K(2*i-1,2*m) + k(1,6);
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*i,2*m-1) = K(2*i,2*m-1) + k(2,5);
K(2*i,2*m) = K(2*i,2*m) + k(2,6);
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-1,2*m-1) = K(2*j-1,2*m-1) + k(3,5);
K(2*j-1,2*m) = K(2*j-1,2*m) + k(3,6);
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);
K(2*j,2*m-1) = K(2*j,2*m-1) + k(4,5);
K(2*j,2*m) = K(2*j,2*m) + k(4,6);
K(2*m-1,2*i-1) = K(2*m-1,2*i-1) + k(5,1);
K(2*m-1,2*i) = K(2*m-1,2*i) + k(5,2);
K(2*m-1,2*j-1) = K(2*m-1,2*j-1) + k(5,3);
K(2*m-1,2*j) = K(2*m-1,2*j) + k(5,4);
K(2*m-1,2*m-1) = K(2*m-1,2*m-1) + k(5,5);
K(2*m-1,2*m) = K(2*m-1,2*m) + k(5,6);
K(2*m,2*i-1) = K(2*m,2*i-1) + k(6,1);
K(2*m,2*i) = K(2*m,2*i) + k(6,2);
K(2*m,2*j-1) = K(2*m,2*j-1) + k(6,3);
K(2*m,2*j) = K(2*m,2*j) + k(6,4);
K(2*m,2*m-1) = K(2*m,2*m-1) + k(6,5);
K(2*m,2*m) = K(2*m,2*m) + k(6,6);
y = K;
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
function y = LinearTriangleElementPStresses(sigma)
%LinearTriangleElementPStresses This function returns the element
% principal stresses and their
% angle given the element
% stress vector.
R = (sigma(1) + sigma(2))/2;
Q = ((sigma(1) - sigma(2))/2)^2 + sigma(3)*sigma(3);
M = 2*sigma(3)/(sigma(1) - sigma(2));
s1 = R + sqrt(Q);
s2 = R - sqrt(Q);
theta = (atan(M)/2)*180/pi;
y = [s1 ; s2 ; theta];
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
function y = LinearTriangleElementStresses(E,NU,t,xi,yi,xj,yj,xm,ym,p,u)
%LinearTriangleElementStresses This function returns the element
% stress vector for a linear
% triangular element with modulus of
% elasticity E, Poisson''s ratio NU,
% thickness t, coordinates of the
% first node (xi,yi), coordinates of
% the second node (xj,yj),
% coordinates of the third node
% (xm,ym), and element displacement
% vector u. Use p = 1 for cases of
% plane stress, and p = 2 for cases
% of plane strain.
% The size of the element stress
% vector is 3 x 1.
A = (xi*(yj-ym) + xj*(ym-yi) + xm*(yi-yj))/2;
betai = yj-ym;
betaj = ym-yi;
betam = yi-yj;
gammai = xm-xj;
gammaj = xi-xm;
gammam = xj-xi;
B = [betai 0 betaj 0 betam 0 ;
0 gammai 0 gammaj 0 gammam ;
gammai betai gammaj betaj gammam betam]/(2*A);
if p == 1
D = (E/(1-NU*NU))*[1 NU 0 ; NU 1 0 ; 0 0 (1-NU)/2];
elseif p == 2
D = (E/(1+NU)/(1-2*NU))*[1-NU NU 0 ; NU 1-NU 0 ; 0 0 (1-2*NU)/2];
end
y = D*B*u;
---------------------------------------------------------------------------------------

Các bác có biết sử dụng MatLab để xử lý số liệu thực nghiệm không ạ. Cụ thể là từ 1 bảng số liệu thực nghiệm, dùng MatLab vẽ đồ thị quan hệ và nội suy ra được phương trình hàm số gần đúng quan hệ giữa các đại lượng. Bác nào biết không ạ, em mong được học hỏi.