www.pudn.com > zhichixiangliangjijiemian.rar > Untitled3.asv

clear;  
% DEFINE BOUNDARIES/PARAMETERS  
Lb=1; 
D=5; 
Q=100; 
Etm=3e7;Ebsb=0.2; 
Stm=2.1e8;Sbsb=0.3; 
 
 
 
%平面应变问题 
Dmat=(tm/(1-2*bsb)/(1+bsb))*[(1-bsb) bsb 0;bsb (1-bsb) 0;0 0 (1-2*bsb)/2]; 
 
% SET UP NODAL COORDINATES 
ndivl=8; 
ndivw=20; 
[x,conn,numcell,numnod] = mesh2(Lb,D,ndivl,ndivw); 
numnod=length(x); 
 
% SET UP QUADRATURE CELLS 
ndivlq = 8;ndivwq = 20; 
[xc,conn,numcell,numq] = mesh2(Lb,D,ndivlq,ndivwq); 
% DETERMINE DOMAINS OF INFLUENCE - UNIFORM NODAL SPACING 
dmax=3.5; 
xspac = Lb/ndivl; 
yspac = D/ndivw; 
dm(1,1:numnod)=dmax*xspac*ones(1,numnod); 
dm(2,1:numnod)=dmax*yspac*ones(1,numnod); 
 
% SET UP GAUSS POINTS, WEIGHTS, AND JACOBIAN FOR EACH CELL 
quado = 4; 
[gauss] = gauss2(quado); 
numq2 = numcell*quado^2; 
gs = zeros(4,numq2); 
[gs] = egauss(xc,conn,gauss,numcell); 
ngs=length(gs); 
 
 
 
% LOOP OVER GAUSS POINTS TO ASSEMBLE DISCRETE EQUATIONS 
k = sparse(numnod*2,numnod*2); 
for gg=gs 
    gpos=gg(1:2,:); 
    weight=gg(3); 
    jac=gg(4); 
% DETERMINE NODES IN NEIGHBORHOOD OF GAUSS POINT 
v = domain(gpos,x,dm,numnod); 
L = length(v); 
en = zeros(1,2*L); 
[phi,dphix,dphiy] = shape(gpos,dmax,x,v,dm); 
Bmat=zeros(3,2*L); 
for j=1:L 
Bmat(1:3,(2*j-1):2*j) = [dphix(j) 0;0 dphiy(j);dphiy(j) dphix(j)]; 
end 
for i=1:L 
en(2*i-1) = 2*v(i)-1; 
en(2*i) = 2*v(i); 
end 
k(en,en) = k(en,en)+sparse((weight*jac)*Bmat'*Dmat*Bmat); 
end 
 
f=21:21:21*9; 
for i=1:length(f) 
    qn(2*i-1)=2*f(i)-1; 
    qn(2*i)=2*f(i); 
    k(qn,qn)=k(qn,qn)*1e005; 
end 
 
 
% DETERMINE NODES ON BOUNDARY, SET UP BC'S 
ind1 = 0;ind2 = 0;ind3=0;ind4=0; 
for j=1:numnod 
if(x(1,j)==-Lb/2) 
ind1=ind1+1; 
nnu(1,ind1) = x(1,j); 
nnu(2,ind1) = x(2,j); 
no1(ind1)=j; 
end 
if(x(1,j)==Lb/2) 
ind2=ind2+1; 
nt(1,ind2) = x(1,j); 
nt(2,ind2) = x(2,j); 
no2(ind2)=j; 
end 
if(x(2,j)==D/2) 
ind3=ind3+1; 
nt2(1,ind3) = x(1,j); 
nt2(2,ind3) = x(2,j); 
no3(ind3)=j; 
end 
if(x(2,j)==-D/2) 
ind4=ind4+1; 
nt1(1,ind4) = x(1,j); 
nt1(2,ind4) = x(2,j); 
no4(ind4)=j; 
end 
end 
lthu = length(nnu); 
ltht = length(nt); 
ltht2=length(nt2); 
ltht1=length(nt1); 
f = zeros(1,numnod*2); 
 
 
%计算zuo节点荷载 
qn=Q/D*yspac; 
qs=1/6*qn*yspac; 
qx=2*qs; 
qk=0.5*qn*yspac; 
 
degao=zeros(3,lthu); 
for i=1:lthu 
    degao(1,i)=qs; 
end 
for i=1:lthu 
    degao(2,i)=qx; 
end 
for i=1:lthu 
    degao(3,i)=qk; 
end 
 
 
ind=0; 
suan=0; 
jisuan=zeros(4,lthu); 
for i=1:lthu 
    ind=ind+1; 
    suan=suan+degao(3,i); 
    jisuan(1,ind)=suan; 
end 
for i=3:lthu-1 
    ind=ind+3; 
    suan=suan+degao(3,i); 
   jisuan(2,ind)=suan; 
end 
for i=1:lthu-1 
    jisuan(3,i)=degao(1,i); 
end 
for i=2:lthu 
    jisuan(4,i)=degao(2,i); 
end 
 
suanmain=zeros(2,lthu); 
for i=1:lthu 
    suanmain(1,i)=jisuan(1,i)+jisuan(2,i)+jisuan(3,i)+jisuan(4,i); 
end 
 
for i=1:lthu 
    f(2*i-1)=suanmain(1,i); 
    f(2*i)=suanmain(2,i); 
end 
f; 
 
%求解节点位移 
d=k\f'; 
u=d(1:2*numnod); 
for i=1:numnod 
    uf(1,i)=u(2*i-1); 
    uf(2,i)=u(2*i); 
end 
uf 
 
 
%求解节点应力 
ind=0; 
enorm=0; 
for gg=x 
ind=ind+1; 
gpos=gg(1:2); 
v = domain(gpos,x,dm,numnod); 
L=length(v); 
en=zeros(1,2); 
[phi,dphix,dphiy]=shape(gpos,dmax,x,v,dm); 
Bmat=zeros(3,2); 
for j=1:L 
Bmat(1:3,(2*j-1):2*j)=[dphix(j) 0;0 dphiy(j);dphiy(j) dphix(j)]; 
en(2*j-1)=2*v(j)-1; 
en(2*j)=2*v(j); 
stress1(1:3,ind)=Dmat*Bmat*u(en); 
end 
end 
x; 
stress1 
 
 
 
%绘制挡土墙体应力云图 
figure(1) 
x1=-Lb/2:xspac:Lb/2; 
x2=linspace(D/2,-D/2,ndivw+1); 
z12=stress1(2,:); 
z12=reshape(z12,ndivw+1,length(x1)); 
pcolor(x1,x2,z12); 
shading interp; 
xlabel('挡土墙体各点X坐标'); 
ylabel('挡土墙体各点Y坐标'); 
title('挡土墙体各点σy应力云图'); 
axis equal; 
 
 
figure(2) 
x3=-Lb/2:xspac:Lb/2; 
x4=linspace(D/2,-D/2,ndivw+1); 
z34=stress1(1,:); 
z34=reshape(z34,ndivw+1,length(x1)); 
pcolor(x3,x4,z34); 
shading interp; 
xlabel('挡土墙体各点X坐标'); 
ylabel('挡土墙体各点Y坐标'); 
title('挡土墙体各点σx应力云图'); 
axis equal; 
 
 
figure(3) 
xb=[-0.5 0.5 0.5 -0.5 -0.5]; 
yb=[-2.5 -2.5 2.5 2.5 -2.5]; 
plot(xb,yb,'k','LineWidth',1.5); 
hold on 
xb1=x(:,no1)+1000*uf(:,no1); 
xb2=x(:,no2)+1000*uf(:,no2); 
xb3=x(:,no3)+1000*uf(:,no3); 
xb4=x(:,no4)+1000*uf(:,no4); 
plot(xb1(1,:),xb1(2,:),'r',xb2(1,:),xb2(2,:),'r',xb3(1,:),xb3(2,:),'r',xb3(1,:),xb3(2,:),'r','LineWidth',1.5); 
axis equal;