Can someone rearrange the code to run

조회 수: 3 (최근 30일)
MINATI
MINATI 2020년 6월 25일
댓글: MINATI 2020년 6월 25일
%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
Rafael Hernandez-Walls
Rafael Hernandez-Walls 2020년 6월 25일
You just have to define the function mass
MINATI
MINATI 2020년 6월 25일
Dear Rafael
Like stiff, mass is also defined in line 31 (I guess) or, please modify if possible

댓글을 달려면 로그인하십시오.

채택된 답변

Rafael Hernandez-Walls
Rafael Hernandez-Walls 2020년 6월 25일
%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
  댓글 수: 1
MINATI
MINATI 2020년 6월 25일
Dear Rafael
it worked
Thank you

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Custom Geometry and PCB Fabrication에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by