计算力学小论文
《网架 matlab 编程》
姓 名：
学 号：
专 业：土木工程
比例：
说明：
本程序为适应正放四角锥、正放抽放四角锥、棋
盘星四角锥三种类型的网架，可根据输入不同的参数
进行变化。前处理比较简单，运用方便。
本程序节点编号采用横向进行，先上层再下层。
这三种都可以运行，并已经附了真实的算例。
最后尝试了非线性的求解，但是结果不收敛，但
在 这 里 附 上 程 序 。
主程序
E=input('输入弹性模量（单位为 Pa）:');
A=input('输入各个单元截面面积（单位为平方米）:');
m=input('输入横向网格个数:');
n=input('输入纵向网格个数:');
Lx=input('输入横向网格尺寸:');
Ly=input('输入纵向网格尺寸:');
h=input('输入网架厚度:');
p=input('向下的节点力大小:');
flag=input('网架的边界约束，0 代表铰接；1 代表刚接:');
t1=input('第一个节点编号:');
t2=input('第二个节点编号:');
M=input('1:正放四角锥；2：正放抽空四角锥；3：棋盘星四角锥:');
if M==1
%定位坐标
for i=1:m+1
for j=1:n+1
x(i+(m+1)*(j-1))=(i-1)*Lx;
y(i+(m+1)*(j-1))=(j-1)*Ly;
z(i+(m+1)*(j-1))=h;
end
end
for i=1:m
for j=1:n
x((m+1)*(n+1)+i+m*(j-1))=Lx/2+(i-1)*Lx;
y((m+1)*(n+1)+i+m*(j-1))=Ly/2+(j-1)*Ly;
z((m+1)*(n+1)+i+m*(j-1))=0;
end
end
%建立总刚矩阵
N=((m+1)*(n+1)+m*n)*3;
u=zeros(N,1);
while e<=10^-6
for i=1:N/3
x(i)=x(i)+u([i],1);
y(i)=y(i)+u([i],1);
z(i)=z(i)+u([i],1);
end
KK=zeros(N,N);
for j=1:n+1
for i=1:m
a=i+(m+1)*(j-1);
b=i+(m+1)*(j-1)+1;
k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b));
KK=SpaceTrussAssemble(KK,k,a,b);
end
end
for j=1:n
for i=1:m
a=i+(m+1)*(j-1);
b=i+(m+1)*j;
k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b));
KK=SpaceTrussAssemble(KK,k,a,b);
a=i+(m+1)*(j-1);
b=(m+1)*(n+1)+i+m*(j-1);
k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b));
KK=SpaceTrussAssemble(KK,k,a,b);
a=i+1+(m+1)*(j-1);
b=(m+1)*(n+1)+i+m*(j-1);
k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b));
KK=SpaceTrussAssemble(KK,k,a,b);
a=i+(m+1)*j;
b=(m+1)*(n+1)+i+m*(j-1);
k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b));
KK=SpaceTrussAssemble(KK,k,a,b);
a=i+1+(m+1)*j;
b=(m+1)*(n+1)+i+m*(j-1);
k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b));
KK=SpaceTrussAssemble(KK,k,a,b);
end
end
for j=1:n
for i=1:m-1
a=(m+1)*(n+1)+i+m*(j-1);
b=(m+1)*(n+1)+i+m*(j-1)+1;
k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b));
KK=SpaceTrussAssemble(KK,k,a,b);
end
end
for j=1:n-1
for i=1:m
a=(m+1)*(n+1)+i+m*(j-1);
b=(m+1)*(n+1)+i+m*j;
k=SpaceTrussElementStiffness(E,A,x(a),y(a),z(a),x(b),y(b),z(b));
KK=SpaceTrussAssemble(KK,k,a,b);
end
K=KK;
end
F=K*u;
P=zeros(((m+1)*(n+1)+m*n)*3,1);
for i=1:(m+1)*(n+1)
P(3*i-1,1)=-p;
end
P=P-F;
Restrain=zeros(1,((m+1)*(n+1)+m*n)*3);
if flag==1
Restrain(1,(m+1)*(n+1)*3+1:((m+1)*(n+1)+m)*3)=1;
Restrain(1,((m+1)*(n+1)+m*n-m)*3+1:((m+1)*(n+1)+m*n)*3)=1;
for j=1:n
Restrain(1,(1+(m+1)*(n+1)+m*(j-1))*3-2:(1+(m+1)*(n+1)+m*(j-1))*3)=1;
Restrain(1,(m+(m+1)*(n+1)+m*(j-1))*3-2:(m+(m+1)*(n+1)+m*(j-1))*3)=1;
end
end
if flag==0
for i=1:m
Restrain(1,((m+1)*(n+1)+i)*3-2:((m+1)*(n+1)+i)*3-1)=1;
Restrain(1,((m+1)*(n+1)+m*n-m+i)*3-2:((m+1)*(n+1)+m*n-m+i)*3-1)=1;
end
for j=1:n
Restrain(1,(1+(m+1)*(n+1)+m*(j-1))*3-2:(1+(m+1)*(n+1)+m*(j-1))*3-1)=1;
Restrain(1,(m+(m+1)*(n+1)+m*(j-1))*3-2:(m+(m+1)*(n+1)+m*(j-1))*3-1)=1;
end
end
end
if M==2
%定位坐标
for i=1:m+1
for j=1:n+1
x(i+(m+1)*(j-1))=(i-1)*Lx;
y(i+(m+1)*(j-1))=(j-1)*Ly;
z(i+(m+1)*(j-1))=h;
end
end
m1=(m-1)/2+1;
n1=(n-1)/2+1;
m2=(m+1)*(n+1);
for j=1:n1
for i=1:m
x(m2+i+(m+m1)*(j-1))=Lx/2+(i-1)*Lx;
y(m2+i+(m+m1)*(j-1))=Ly/2+(j-1)*2*Ly;
z(m2+i+(m+m1)*(j-1))=0;
