如何將代碼更改爲DYNAMIC LOAD？

Question

我幾乎是MATLAB新手。 我完全不知道如何更改我的代碼以正確執行。 這裏是關於分析一個3d桁架的代碼，並增加了一個特殊點的負載。

function D=DataT3D 
%m number of elements 
%n number of nodes 
m=25;n=10; 
%coordinates of nodes [(X Y Z) for each node] 
Coord=[-37.5 0 200;37.5 0 200;-37.5 37.5 100;37.5 37.5 100;37.5 -37.5 
100;-37.5 -37.5 100; 
-100 100 0;100 100 0;100 -100 0;-100 -100 0]; 
%conection of the nodes [first in coordinates is the first node and ...] 
Con=[1 2;1 4;2 3;1 5;2 6;2 4;2 5;1 3;1 6;3 6;4 5;3 4;5 6;3 10;6 7;4 9;5 8;4 
7;3 8;5 10;6 9; 
6 10;3 7;4 8;5 9]; Con(:,3:4)=0; 
%Re degrees of freedom for each node (free=0 & fixed=1) 
Re=ones(n,6); 
Re(1:6,1:3)=zeros(6,3); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% concentrated loads on nodes 
Load=zeros(n,6); 
Load([1,2,3,6],1:3)=[1 -10 -10;0 -10 -10;0.5 0 0;0.6 0 0]; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% uniform loads in local coordinate system 
w=zeros(m,3); 
% E: material elastic modules G:shear elastic modules J:torsional constant 
E=ones(1,m)*1e4;nu=0.3;G=E/(2*(1+nu)); 
% A:cross sectional area and Iy Iz: moment of inertia 
A=ones(1,m)*0.5;Iz=ones(1,m);Iy=ones(1,m);J=ones(1,m); 
%St: settlement of supports & displacements of free nodes 
St=zeros(n,6); be=zeros(1,m); 
% All of the variables are transposed and stored in a structure array in the 
name of D 
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'); 
end 

此代碼作爲函數運行另一個代碼。

function [Q,V,R]=MSA(D); 
m=D.m; 
n=D.n; 
% the matrix to store K*T for each member 12*12*m 
Ni=zeros(12,12,m); 
% global stiffness matrix of the structure 6n*6n 
S=zeros(6*n); 
% element fixed end forces in global coordinate 6n*1 
Pf=S(:,1); 
% internal forces and moments in local coordinate system for each member 
% 12*m 
Q=zeros(12,m); 
% element fixed end forces in local coordinate for each member 12*m 
Qfi=Q; 
% member code numbers* (mcn) in global stiffness matrix for each member 
% 12*m 
Ei=Q; 
for i=1:m 
% connectivity and release of the both member ends 4*1 
H=D.Con(:,i); 
% difference of beginning and end nodes coordinates 3*1 
C=D.Coord(:,H(2))-D.Coord(:,H(1)); 
% member code numbers (mcn) in global stiffness matrix for a member 
% 1*12 
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]; 
% transformation matrix related to the 
% coordinate transformation which in considering member orientation 
% 12*12 
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); 
% local stiffness matrix for each member 
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]; 
% uniform loads in local coordinate system for each member 1*3 
w=D.w(:,i)'; 
% local fixed-end force vector for a member, corresponding to external 
% loads 12*1 
Qf=-L^2/12*[6*w/L 0 -w(3) w(2) 6*w/L 0 w(3) -w(2)]'; 
% local fixed-end force vector for a member, corresponding to support 
% displacements 12*1 
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; 
% type of member release* 0 1 2 3 
% M: A matrix for modifying stiffness matrix and local fixed-end force 
vector of a released member ends such K=M*K , Qf=M*Qf and Qfs=M*Qfs 
switch 2*H(3)+H(4) 
case 0;B=2*B/3; % released at both ends 
    M(:,[p,q])=[-B -B;W Z;B B;Z W]; 
case 1; % released at the beginning 
    M(:,p)=[-B;W;B;A]; 
case 2; % released at the end 
    M(:,q)=[-B;A;B;W]; 
end 
K=M*K;Ni(:,:,i)=K*T; 
% global stiffness matrix of the structure 6n*6n 
S(e,e)=S(e,e)+T'*Ni(:,:,i); 
Qfi(:,i)=M*Qf; 
% element fixed end forces in global coordinate 6n*1 
Pf(e)=Pf(e)+T'*M*(Qf+Qfs); 
% member code numbers* (mcn) in global stiffness matrix for each member 
% 12*m 
Ei(:,i)=e; 
end 
% Deflections in global coordinate syste 6*n 
V=1-(D.Re|D.St); 
% f: A vector that indicates the number of degree of freedom ndof*1 
f=find(V); 
V(f)=S(f,f)\(D.Load(f)-Pf(f)); 
% Supports reactions in global coordinate system 6*n 
R=reshape(S*V(:)+Pf,6,n); 
R(f)=0;V=V+D.St; 
for i=1:m 
% internal forces and moments in local coordinate system 12*m 
Q(:,i)=Ni(:,:,i)*V(Ei(:,i))+Qfi(:,i); 
end 

因此，這裏的問題 如何負載變化的動態形式？現在代碼是靜態的。它意味着我們現在將一個集中負載放到一個節點上，然後我們得到答案，例如根據所添加的集中負載來移動其他節點。但是我們怎樣才能將它改變成一種動態的加載形式，這意味着我們在一個節點上添加一個具有不同時間值（例如5秒）的負載，並以不同的時間步驟得到其他節點的答案。

可以使用運行第二個代碼：

d = DataT3D; [Q，V，R] = MSA（d）; **

幫助需要。

附代碼： DataT3D MSA

Answer 1

試試這個：

DataT3D

function D=DataT3D 
%m number of elements 
%n number of nodes 
m=25;n=10; 
%coordinates of nodes [(X Y Z) for each node] 
Coord=[-37.5 0 200;37.5 0 200;-37.5 37.5 100;37.5 37.5 100;37.5 -37.5 100;-37.5 -37.5 100; 
-100 100 0;100 100 0;100 -100 0;-100 -100 0]; 
%conection of the nodes [first in coordinates is the first node and ...] 
Con=[1 2;1 4;2 3;1 5;2 6;2 4;2 5;1 3;1 6;3 6;4 5;3 4;5 6;3 10;6 7;4 9;5 8;4 7;3 8;5 10;6 9; 
6 10;3 7;4 8;5 9]; 
Con(:,3:4)=0; 
%Re degrees of freedom for each node (free=0 & fixed=1) 
Re=ones(n,6); 
Re(1:6,1:3)=zeros(6,3); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% concentrated loads on nodes 

record = [1 2 3 4 5 6 7 8 9 10]; % Example forces on node 1 on different time intervals 
% Create storage for Q, V and R 
allQ = cell(2,1); 
allV = cell(2,1); 
allR = cell(2,1); 


for t=1:length(record) 

    Load = [record(t) 0 0 0 0 0; zeros(n*1,6)]; % Load has only a Fx and all other forces and moments are zero 

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
    % uniform loads in local coordinate system 
    w=zeros(m,3); 
    % E: material elastic modules G:shear elastic modules J:torsional constant 
    E=ones(1,m)*1e4;nu=0.3;G=E/(2*(1+nu)); 
    % A:cross sectional area and Iy Iz: moment of inertia 
    A=ones(1,m)*0.5;Iz=ones(1,m);Iy=ones(1,m);J=ones(1,m); 
    %St: settlement of supports & displacements of free nodes 
    St=zeros(n,6); be=zeros(1,m); 
    % All of the variables are transposed and stored in a structure array in the 
    %name of D 
    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'); 

    [allQ{t},allV{t},allR{t}]=MSA(D); % Save the results for Q, V and R  

end 

allQ = cell2mat(allQ) 
allV = cell2mat(allV) 
allR = cell2mat(allR) 

end 

MSA

function [Q,V,R]=MSA(D) 
m=D.m; 
n=D.n; 
% the matrix to store K*T for each member 12*12*m 
Ni=zeros(12,12,m); 
% global stiffness matrix of the structure 6n*6n 
S=zeros(6*n); 
% element fixed end forces in global coordinate 6n*1 
Pf=S(:,1); 
% internal forces and moments in local coordinate system for each member 
% 12*m 
Q=zeros(12,m); 
% element fixed end forces in local coordinate for each member 12*m 
Qfi=Q; 
% member code numbers* (mcn) in global stiffness matrix for each member 
% 12*m 
Ei=Q; 
for i=1:m 
% connectivity and release of the both member ends 4*1 
H=D.Con(:,i); 
% difference of beginning and end nodes coordinates 3*1 
C=D.Coord(:,H(2))-D.Coord(:,H(1)); 
% member code numbers (mcn) in global stiffness matrix for a member 
% 1*12 
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]; 
% transformation matrix related to the 
% coordinate transformation which in considering member orientation 
% 12*12 
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); 
% local stiffness matrix for each member 
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]; 
% uniform loads in local coordinate system for each member 1*3 
w=D.w(:,i)'; 
% local fixed-end force vector for a member, corresponding to external 
% loads 12*1 
Qf=-L^2/12*[6*w/L 0 -w(3) w(2) 6*w/L 0 w(3) -w(2)]'; 
% local fixed-end force vector for a member, corresponding to support 
% displacements 12*1 
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; 
% type of member release* 0 1 2 3 
% M: A matrix for modifying stiffness matrix and local fixed-end force 
%vector of a released member ends such K=M*K , Qf=M*Qf and Qfs=M*Qfs 
switch 2*H(3)+H(4) 
case 0;B=2*B/3; % released at both ends 
    M(:,[p,q])=[-B -B;W Z;B B;Z W]; 
case 1; % released at the beginning 
    M(:,p)=[-B;W;B;A]; 
case 2; % released at the end 
    M(:,q)=[-B;A;B;W]; 
end 
K=M*K;Ni(:,:,i)=K*T; 
% global stiffness matrix of the structure 6n*6n 
S(e,e)=S(e,e)+T'*Ni(:,:,i); 
Qfi(:,i)=M*Qf; 
% element fixed end forces in global coordinate 6n*1 
Pf(e)=Pf(e)+T'*M*(Qf+Qfs); 
% member code numbers* (mcn) in global stiffness matrix for each member 
% 12*m 
Ei(:,i)=e; 
end 
% Deflections in global coordinate syste 6*n 
V=1-(D.Re|D.St); 
% f: A vector that indicates the number of degree of freedom ndof*1 
f=find(V); 
D.Load(f); 
V(f)=S(f,f)\(D.Load(f)-Pf(f)); 
% Supports reactions in global coordinate system 6*n 
R=reshape(S*V(:)+Pf,6,n); 
R(f)=0;V=V+D.St; 
for i=1:m 
% internal forces and moments in local coordinate system 12*m 
Q(:,i)=Ni(:,:,i)*V(Ei(:,i))+Qfi(:,i); 
end