bridge max load error--Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.204862e-18.

1 view (last 30 days)
% Units: Kg & mm
m=21; % Number of elements
n=12; % Number of nodes
Coord(1,1:2)=[0 0];
Coord(2,1:2)=[27.5 0];
Coord(3,1:2)=[27.5+87.5 0];
Coord(4,1:2)=[27.5+87.5*2 0];
Coord(5,1:2)=[27.5+87.5*3 0];
Coord(6,1:2)=[27.5+87.5*4 0];
Coord(7,1:2)=[27.5*2+87.5*4 0];
Coord(8,1:2)=[27.5+87.5*4 87.5];
Coord(9,1:2)=[27.5+87.5*3 87.5];
Coord(10,1:2)=[27.5+87.5*2 87.5];
Coord(11,1:2)=[27.5+87.5 87.5];
Coord(12,1:2)=[27.5 87.5];
Coord(:,3)=0;
Con=[1 2;1 12;2 3;2 12;3 4;3 12;3 11;3 10;4 5;4 10;5 6;5 10;5 9;5 8;6 7;6 8;7 8;8 9;9 10;10 11;11 12];
Con(:,3:4)=0;
Re=ones(n,6);
Re(:,1:2)=0;
Re(1,2)=1;
Re(7,2)=1;
Load=zeros(n,6);Load(10,2)=-200;
w=zeros(m,3);
E=ones(1,m)*2.1e6;nu=0.3;G=E/(2*(1+nu));
A=ones(1,m)*100;Iz=ones(1,m)*1000;Iy=ones(1,m)*1000;J=ones(1,m)*50;
St=zeros(n,6);
be=zeros(1,m);
D=struct('m',m,'n',n,'Coord',Coord','Con',Con','Re',Re','Load',Load',...
'w',w','E',E','G',G','A',A','Iz',Iz','Iy',Iy','J',J','St',St','be',be'); % A transpose is done at this point
[Q,V,R]=MSA(D);
for mm=1:m
TrussForce(D,Q,mm)
end
MSAPlot2(D,Q,V,R)
function [Q,V,R]=MSA(D)
m=D.m;n=D.n;Ni=zeros(12,12,m);S=zeros(6*n);Pf=S(:,1);Q=zeros(12,m);Qfi=Q;Ei=Q;
for i=1:m
H=D.Con(:,i);C=D.Coord(:,H(2))-D.Coord(:,H(1));e=[6*H(1)-5:6*H(1),6*H(2)-5:6*H(2)];c=D.be(i);
[a,b,L]=cart2sph(C(1),C(3),C(2));ca=cos(a);sa=sin(a);cb=cos(b);sb=sin(b);cc=cos(c);sc=sin(c);
r=[1 0 0;0 cc sc;0 -sc cc]*[cb sb 0;-sb cb 0;0 0 1]*[ca 0 sa;0 1 0;-sa 0 ca];T=kron(eye(4),r);
co=2*L*[6/L 3 2*L L];x=D.A(i)*L^2;y=D.Iy(i)*co;z=D.Iz(i)*co;g=D.G(i)*D.J(i)*L^2/D.E(i);
K1=diag([x,z(1),y(1)]);K2=[0 0 0;0 0 z(2);0 -y(2) 0];K3=diag([g,y(3),z(3)]);K4=diag([-g,y(4),z(4)]);
K=D.E(i)/L^3*[K1 K2 -K1 K2;K2' K3 -K2' K4;-K1 -K2 K1 -K2;K2' K4 -K2' K3];
w=D.w(:,i)';Qf=-L^2/12*[6*w/L 0 -w(3) w(2) 6*w/L 0 w(3) -w(2)]';Qfs=K*T*D.St(e)';
A=diag([0 -0.5 -0.5]);B(2,3)=1.5/L;B(3,2)=-1.5/L;W=diag([1,0,0]);Z=zeros(3);M=eye(12);p=4:6;q=10:12;
switch 2*H(3)+H(4)
case 0;B=2*B/3;M(:,[p,q])=[-B -B;W Z;B B;Z W];case 1;M(:,p)=[-B;W;B;A];case 2;M(:,q)=[-B;A;B;W];
end
K=M*K;Ni(:,:,i)=K*T;S(e,e)=S(e,e)+T'*Ni(:,:,i);Qfi(:,i)=M*Qf;Pf(e)=Pf(e)+T'*M*(Qf+Qfs);Ei(:,i)=e;
end
V=1-(D.Re|D.St);f=find(V);V(f)=S(f,f)\(D.Load(f)-Pf(f));R=reshape(S*V(:)+Pf,6,n);R(f)=0;V=V+D.St;
for i=1:m
Q(:,i)=Ni(:,:,i)*V(Ei(:,i))+Qfi(:,i);
end

Answers (0)

Categories

Find more on 线性代数 in Help Center and File Exchange

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!