www.pudn.com > zhichixiangliangjijiemian.rar > jiacondition.M
clear;
% DEFINE BOUNDARIES/PARAMETERS
Lb=5;
D=4;
young = 3e7;nu=0.167;
Q1=0;
Q2=0;
Q=100;
% PLANE STRESS DMATRIX
Dmat = (young/(1-nu^2))*[1 nu 0;nu 1 0;0 0 (1-nu)/2];
% SET UP NODAL COORDINATES
ndivl=20;
ndivw=16;
[x,conn1,conn2,numnod] = mesh2(Lb,D,ndivl,ndivw);
%get away points in the hole
ind=0;
for i=1:length(x)
if( x(1,i)<=-1.25)|(x(1,i)>=1.25)|(x(2,i)>=1.5)|(x(2,i)<=-0.25)
ind=ind+1;
xguodu(:,ind)=x(:,i);
end
end
x=xguodu;
numnod=length(x);
% 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 QUADRATURE CELLS
ndivlq = 20;ndivwq = 16;
[xc,conn,numcell,numq] = mesh2(Lb,D,ndivlq,ndivwq);
% 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);
%get away Gauss points in the hole
ind=0;
for i=1:length(gs)
if( gs(1,i)<=-1.25)|(gs(1,i)>=1.25)|(gs(2,i)>=1.5)|(gs(2,i)<=-0.25)
ind=ind+1;
gsguodu(:,ind)=gs(:,i);
end
end
gs=gsguodu;
% 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
% 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(numnod*2,1);
%SET UP GAUSS POINTS ALONG TRACTION BOUNDARY
xs=linspace(-Lb/2,Lb/2,17);
xz=linspace(D/2,-D/2,13);
xshang=[xs;(D/2)*ones(1,length(xs))];
xxia=[xs;(-D/2)*ones(1,length(xs))];
xzuo=[(-Lb/2)*ones(1,length(xz));xz];
xyou=[(Lb/2)*ones(1,length(xz));xz];
lxs=length(xshang);
lxx=length(xxia);
lxz=length(xzuo);
lxy=length(xyou);
ind=0;
gauss=gauss2(quado);
for i=1:(lxy-1)
ycen=(xyou(2,i+1)+xyou(2,i))/2;
jcob = abs((xyou(2,i+1)-xyou(2,i))/2);
for j=1:quado
mark(j) = ycen-gauss(1,j)*jcob;
ind = ind+1;
gst(1,ind)=xyou(1,i);
gst(2,ind)=mark(j);
gst(3,ind)=gauss(2,j);
gst(4,ind)=jcob;
end
end
%SET UP GAUSS POINTS ALONG DISPLACEMENT BOUNDARY
gsu=gst;
gsu(1,1:ind)=ones(1,ind)*(-Lb/2);
%INTEGRATE FORCES ALONG BOUNDARY
t1=1:4:4*(lxy-1);
t2=4:4:4*(lxy-1);
t3=2:4:4*(lxy-1);
t4=3:4:4*(lxy-1);
gst1=zeros(4,length(t1));
gst2=zeros(4,length(t2));
gst3=zeros(4,length(t3));
gst4=zeros(4,length(t4));
for i=1:length(t1)
gst1(:,i)=gst(:,t1(i));
end
for i=1:length(t2)
gst2(:,i)=gst(:,t2(i));
end
for i=1:length(t3)
gst3(:,i)=gst(:,t3(i));
end
for i=1:length(t4)
gst4(:,i)=gst(:,t4(i));
end
gst12=[gst1,gst2];
gst34=[gst3,gst4];
for gt=gst12
gpos = gt(1:2);
weight=gt(3);
jac = gt(4);
v = domain(gpos,x,dm,numnod);
L = length(v);
en = zeros(1,2*L);
force=zeros(1,2*L);
[phi,dphix,dphiy] = shape(gpos,dmax,x,v,dm);
tx=Q1*0.1997*abs(xyou(2,1)-xyou(2,2))/(weight*jcob);
ty=0;
for i=1:L
en(2*i-1) = 2*v(i)-1;
en(2*i) = 2*v(i);
force(2*i-1)=tx*phi(i);
force(2*i) = ty*phi(i);
end
f(en) = f(en) + jac*weight*force';
end
for gt=gst34
gpos = gt(1:2);
weight=gt(3);
jac = gt(4);
v = domain(gpos,x,dm,numnod);
L = length(v);
en = zeros(1,2*L);
force=zeros(1,2*L);
[phi,dphix,dphiy] = shape(gpos,dmax,x,v,dm);
tx=Q1*0.3003*abs(xyou(2,1)-xyou(2,2))/(weight*jcob);
ty=0;
for i=1:L
en(2*i-1) = 2*v(i)-1;
en(2*i) = 2*v(i);
force(2*i-1)=tx*phi(i);
force(2*i) = ty*phi(i);
end
f(en) = f(en) + jac*weight*force';
end
gsu1=zeros(4,length(t1));
gsu2=zeros(4,length(t2));
gsu3=zeros(4,length(t3));
gsu4=zeros(4,length(t4));
for i=1:length(t1)
gsu1(:,i)=gsu(:,t1(i));
end
for i=1:length(t2)
gsu2(:,i)=gsu(:,t2(i));
end
for i=1:length(t3)
gsu3(:,i)=gsu(:,t3(i));
end
for i=1:length(t4)
gsu4(:,i)=gsu(:,t4(i));
end
gsu12=[gsu1,gsu2];
gsu34=[gsu3,gsu4];
for gt=gsu12
gpos = gt(1:2);
weight=gt(3);
jac = gt(4);
v = domain(gpos,x,dm,numnod);
L = length(v);
en = zeros(1,2*L);
force=zeros(1,2*L);
[phi,dphix,dphiy] = shape(gpos,dmax,x,v,dm);
tx=-Q1*0.1997*abs(xzuo(2,1)-xzuo(2,2))/(weight*jcob);
ty=0;
for i=1:L
en(2*i-1) = 2*v(i)-1;
en(2*i) = 2*v(i);
force(2*i-1)=tx*phi(i);
force(2*i) = ty*phi(i);
end
f(en) = f(en) + jac*weight*force';
end
for gt=gsu34
gpos = gt(1:2);
weight=gt(3);
jac = gt(4);
v = domain(gpos,x,dm,numnod);
L = length(v);
en = zeros(1,2*L);
force=zeros(1,2*L);
[phi,dphix,dphiy] = shape(gpos,dmax,x,v,dm);
tx=-Q1*0.3003*abs(xzuo(2,1)-xzuo(2,2))/(weight*jcob);
ty=0;
for i=1:L
en(2*i-1) = 2*v(i)-1;
en(2*i) = 2*v(i);
force(2*i-1)=tx*phi(i);
force(2*i) = ty*phi(i);
end
f(en) = f(en) + jac*weight*force';
end
%shijia sanjiao hezai
gstsh=gsu;
ind=0;
for i=1:length(gstsh)
ind=ind+1;
gstshang(:,ind)=gstsh(:,length(gstsh)+1-ind);
end
%INTEGRATE FORCES ALONG BOUNDARY
t1=1:4:length(gstshang);
t2=4:4:length(gstshang);
t3=2:4:length(gstshang);
t4=3:4:length(gstshang);
degao=zeros(3,length(gstshang));
degao(1,t1)=0.1997*ones(1,length(t1))*abs(xyou(2,1)-xyou(2,2));
degao(1,t2)=0.1997*ones(1,length(t1))*abs(xyou(2,1)-xyou(2,2));
degao(1,t3)=0.3003*ones(1,length(t1))*abs(xyou(2,1)-xyou(2,2));
degao(1,t4)=0.3003*ones(1,length(t1))*abs(xyou(2,1)-xyou(2,2));
ind=0;
suan=0;
for i=1:length(gstshang)
ind=ind+1;
suan=suan+degao(1,i);
jisuan(ind)=suan;
end
ind=0;
for i=1:length(gstshang)
ind=ind+1;
suan=Q*jisuan(i)/D;
degao(2,ind)=suan;
end
ind=0;
suanmian(1)=0;
suanmian(2:length(gstshang)+1)=degao(2,:);
for i=1:length(gstshang)
ind=ind+1;
suan=(suanmian(i)+suanmian(i+1))*degao(1,i)/2;
degao(3,ind)=suan;
end
ind=0;
for gt=gstshang
ind=ind+1;
gpos = gt(1:2);
weight=gt(3);
jac = gt(4);
v = domain(gpos,x,dm,numnod);
L = length(v);
en = zeros(1,2*L);
force=zeros(1,2*L);
[phi,dphix,dphiy] = shape(gpos,dmax,x,v,dm);
tx=degao(3,ind)/(weight*jcob);
ty=0;
for i=1:L
en(2*i-1) = 2*v(i)-1;
en(2*i) = 2*v(i);
force(2*i-1)=tx*phi(i);
force(2*i) = ty*phi(i);
end
f(en) = f(en) + jac*weight*force';
end
%SET UP GAUSS POINTS ALONG TRACTION BOUNDARY
ind=0;
gauss=gauss2(quado);
gaussang1=gauss(1,:);
gaussang1=(-1)*gauss(1,:);
gaussang=[gaussang1;gauss(2,:)];
for i=1:(lxs-1)
ycen=(xshang(1,i+1)+xshang(1,i))/2;
jcob = abs((xshang(1,i+1)-xshang(1,i))/2);
for j=1:quado
mark(j) = ycen-gaussang(1,j)*jcob;
ind = ind+1;
gsth(1,ind)=mark(j);
gsth(2,ind)=xshang(2,i);
gsth(3,ind)=gaussang(2,j);
gsth(4,ind)=jcob;
end
end
gsuh=gsth;
gsuh(2,:)=-1*gsth(2,:);
%INTEGRATE FORCES ALONG BOUNDARY
t1=1:4:4*(lxs-1);
t2=4:4:4*(lxs-1);
t3=2:4:4*(lxs-1);
t4=3:4:4*(lxs-1);
gst1=zeros(4,length(t1));
gst2=zeros(4,length(t2));
gst3=zeros(4,length(t3));
gst4=zeros(4,length(t4));
for i=1:length(t1)
gst1(:,i)=gsth(:,t1(i));
end
for i=1:length(t2)
gst2(:,i)=gsth(:,t2(i));
end
for i=1:length(t3)
gst3(:,i)=gsth(:,t3(i));
end
for i=1:length(t4)
gst4(:,i)=gsth(:,t4(i));
end
gst12=[gst1,gst2];
gst34=[gst3,gst4];
for gt=gst12
gpos = gt(1:2);
weight=gt(3);
jac = gt(4);
v = domain(gpos,x,dm,numnod);
L = length(v);
en = zeros(1,2*L);
force=zeros(1,2*L);
[phi,dphix,dphiy] = shape(gpos,dmax,x,v,dm);
tx=0;
ty=-Q2*0.1997*abs(xshang(1,1)-xshang(1,2))/(weight*jcob);
for i=1:L
en(2*i-1) = 2*v(i)-1;
en(2*i) = 2*v(i);
force(2*i-1)=tx*phi(i);
force(2*i) = ty*phi(i);
end
f(en) = f(en) + jac*weight*force';
end
for gt=gst34
gpos = gt(1:2);
weight=gt(3);
jac = gt(4);
v = domain(gpos,x,dm,numnod);
L = length(v);
en = zeros(1,2*L);
force=zeros(1,2*L);
[phi,dphix,dphiy] = shape(gpos,dmax,x,v,dm);
tx=0;
ty=-Q2*0.3003*abs(xshang(1,1)-xshang(1,2))/(weight*jcob);
for i=1:L
en(2*i-1) = 2*v(i)-1;
en(2*i) = 2*v(i);
force(2*i-1)=tx*phi(i);
force(2*i) = ty*phi(i);
end
f(en) = f(en) + jac*weight*force';
end
% INTEGRATE G MATRIX AND Q VECTOR ALONG DISPLACEMENT BOUNDARY
qk = zeros(1,2*lxx);
GG = zeros(numnod*2,lxx*2);
ind1=0;ind2=0;
for i=1:(lxx-1)
ind1=ind1+1;
m1 = ind1; m2 = m1+1;
x1 = xxia(1,m1); x2 =xxia(1,m2);
len = x1-x2;
for j=1:quado
ind2=ind2+1;
gpos = gsuh(1:2,ind2);
weight = gsuh(3,j);
jac = gsuh(4,j);
xp1 = gpos(1,1);
yp1 = gpos(2,1);
uyex1 = 0;
uxex1=0;
N1 = (gpos(1,1)-x2)/len; N2 = 1-N1;
qk(2*m1-1) = qk(2*m1-1)-weight*jac*N1*uxex1;
qk(2*m1) = qk(2*m1) - weight*jac*N1*uyex1;
qk(2*m2-1) = qk(2*m2-1) -weight*jac*N2*uxex1;
qk(2*m2) = qk(2*m2) - weight*jac*N2*uyex1;
v = domain(gpos,x,dm,numnod);
[phi,dphix,dphiy] = shape(gpos,dmax,x,v,dm);
L = length(v);
for n=1:L
G1 = -weight*jac*phi(n)*[N1 0;0 N1];
G2 = -weight*jac*phi(n)*[N2 0;0 N2];
c1=2*v(n)-1;c2=2*v(n);c3=2*m1-1;c4=2*m1;
c5=2*m2-1;c6=2*m2;
GG(c1:c2,c3:c4)=GG(c1:c2,c3:c4)+ G1;
GG(c1:c2,c5:c6)=GG(c1:c2,c5:c6)+G2;
end
end
end
% ENFORCE BC'S USING LAGRANGE MULTIPLIERS
f = [f;zeros(lxx*2,1)];
f(numnod*2+1:numnod*2+2*lxx,1) = qk';
m = sparse([k GG;GG' zeros(lxx*2)]);
d=m\f;
u=d(1:2*numnod);
for i=1:numnod
u2(1,i) = u(2*i-1);
u2(2,i) = u(2*i);
end
% SOLVE FOR OUTPUT VARIABLES - DISPLACEMENTS
displ=zeros(1,2*numnod);
ind=0;
for gg=x
ind = ind+1;
gpos = gg(1:2);
v = domain(gpos,x,dm,numnod);
[phi,dphix,dphiy] = shape(gpos,dmax,x,v,dm);
displ(2*ind-1) = phi*u2(1,v)';
displ(2*ind) = phi*u2(2,v)';
end
for i=1:numnod
disp2(1,i)=displ(2*i-1);
disp2(2,i)=displ(2*i);
end
% SOLVE FOR STRESSES AT GAUSS POINTS
ind = 0;
enorm=0;
for gg=gs
ind = ind+1;
gpos = gg(1:2);
weight = gg(3);
jac = gg(4);
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
stress(1:3,ind) = Dmat*Bmat*u(en);
end
%求解节点应力
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;
disp2;
%绘制开洞墙体应力云图
figure(1)
x1=-Lb/2:0.25:-1.2;
x2=linspace(D/2,-D/2,ndivw+1);
z12=stress1(2,1:102);
z12=reshape(z12,ndivw+1,length(x1));
x3=1.25:0.25:2.5;
x4=linspace(D/2,-D/2,ndivw+1);
z34=stress1(2,202:303);
z34=reshape(z34,ndivw+1,length(x3));
z5678=stress1(2,103:201);
z5678=reshape(z5678,11,9);
x5=-1.25:0.25:1.25;
x6=2:-0.25:1.5;
z56=z5678(1:3,:);
z56=[stress1(2,86:88)' z56 stress1(2,202:204)'];
x7=-1.25:0.25:1.25;
x8=-0.25:-0.25:-2;
z78=z5678(4:11,:);
z78=[stress1(2,95:102)' z78 stress1(2,211:218)'];
pcolor(x1,x2,z12);
hold on
pcolor(x3,x4,z34);
pcolor(x5,x6,z56);
pcolor(x7,x8,z78);
shading interp;
xlabel('开洞墙体各点X坐标');
ylabel('开洞墙体各点Y坐标');
title('开洞墙体各点σy应力云图');
axis equal;
grid on;
hold off
%
figure(2)
x1=-Lb/2:0.25:-1.2;
x2=linspace(D/2,-D/2,ndivw+1);
z12=stress1(1,1:102);
z12=reshape(z12,ndivw+1,length(x1));
x3=1.25:0.25:2.5;
x4=linspace(D/2,-D/2,ndivw+1);
z34=stress1(1,202:303);
z34=reshape(z34,ndivw+1,length(x3));
z5678=stress1(1,103:201);
z5678=reshape(z5678,11,9);
x5=-1.25:0.25:1.25;
x6=2:-0.25:1.5;
z56=z5678(1:3,:);
z56=[stress1(1,86:88)' z56 stress1(1,202:204)'];
x7=-1.25:0.25:1.25;
x8=-0.25:-0.25:-2;
z78=z5678(4:11,:);
z78=[stress1(1,95:102)' z78 stress1(1,211:218)'];
pcolor(x1,x2,z12);
hold on
pcolor(x3,x4,z34);
pcolor(x5,x6,z56);
pcolor(x7,x8,z78);
shading interp;
xlabel('开洞墙体各点X坐标');
ylabel('开洞墙体各点Y坐标');
title('开洞墙体各点σx应力云图');
axis equal;
grid on;
hold off
figure(4)
xb=[-2.5 2.5 2.5 -2.5 -2.5];
yb=[-2 -2 2 2 -2];
plot(xb,yb,'k','LineWidth',1.5);
hold on
xbkk=[-1.25 1.25 1.25 -1.25 -1.25];
ybkk=[-0.25 -0.25 1.5 1.5 -0.25];
plot(xbkk,ybkk,'k','LineWidth',1.5);
axis([-4,4,-3,3]);
xb1=x(:,no1)+10000*disp2(:,no1);
xb2=x(:,no2)+10000*disp2(:,no2);
xb3=x(:,no3)+10000*disp2(:,no3);
xb4=x(:,no4)+10000*disp2(:,no4);
plot(xb1(1,:),xb1(2,:),'r',xb2(1,:),xb2(2,:),'r',xb3(1,:),xb3(2,:),'r',xb4(1,:),xb4(2,:),'r','LineWidth',1.5);
indk3=0;indk4=0;
for j=1:numnod
if x(2,j)==1.5
indk3=indk3+1;
nok3(indk3)=j;
end
if x(2,j)==-0.25
indk4=indk4+1;
nok4(indk4)=j;
end
end
xbkk1=x(:,88:95)+10000*disp2(:,88:95);
xbkk2=x(:,204:211)+10000*disp2(:,204:211);
xbkk3=x(:,nok3(6:16))+10000*disp2(:,nok3(6:16));
xbkk4=x(:,nok4(6:16))+10000*disp2(:,nok4(6:16));
plot(xbkk1(1,:),xbkk1(2,:),'r',xbkk2(1,:),xbkk2(2,:),'r',xbkk3(1,:),xbkk3(2,:),'r',xbkk4(1,:),xbkk4(2,:),'r','LineWidth',1.5);
legend('\fontsize{10}\bf变形前','\fontsize{10}\bf变形后');
title('\fontsize{12}\bf开洞墙体变形图')