Can someone rearrange the code to run
3 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
MINATI
am 25 Jun. 2020
Kommentiert: MINATI
am 25 Jun. 2020
%subdivisions space /time
ht=0.01; Tmax=1.2; nx=33; hx=1/(nx-1); x=[0:hx:1]';
%matrices
K=stiff(1/pi^2,hx,nx); M=mass(1/ht,hx,nx); A=M+K;
%boundary conditions by deleting rows
A(nx,:)=[]; A(:,nx)=[]; A(1,:)=[]; A(:,1)=[];
%creation
R=chol(A);
%initial condition
u=cos(pi*(x-1/2));
%time step loop
k=0;
while (k*ht<Tmax)
k=k+1;
b=M*u;
b(nx)=[];b(1)=[];
%solve
u=zeros(nx,1);
u(2:nx-1)=R\(R'\b);
end
function R=stiff(nu,h,n)
%
[m1,m2]=size(nu);
if (length(nu)==1)
ee=nu*ones(n-1,1)/h; e1=[ee;0]; e2=[0;ee]; e=e1+e2;
else
ee=.5*(nu(1:n-1)+nu(2:n))/h; e1=[ee;0]; e2=[0;ee]; e=e1+e2;
end
R=spdiags([-e1 e -e2],-1:1,n,n);
return
function M=mass(alpha,h,n)
[m1,m2]=size(alpha);
if(length(alpha)==1)
ee=alpha*h*ones(n-1,1)/6;
e1=[ee;0]; e2=[0;ee]; e=e1+e2;
else
ee=h*(alpha(1:n-1)+alpha(2:n))/12;
e1=[ee;0]; e2=[0;ee]; e=e1+e2;
end
M=spdiags([e1 2*e e2],-1:1,n,n);
return
% Set the matrices to zero
k=zeros(6,6,2); K=zeros(6,6,2); Gamma=zeros(6,6,2);
% Enter parameter values, in units: mm^2, mm^4, and MPa(10^6 N/m)
A=6000; II=200*10^6; EE=200000;
% Convert units into meters and kN
A=A/10^6; II = II/10^12; EE =EE*1000;
% Element 1
i=[0,0]; j=[7.416,3];
[k(:,:,1),K(:,:,1),Gamma(:,:,1)]=stiff(EE,II,A,i,j);
% Element 2
i=j; j=[15.416,3];
[k(:,:,2),K(:,:,2),Gamma(:,:,2)]=stiff(EE,II,A,i,j);
% Define the elementary stiffness matrix
ID=[-4 1 -7;-5 2 -8;-6 3 -9];
% Define the connection matrix
LM=[-4 -5 -6 1 2 3; 1 2 3 -7 -8 -9];
% Assemble the augmented stiffness matrix
Kaug = zeros(9);
for elem=1:2
for r=1:6
lr=abs(LM(elem,r));
for c=1:6
lc=abs(LM(elem,c));
Kaug(lr,lc)=Kaug(lr,lc)+K(r,c,elem);
end
end
end
% Extract the stiffness matrix
Ktt=Kaug(1:3,1:3);
% Determine the reactions at the nodes in local coordinates
fea(1:6,1)=0;
fea(1:6,2)=[0,8*4/2,4*8^2/12,0,8*4/2,-4*8^2/12]';
% Determine the reactions in the global coordinate system
FEA(1:6,1)=Gamma(:,:,1)'*fea(1:6,1);
FEA(1:6,2)=Gamma(:,:,2)'*fea(1:6,2);
% FEA_Rest for constrained nodes
FEA_Rest=[0,0,0,FEA(4:6,2)'];
% Assemble the right-hand side for non-constrained nodes
P(1)=50*3/8; P(2)=-50*7.416/8-FEA(2,2); P(3)=-FEA(3,2);
% Solve to find the displacements in meters and in radians
Displacements=inv(Ktt)*P'
% Extract Kut
Kut = Kaug(4:9,1:3);
% Compute the reactions and introduce boundary conditions
Reactions=Kut*Displacements+FEA_Rest'
% Solve to find the internal forces, excluding fixed points
dis_global(:,:,1)=[0,0,0,Displacements(1:3)'];
dis_global(:,:,2)=[Displacements(1:3)',0,0,0];
for elem=1:2
dis_local= Gamma(:,:,elem)*dis_global(:,:,elem)';
int_forces= k(:,:,elem)*dis_local+fea(1:6,elem)
end
%The above script calls the stiff function, which can be implemented as follows:
function [k,K,Gamma] = stiff( EE,II,A,i,j )
% Find the length
L=sqrt((j(2)-i(2))^2+(j(1)-i(1))^2);
% Compute the angle theta
if(j(1)-i(1))~=2
alpha=atan((j(2)-i(2))/(j(1)-i(1)))
else
alpha=-pi/2;
end
% Form the rotation matrix Gamma
Gamma =[cos(alpha) sin(alpha) 0 0 0 0;
-sin(alpha) cos(alpha) 0 0 0 0;
0 0 1 0 0 0;
0 0 0 cos(alpha) sin(alpha) 0;
0 0 0 -sin(alpha) cos(alpha) 0;
0 0 0 0 0 1];
% Form the elementary stiffness matrix in local coordinates
EI=EE*II; EA=EE*A;
k=[EA/L, 0, 0, -EA/L, 0, 0;
0, 12*EI/L^3, 6*EI/L^2, 0, -12*EI/L^3,6*EI/L^2;
0, 6*EI/L^2, 4*EI/L, 0 -6*EI/L^2, 2*EI/L;
-EA/L, 0 ,0 , EA/L, 0, 0;
0, -12*EI/L^3, -6*EI/L^2, 0, 12*EI/L^3, -6*EI/L^2;
0, 6*EI/L^2, 2*EI/L, 0, -6*EI/L^2, 4*EI/L];
% Elementary matrix in global coordinates
K=Gamma'*k*Gamma;
end
end
end
%%%% I got this code from a book to solve heat eqn but unable to arrange this.
%% Please someone arrange this so that it could be of use
2 Kommentare
Akzeptierte Antwort
Rafael Hernandez-Walls
am 25 Jun. 2020
%subdivisions space /time
ht=0.01; Tmax=1.2; nx=33; hx=1/(nx-1); x=[0:hx:1]';
%matrices
K=stiff2(1/pi^2,hx,nx); M=mass(1/ht,hx,nx);
A=M+K;
%boundary conditions by deleting rows
A(nx,:)=[]; A(:,nx)=[]; A(1,:)=[]; A(:,1)=[];
%creation
R=chol(A);
%initial condition
u=cos(pi*(x-1/2));
%time step loop
k=0;
while (k*ht<Tmax)
k=k+1;
b=M*u;
b(nx)=[];b(1)=[];
%solve
u=zeros(nx,1);
u(2:nx-1)=R\(R'\b);
end
% Set the matrices to zero
k=zeros(6,6,2); K=zeros(6,6,2); Gamma=zeros(6,6,2);
% Enter parameter values, in units: mm^2, mm^4, and MPa(10^6 N/m)
A=6000; II=200*10^6; EE=200000;
% Convert units into meters and kN
A=A/10^6; II = II/10^12; EE =EE*1000;
% Element 1
i=[0,0]; j=[7.416,3];
[k(:,:,1),K(:,:,1),Gamma(:,:,1)]=stiff(EE,II,A,i,j);
% Element 2
i=j; j=[15.416,3];
[k(:,:,2),K(:,:,2),Gamma(:,:,2)]=stiff(EE,II,A,i,j);
% Define the elementary stiffness matrix
ID=[-4 1 -7;-5 2 -8;-6 3 -9];
% Define the connection matrix
LM=[-4 -5 -6 1 2 3; 1 2 3 -7 -8 -9];
% Assemble the augmented stiffness matrix
Kaug = zeros(9);
for elem=1:2
for r=1:6
lr=abs(LM(elem,r));
for c=1:6
lc=abs(LM(elem,c));
Kaug(lr,lc)=Kaug(lr,lc)+K(r,c,elem);
end
end
end
% Extract the stiffness matrix
Ktt=Kaug(1:3,1:3);
% Determine the reactions at the nodes in local coordinates
fea(1:6,1)=0;
fea(1:6,2)=[0,8*4/2,4*8^2/12,0,8*4/2,-4*8^2/12]';
% Determine the reactions in the global coordinate system
FEA(1:6,1)=Gamma(:,:,1)'*fea(1:6,1);
FEA(1:6,2)=Gamma(:,:,2)'*fea(1:6,2);
% FEA_Rest for constrained nodes
FEA_Rest=[0,0,0,FEA(4:6,2)'];
% Assemble the right-hand side for non-constrained nodes
P(1)=50*3/8; P(2)=-50*7.416/8-FEA(2,2); P(3)=-FEA(3,2);
% Solve to find the displacements in meters and in radians
Displacements=inv(Ktt)*P'
% Extract Kut
Kut = Kaug(4:9,1:3);
% Compute the reactions and introduce boundary conditions
Reactions=Kut*Displacements+FEA_Rest'
% Solve to find the internal forces, excluding fixed points
dis_global(:,:,1)=[0,0,0,Displacements(1:3)'];
dis_global(:,:,2)=[Displacements(1:3)',0,0,0];
for elem=1:2
dis_local= Gamma(:,:,elem)*dis_global(:,:,elem)';
int_forces= k(:,:,elem)*dis_local+fea(1:6,elem)
end
%The above script calls the stiff function, which can be implemented as follows:
function [k,K,Gamma] = stiff( EE,II,A,i,j )
% Find the length
L=sqrt((j(2)-i(2))^2+(j(1)-i(1))^2);
% Compute the angle theta
if(j(1)-i(1))~=2
alpha=atan((j(2)-i(2))/(j(1)-i(1)))
else
alpha=-pi/2;
end
% Form the rotation matrix Gamma
Gamma =[cos(alpha) sin(alpha) 0 0 0 0;
-sin(alpha) cos(alpha) 0 0 0 0;
0 0 1 0 0 0;
0 0 0 cos(alpha) sin(alpha) 0;
0 0 0 -sin(alpha) cos(alpha) 0;
0 0 0 0 0 1];
% Form the elementary stiffness matrix in local coordinates
EI=EE*II; EA=EE*A;
k=[EA/L, 0, 0, -EA/L, 0, 0;
0, 12*EI/L^3, 6*EI/L^2, 0, -12*EI/L^3,6*EI/L^2;
0, 6*EI/L^2, 4*EI/L, 0 -6*EI/L^2, 2*EI/L;
-EA/L, 0 ,0 , EA/L, 0, 0;
0, -12*EI/L^3, -6*EI/L^2, 0, 12*EI/L^3, -6*EI/L^2;
0, 6*EI/L^2, 2*EI/L, 0, -6*EI/L^2, 4*EI/L];
% Elementary matrix in global coordinates
K=Gamma'*k*Gamma;
end
function R=stiff2(nu,h,n)
%
[m1,m2]=size(nu);
if (length(nu)==1)
ee=nu*ones(n-1,1)/h; e1=[ee;0]; e2=[0;ee]; e=e1+e2;
else
ee=.5*(nu(1:n-1)+nu(2:n))/h; e1=[ee;0]; e2=[0;ee]; e=e1+e2;
end
R=spdiags([-e1 e -e2],-1:1,n,n);
return
end
function M=mass(alpha,h,n)
[m1,m2]=size(alpha);
if(length(alpha)==1)
ee=alpha*h*ones(n-1,1)/6;
e1=[ee;0]; e2=[0;ee]; e=e1+e2;
else
ee=h*(alpha(1:n-1)+alpha(2:n))/12;
e1=[ee;0]; e2=[0;ee]; e=e1+e2;
end
M=spdiags([e1 2*e e2],-1:1,n,n);
return
end
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Custom Geometry and PCB Fabrication finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!