MATLAB Answers

finite element method with matrix error

조회 수: 2(최근 30일)
zina shadidi
zina shadidi 2020년 10월 26일
댓글: zina shadidi 2020년 10월 29일
I am trying to calculate energy using the folowing code , the matrices elements2s, nodes2s was the results of another script.
clc
clear all
% Initial data:
TubeLength = 3.38e-9; %m
TubeDiameter = 1.011e-9; %m
D = 3; %degrees of freedom at each node
presc_disp=1e-9; %m
L=0.141e-9; %m Length of a C-C bond
numIterations=1; %number of iterations
numNaN=0;
% %%%%% the ref .
%%Force constants from Tserpes et al. (2005):
% Kr is EA=6.52e-7 %N·nm^-1
% ktheta is EI=8.76e-10 %N·nm·rad^-2
% ktao is GJ=2.78e-10 %N·nm·rad^-2
EA=L*652 %N
EI=L*8.76e-19 %N·m^2·rad^-2
GJ=L*2.78e-19 %N·m^2·rad^-2
% Transformation to nondimensional space:
% [F]=EI/L^2; [L]=L. Then, multiply the resulting forces with
% EI/L^2 and displacements with L to get forces in MN and
% displacements in nm.
Lreal=L;
EIreal=EI;
presc_disp_real=presc_disp;
EA=EA/EI*L^2;
GJ=GJ/EI;
presc_disp=presc_disp/L;
L=1;
EI=1;
% Stiffness matrix of a C-C bond (adapted to 3 DOF):
KBond=[EA/L 0 0 -EA/L 0 0
0 12*EI/L^3 0 0 -12*EI/L^3 0 ;
0 0 12*EI/L^3 0 0 -12*EI/L^3;
-EA/L 0 0 EA/L 0 0 ;
0 -12*EI/L^3 0 0 12*EI/L^3 0 ;
0 0 -12*EI/L^3 0 0 12*EI/L^3]
% Calculations:
it=1;
while it<numIterations+1
load elements2s.txt
load nodes2s.txt
nNodes2s= 199;
nElements2s = 199;
nodesDef=zeros(1,nNodes2s);
k=zeros(nNodes2s*D,nNodes2s*D);
% Defects:
defectsPercentage=1; %per cent
numDef=round((nElements2s*defectsPercentage)/100);
for i=1:numDef
def=round((nElements2s-1)*rand)+1; %The addition of 1 and -1 is to avoid zero as the resulting random number. elements(def,:)=[];
nElements2s=nElements2s-1;
end
zeroMatrix = [0 0 0; 0 0 0; 0 0 0];
vi = [1 0 0];
vj = [0 1 0];
vk = [0 0 1];
for i=1:nElements2s-1
n = elements2s(i,1);
m = elements2s(i,2);
localZ = [1-2, 3-2, 3-4];
localZ = localZ/sqrt(localZ*localZ');
% Find local y
if not(localZ(3)==0)
y3 = (-localZ(1)-2*localZ(2))/localZ(3);
localY = [1 2 y3];
elseif not(localZ(2)==0)
y2 = (-localZ(1)-2*localZ(3))/localZ(2);
localY = [1 y2 2];
elseif not(localZ(1)==0)
y1 = (-localZ(2)-2*localZ(3))/localZ(1);
localY = [y1 1 2];
end
localY = localY/sqrt(localY*localY');
% Find local z
localX = cross(localY, localZ);
localX = localX/sqrt(localX*localX');
% Values for rotation matrix
cosZI = localZ*vi';
cosZJ = localZ*vj';
cosZK = localZ*vk';
cosYI = localY*vi';
cosYJ = localY*vj';
cosYK = localY*vk';
cosXI = localX*vi';
cosXJ = localX*vj';
cosXK = localX*vk';
rotationSmall = [
cosZI cosZJ cosZK;
cosYI cosYJ cosYK;
cosXI cosXJ cosXK];
RR = [rotationSmall zeroMatrix;zeroMatrix rotationSmall];
KElement = RR'*KBond*RR;
n = elements2s(i,1);
m = elements2s(i,2);
k((n-1)*D+1:n*D,(n-1)*D+1:n*D)=k((n-1)*D+1:n*D, (n-1)*D+1:n*D)+KElement(1:D, 1:D);
k((m-1)*D+1:m*D,(m-1)*D+1:m*D)=k((m-1)*D+1:m*D, (m-1)*D+1:m*D)+KElement(1:D, 1:D);
k((n-1)*D+1:n*D,(m-1)*D+1:m*D)=k((n-1)*D+1:n*D, (m-1)*D+1:m*D)+KElement(1:D, D+1:D*2);
k((m-1)*D+1:m*D,(n-1)*D+1:n*D)=k((m-1)*D+1:m*D, (n-1)*D+1:n*D)+KElement(D+1:D*2, 1:D);
end
head = [1 4 5 16 17 34 35 51 52 72 73 86 87 476 477 478 479 482 483 488 489 490 491 492 493 494;
1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 ;
1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 ;
1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2
];
fload = [1 4 5 16 17 34 35 51 52 72 73 86 87 476 477 478 479 482 483 488 489 490 491 492 493 494;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ;
0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1
]*presc_disp;
save k_real.mat k;
loads=zeros(1,D*nNodes);
for i=1:size(head(1,:),2)
for j=1:D
if head(j+1,i)~=0 % displacement prescribed
for m=1:D*nNodes
loads(m)=loads(m)-k(m,(head(1,i)-1)*D+j)*fload(j+1,i);
end
k((head(1,i)-1)*D+j,:)=zeros(1,D*nNodes);
k(:,(head(1,i)-1)*D+j)=zeros(D*nNodes,1);
k((head(1,i)-1)*D+j,(head(1,i)-1)*D+j)=1.0;
loads(1,(head(1,i)-1)*D+j)=fload(j+1,i);
if fload(j+1,i) ~= 0
kdof=j;
end
else % force prescribed
loads(1,(head(1,i)-1)*D+j)=fload(j+1,i);
if fload(j+1,i) ~= 0
kdof=j;
end
end
end
end
clear RR KEelement;
clear elements;
disps=loads/k;
clear k;
load k_real.mat
forces=k*disps';
clear k
tip=[476 477 478 479 482 483 488 489 490 491 492 493 494];
% 'tip'contains the node numbers at one tip
avdisp=0;
for i=1:size(tip,2)
tip_disps(i)=disps((tip(i)-1)*D+kdof);
avdisp=avdisp+disps((tip(i)-1)*D+kdof);
end
fixed=[1 4 5 16 17 34 35 51 52 72 73 86 87];
% 'fixed'contains the node numbers at the other tip
totforce=0;
for i=1:size(tip,2)
totforce=totforce+forces((tip(i)-1)*D+kdof);
end
avdisp=avdisp/size(tip,2);
avdisp=avdisp*Lreal; % in meters
tip_disps=tip_disps*Lreal; % in meters
totforce = totforce * EIreal/Lreal^2; % in Newtons
YY=totforce/(TubeDiameter*avdisp/TubeLength*pi) % N/m=Pa*m
Young(it)=YY;
it=it+1;
end
matlab writng an error which is
{Index in position 1 is invalid. Array indices must be positive integers or logical values.
Error in Untitled (line 99)
k((m-1)*dof+1:m*dof,(m-1)*dof+1:m*dof)=k((m-1)*dof+1:m*dof, (m-1)*dof+1:m*dof)+KElement(1:dof, 1:dof);
>> }
can anyone help me please
many thanks in advance
  댓글 수: 15
zina shadidi
zina shadidi 2020년 10월 29일
Walter Roberson, I can understand you but I have Null matrices as a results of young modulus , and other variables , and this is a wrong results in physics . In fact I dont know what I will do , may be I should start with another code.
I do appreciate your help and thank you very much.

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

답변(1개)

Asad (Mehrzad)  Khoddam
Asad (Mehrzad) Khoddam 2020년 10월 26일
Can you show the values of the nodes and elements data?

Community Treasure Hunt

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

Start Hunting!

Translated by