《有限单元法》--王勖成，习题2.12 MATLAB 程序

%% 《有限单元法》--王勖成，习题2.12 MATLAB 程序

% 仅供参考

function fem()

P=[-10 0 10 0 10 0 -10 0]';%结构节点载荷列向量

[K,P]=xQiYi(P,make_K());

solu(K,P);

%% 存储各单元的节点 x、y 值

function [a,b]=get_xy(m)

x=[0 12 0;

0 12 12;];

y=[4 0 0;

4 4 0];

a=x(m,:);

b=y(m,:);

end

%% 求解并显示

function solu(K,P)

n=1000;

JieDianWeiYi=inv(K)*P

% YingLi=WeiYi(3,1)/12*2e5 %x向拉、压时可用

plot([0 12 12 0],[4 4 0 0],'b');hold on;xlim([-1 14]);ylim([-4 8]);

plot(JieDianWeiYi([1 3 5 7],1)'*n+[0 12 12 0],JieDianWeiYi([2 4 6 8],1)'*n+[4 4 0 0],'r');

end

%% 用“对角元素改 1 法” ，引入位移约束并使总的K 不奇异

function [b,a]=xQiYi(P,K)

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

d=[0 0 1 1 1 1 0 0]';%左侧两结点铰支

% 共 4 个节点，顺序为：左上角开始1 ，顺时针排序

% 共 8个位移分量，[dx1 dy1 dx2 dy2 ... ]

% 位移为 0的分量，对应项填 0，其他填非 0 值

% 例：右侧两点铰支，则为 d=[ 1 1 0 0 0 0 1 1]'

for i=1:1:8

if d(i,1)==0

for j=1:1:8

K(i,j)=0;

K(j,i)=0;

end

P(i,1)=0;

K(i,i)=1;

end

a=P;

b=K;

end

%% 组装总的结构刚度矩阵

function K=make_K()

K=[su([1 1],[1 1]) su(0,[1 2]) su([1 2],[1 3]) su([1 3],0);

su(0,[2 1]) su(0,[2 2]) su(0,[2 3]) su(0,0);

su([2 1],[3 1]) su(0,[3 2]) su([2 2],[3 3]) su([2 3],0);

su([3 1],0) su(0,0) su([3 2],0) su([3 3],0)];

end

%% 分块计算矩阵的和

function b=su(a,b)

K1=get_K(1);

K2=get_K(2);

if a==0

K1=zeros(6,6);

a=[1,1];

end

if b==0

K2=zeros(6,6);

b=[1,1];

end

h11=a(1)*2-1; h12=a(1)*2; l11=a(2)*2-1; l12=a(2)*2;

h21=b(1)*2-1; h22=b(1)*2; l21=b(2)*2-1; l22=b(2)*2;

b=[K1(h11,l11)+K2(h21,l21) K1(h11,l12)+K2(h21,l22);

K1(h12,l11)+K2(h22,l21) K1(h12,l12)+K2(h22,l22)];

end

%% 计算 m 单元的刚度矩阵

function K=get_K(m)

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

v=0.3;

D=2e5/(1-v^2)*[ 1 ,v ,0 ;

v ,1 , 0 ;

0 ,0 ,(1-v)/2];

%弹性矩阵

[x,y]=get_xy(m);

B=get_B(x,y);

A=get_A(x,y);

K=B'*D*B*A;

end

%% 计算单元应变矩阵 B

function B=get_B(x,y)

[~,b,c]=get_abc(x,y);

B=zeros(3,6);

[x,y]=get_xy(1);

A=get_A(x,y);

for i=1:1:3

B(:,i*2-1:i*2)=[b(i) 0;

0 c(i);

c(i) b(i)];

end

B=B/(2*A);

end

%% 计算 ai bi ci

function [a,b,c]=get_abc(x , y)

a=zeros(2,3);

b=a;

c=a;

for l=0:1:2

i=l+1;

j=rem(l+1,3)+1;

m=rem(l+2,3)+1;

a(i)=det([x(j) y(j);

x(m) y(m)]);

b(i)= y(j) - y(m);

c(i)= -x(j) + x(m);

end

%% 计算单元面积

function A=get_A(x,y)

A=-det([1 x(1) y(1); 1 x(2) y(2) ; 1 x(3) y(3)])/2;

end