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

1 Ansicht (letzte 30 Tage)
% 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

Antworten (0)

Kategorien

Mehr zu 数学 finden Sie in Help Center und File Exchange

Produkte

Community Treasure Hunt

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

Start Hunting!