165 lines
4.9 KiB
Matlab
165 lines
4.9 KiB
Matlab
clear all;clc;close all;
|
||
% global Elements gKA
|
||
%标准国际单位制
|
||
Length=2;%2m
|
||
min_H=0.2;max_H=0.6;%截面高度
|
||
Width=1;%截面宽度
|
||
E=3e10;%Pa 3e4MPa
|
||
u=0.2;% 泊松比
|
||
q0=1000;%均布力；1kN/m
|
||
F_jizhong=2e3;%集中力2kN
|
||
%理论解
|
||
% EI=E*Width*min_H^3/12;disp=[q0*Length^4/(8*EI)+F_jizhong*Length^3/(3*EI)]*1000;
|
||
%定义节点和单元
|
||
Ne_x=16;%单元列数4
|
||
Ne_y=8;%单元行数2
|
||
D=[E/(1-u^2), u*E/(1-u^2), 0;
|
||
u*E/(1-u^2), E/(1-u^2), 0;
|
||
0, 0, E/(2*(1+u))];
|
||
|
||
Ne=Ne_x*Ne_y; % 单元总数
|
||
Nnod=(Ne_x+1)*(2*Ne_y+1)+Ne_x*(Ne_y+1); % 节点总数
|
||
%节点坐标
|
||
x=linspace(-Length/2,Length/2,2*Ne_x+1); %节点从左到右的横坐标
|
||
rt=[repmat([2*Ne_y+1 Ne_y+1],1,Ne_x) 2*Ne_y+1]; % 每列节点数 5 3 5 3 5 3 ……
|
||
x=repelem(x,rt); % x坐标扩展为所有节点的坐标
|
||
y=[];
|
||
dy=(max_H-min_H)/Ne_x/4;
|
||
for i=1:Ne_x
|
||
y=[y,linspace(-max_H/2+(i-1)*2*dy,max_H/2-(i-1)*2*dy,2*Ne_y+1),linspace(-max_H/2+[(i-1)*2+1]*dy,max_H/2-[(i-1)*2+1]*dy,Ne_y+1)];
|
||
end
|
||
y=[y,linspace(-min_H/2,min_H/2,2*Ne_y+1)];
|
||
Nodes=[x;y]';
|
||
|
||
% 单元节点编号
|
||
Elements=zeros(Ne,8);k=0;
|
||
for i=1:Ne_y:Ne
|
||
Elements(i:i+Ne_y-1,1)=(1+k*(3*Ne_y+2)):2:(1+k*(3*Ne_y+2)+2*(Ne_y-1));
|
||
Elements(i:i+Ne_y-1,2)=Elements(i:i+Ne_y-1,1)+2;
|
||
Elements(i:i+Ne_y-1,3)=Elements(i:i+Ne_y-1,2)+3*Ne_y+2;
|
||
Elements(i:i+Ne_y-1,4)=Elements(i:i+Ne_y-1,3)-2;
|
||
Elements(i:i+Ne_y-1,5)=2+k*(3*Ne_y+2):2:2+k*(3*Ne_y+2)+2*(Ne_y-1);
|
||
Elements(i:i+Ne_y-1,6)=2*Ne_y+3+k*(3*Ne_y+2):2*Ne_y+3+k*(3*Ne_y+2)+Ne_y-1;
|
||
Elements(i:i+Ne_y-1,7)=Elements(i:i+Ne_y-1,3)-1;
|
||
Elements(i:i+Ne_y-1,8)=Elements(i:i+Ne_y-1,6)-1;
|
||
k=k+1;
|
||
end
|
||
%定义分布荷载
|
||
E_ID_p1 = Ne_y:Ne_y:Ne_x*Ne_y; % 分布荷载作用单元编号
|
||
N_ID_p1 = zeros(size(E_ID_p1,2),size(Elements,2));%分布荷载所在单元节点编号
|
||
j = 1;
|
||
for i_q = E_ID_p1
|
||
N_ID_p1(j,:) = Elements(i_q,:); %分布荷载所在单元节点编号
|
||
j = j+1;
|
||
end
|
||
|
||
%施加集中荷载的节点编号
|
||
E_ID_p2 = max(Ne_y:Ne_y:Ne); % 集中荷载作用单元编号***
|
||
Nod_1=Elements(E_ID_p2,4); % 集中荷载节点编号***
|
||
|
||
%构件跨中节点编号
|
||
Elt_p3=max(Ne_y:Ne_y:Ne/2); % 构件跨中单元编号
|
||
Nod_mid=Elements(Elt_p3,3) ; % 构件跨中节点编号
|
||
|
||
%构件左端节点编号
|
||
Elt_p4=Ne_y ; % 构件左端单元编号
|
||
Nod_left=Elements(Elt_p4,2) ; % 构件左端节点编号
|
||
|
||
%边界条件：type = 1 -- fixed disp [ 1--节点固定; 0--节点自由 ];
|
||
% type = 2 -- point load
|
||
% node, x, y, z, rx, ry, rz, type
|
||
|
||
Boundary = [1 1.0 1.0 0.0 0.0 0.0 0.0 1;
|
||
Nod_left 1.0 1.0 0.0 0.0 0.0 0.0 1;
|
||
Nod_1 0.0 F_jizhong 0.0 0.0 0.0 0.0 2 ];
|
||
|
||
% 初始化结果存储空间
|
||
K = zeros( 2 * Nnod , 2 * Nnod );
|
||
F = zeros( 2 * Nnod , 1 );
|
||
% gKA = zeros( 2 * Nnod , 2 * Nnod );
|
||
%计算刚度矩阵
|
||
for i=1:Ne
|
||
ke=ele_stiff_matrix(i,Nodes,Elements,E,u,Width);
|
||
ENodes( : , 1 ) = Elements( i , : );
|
||
% 组装为整体刚度阵
|
||
iEnum = size( ENodes , 1 ); % 单元的结点总数
|
||
iDOFnum = 2 * iEnum; % 单元的总自由度数
|
||
iDof = zeros( 2 * iEnum , 1 ); % 各结点自由度编号
|
||
|
||
for il=1:iEnum
|
||
inod = ENodes( il );
|
||
iDof( 2 * il-1 : 2 * il) = 2 * inod-1 : 2 * inod;
|
||
end
|
||
% 拼装：ir为行局部编号, irg为对应整体编号
|
||
for ir = 1 : iDOFnum
|
||
irg = iDof( ir );
|
||
% ic为列局部编号, icg为对应整体编号
|
||
for ic = 1 : iDOFnum
|
||
icg = iDof( ic );
|
||
K( irg , icg ) = K( irg , icg ) + ke( ir , ic );
|
||
end
|
||
end
|
||
% AssembleEleStif( i , ke ); % 3绾у瓙绋嬪簭锛屾暣浣撳畾浣?
|
||
|
||
%%%%%%%%%%%%
|
||
end
|
||
% aa=gKA-K;
|
||
%施加荷载边界条件
|
||
%均布荷载
|
||
iLenElt = size(N_ID_p1,1); %均布荷载所在单元编号
|
||
for ie = 1:iLenElt
|
||
for ind = 1:2 % ind结点自由度编号
|
||
iNode = N_ID_p1(ie,:); %获取单元结点
|
||
irow = (iNode-1)*2+ind;
|
||
Pe = UniLoad(ie,N_ID_p1,q0,Nodes,Elements); % 3级子程序：EquivalentLoad(ie)
|
||
F(irow) = F(irow)+Pe(ind:2:end,:);
|
||
end
|
||
end
|
||
|
||
|
||
% 集中荷载
|
||
iBct = size( Boundary , 1 );
|
||
for ib = 1 : iBct
|
||
iNode = Boundary( ib , 1 );
|
||
iType = Boundary( ib , 8 );
|
||
if iType == 2 % 点荷载边界条件
|
||
for ind = 1 : 2
|
||
irow = ( iNode-1 ) * 2 + ind;
|
||
F( irow ) = F( irow ) + Boundary( ib , ind + 1 );
|
||
end
|
||
end
|
||
end
|
||
|
||
|
||
|
||
|
||
%施加支座边界条件
|
||
for ib = 1 : iBct
|
||
iNode = Boundary( ib , 1 ); % 结点编号
|
||
iType = Boundary( ib , 8 );
|
||
if iType == 1 % 固定点边界条件
|
||
for ind = 1 : 2
|
||
irow = ( iNode-1 ) * 2 + ind;
|
||
if ( Boundary( ib , ind + 1 ) == 1 )
|
||
K( irow , : ) = 0.0;
|
||
K( : , irow ) = 0.0;
|
||
K( irow , irow ) = 1.0; % 划零置一法
|
||
F(irow)=0.0;
|
||
% gKA( irow , irow ) = 1.0e+20; % 乘大数法
|
||
end
|
||
end
|
||
end
|
||
end
|
||
|
||
Disp = K \ F;
|
||
% 存储结点位移
|
||
gNTu = zeros( Nnod , 2 );
|
||
for id = 1 : Nnod
|
||
gNTu( id , 1 : 2 ) = Disp( 2 * id - 1 : 2 * id );
|
||
end
|
||
% 绘制变形图
|
||
SF=1000;
|
||
PlotResults2(Nodes,Elements,gNTu,SF);
|
||
fprintf('自由端最大挠度%fmm\n',max(abs(gNTu(Nnod,2)))*1000);
|
||
% [NodeStrain,NodeStress]=CalculateStrainAndStress2(Disp,D,Nodes,Elements);
|
||
% PlotContour(Nodes,Elements,Disp,NodeStress(1,:)) |