finite element method with matrix error

조회 수: 5 (최근 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
이전 댓글 13개 표시
Walter Roberson
Walter Roberson 2020년 10월 29일
sz = size(k);
if sz(1) < n*d
k(sz(1)+1:n*d, :) = k(sz(1),:) + (1:n*d-sz(1));
end
if sz(2) < n*d
k(:, sz(2)+1:n*d) = (1:n*d-sz(2)).' + k(:, sz(2));
end
You only talked about extending according to the last value on each "line", but your matrices are short on both sides, so this code also extends by adding 1 along columns as needed.
I think you are talking the wrong approach: I think you should be initializing to the maximum size.
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?

카테고리

Help Center 및 File Exchange에서 Logical에 대해 자세히 알아보기

태그

Community Treasure Hunt

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

Start Hunting!

Translated by