I have written a code for pressure calculation in bearing. i have found this error. Index in position 1 exceeds array bounds (must not exceed 142). Error in (line 680)

7 views (last 30 days)
%*****************************Input Data***********************************
B= 1; % value of L/D ratio or lamda
NX= 890; % Number of element in x-direction
NY= 141; % Number of element in y-direction
TL= 0.0001; % tolerance
AMU= 1; % dimensionless viscosity
AKY= 12; % Turbulence coefficient
OMEGAR= 1; % Angular velocity
SOMinput= 0.1; % Sommerfeld number
WOinput= 11.4740; % External non-dimensional load
SH= 270; % Angle at which load is applied with respect to x-axis
ISW = 0; % If ISW is zero, Sommerfeld number to be read
% If ISW is one, External non-dimensional load to be read
TOA = 0; % If TOA is 0, Analysis of Plain Journal Bearing
% If TOA is 1, Analysis of Bearing with Circular Texture
% If TOA is 2, Analysis of Bearing with Square Texture
% If TOA is 3, Analysis of Bearing with Densely Distributed Square Texture
alphaT1 = 290; % Start of texture in alpha direction in degree
alphaT2 = 330; % End of texture in alpha direction in degree
beta1 = -0.5; % Start of texture in beta direction
beta2 = 0.5; % End of texture in beta direction
deltah = 1.0; % Depth of texture
R = 0.0286; % Dimensionless Radius of Circular Texture
DBTH = 0.0134; % Dimensionless Distance between Circular texture in Horizontal direction
DBTV = 0.0138; % Dimensionless Distance between Circular texture in Vertical direction
LSTH = 0.0706; % Dimensionless Length of Square Texture in Horizontal direction
LSTV = 0.0710; % Dimensionless Length of Square Texture in Vertical direction
PI= 3.14159265; % value of pi
DVXJ= 0; % disturbance in XJ
DVZJ= 0; % disturbance in ZJ
XOJ = 0.5; % Arbitrary x-coordinate of journal centre
ZOJ = -0.5; % Arbitrary Z-coordinate of journal centre
ITER=100; % Number of iteration
n = 1; % Number of Gauss point
%*****************Data calculation on the basis of input data**************
dalpha = (2.0*PI)/NX;
dbeta = (2*B)/NY;
SHI = (SH*PI)/180; % converting angle into radian
NG = n*n; % Total number of Gauss points inside elements
if (ISW==0)
SOM = SOMinput;
WO = (2.0*B)/(PI*SOM); % calculating load on the basis of sommerfeld number
WOX = WO*cos(SHI); % component of load in horizontal direction
WOZ = WO*sin(SHI); % component of load in vertical direction
elseif (ISW==1)
WO = WOinput; % calculating load on the basis of sommerfeld number
WOX = WO*cos(SHI); % component of load in horizontal direction
WOZ = WO*sin(SHI); % component of load in vertical direction
SOM = (2.0*B)/(PI*WO);
end
alpha1 = (alphaT1*PI)/180; % converting texture start angle into radian
alpha2 = (alphaT2*PI)/180; % converting texture end angle into radian
HDC = DBTH+(2*R); % Horizontal Distance between centres in circular texture
VDC = DBTV+(2*R); % Vertical Distance between centre in circular texture
CX1 = alpha1+(HDC/2); % First centre of cicular texture
CY1 = beta1+(VDC/2);
if (TOA==1)
NTH = round((alpha2-alpha1)/HDC); % Number of texture in horizontal direction
NTV = round((beta2-beta1)/VDC); % Number of texture in vertical direction
elseif (TOA==2)||(TOA==3)
NTH = round((alpha2-alpha1)/LSTH);
NTV = round((beta2-beta1)/LSTV);
end
%***********************Meshing and grading of elements********************
% Meshing the bearing surface into elements
N=NX*NY; % calling is done as M(Node, X=1 or Y=2, Element number)
M = zeros (4,2,N);
for J=1:NX
E=J;
strtdb=-B;
for K=1:NY
if(mod(K,2)==0)
IE=1+((J-1)*2);
end
if(mod(K,2)~=0)
IE=((2*NX)-1)-((J-1)*2);
end
strtdalpha=(J-1)*dalpha;
M(1,1,E)=strtdalpha;
M(1,2,E)=strtdb;
M(2,1,E)=strtdalpha+dalpha;
M(2,2,E)=strtdb;
M(3,1,E)=strtdalpha+dalpha;
M(3,2,E)=strtdb+dbeta;
M(4,1,E)=strtdalpha;
M(4,2,E)=strtdb+dbeta;
E=E+IE;
strtdb=strtdb+dbeta;
end
end
% Meshing the bearing surface into Nodes
T = (NX+1)*(NY+1); % calling is done as Q(X=1 or Y=2, Node number)
Q = zeros(2,T);
for J=1:NX+1
E=J;
strtdb=-B;
for K=1:NY+1
if(mod(K,2)==0)
IE=1+((J-1)*2);
end
if(mod(K,2)~=0)
IE=((2*(NX+1))-1)-((J-1)*2);
end
strtdalpha=(J-1)*dalpha;
Q(1,E)=strtdalpha;
Q(2,E)=strtdb;
E=E+IE;
strtdb=strtdb+dbeta;
end
end
%*****Divide into matrix of Texture elements and Non-texture elements******
% Circular Texture(CT)
if (TOA==1)
% Finding all the centres of textures
NT = NTH*NTV;
C = zeros(2,NT);
for J=1:NTH
E=J;
strtdb=CY1;
for K=1:NTV
if(mod(K,2)==0)
IE=1+((J-1)*2);
end
if(mod(K,2)~=0)
IE=((2*NTH)-1)-((J-1)*2);
end
strtdalpha=CX1+(J-1)*HDC;
C(1,E)=strtdalpha;
C(2,E)=strtdb;
E=E+IE;
strtdb=strtdb+VDC;
end
end
% Dividing elements into matrix of Texture and Non Texture
i=1;
j=1;
for I=1:N
% x-coordinates of nodes of an elements
x1=M(1,1,I);
x2=M(2,1,I);
x3=M(3,1,I);
x4=M(4,1,I);
% y-coordinates of nodes of an elements
y1=M(1,2,I);
y2=M(2,2,I);
y3=M(3,2,I);
y4=M(4,2,I);
if (x1>=alpha1)&&(x2<=alpha2)&&(y1>=beta1)&&(y2<=beta2)
R1= sqrt(((C(1,:)-x1).^2)+((C(2,:)-y1).^2));
if (any(R1<=R))
% x-coordinates of nodes of an elements
MTE(1,1,i)=x1; % MTE(Matrix of Texture elements)
MTE(2,1,i)=x2;
MTE(3,1,i)=x3;
MTE(4,1,i)=x4;
% y-coordinates of nodes of an elements
MTE(1,2,i)=y1;
MTE(2,2,i)=y2;
MTE(3,2,i)=y3;
MTE(4,2,i)=y4;
i=i+1;
else
% x-coordinates of nodes of an elements
MNTE(1,1,j)=x1; % MNTE(Matrix of Non Texture elements)
MNTE(2,1,j)=x2;
MNTE(3,1,j)=x3;
MNTE(4,1,j)=x4;
% y-coordinates of nodes of an elements
MNTE(1,2,j)=y1;
MNTE(2,2,j)=y2;
MNTE(3,2,j)=y3;
MNTE(4,2,j)=y4;
j=j+1;
end
else
% x-coordinates of nodes of an elements
MNTE(1,1,j)=x1; % MNTE(Matrix of Non Texture elements)
MNTE(2,1,j)=x2;
MNTE(3,1,j)=x3;
MNTE(4,1,j)=x4;
% y-coordinates of nodes of an elements
MNTE(1,2,j)=y1;
MNTE(2,2,j)=y2;
MNTE(3,2,j)=y3;
MNTE(4,2,j)=y4;
j=j+1;
end
end
% Dividing Nodes into matrix of Texture and Non Texture
i=1;
j=1;
for I=1:T
x=Q(1,I);
y=Q(2,I);
if (x>=alpha1)&&(x<=alpha2)&&(y>=beta1)&&(y<=beta2)
R1= sqrt(((C(1,:)-x).^2)+((C(2,:)-y).^2));
if (any(R1<=R))
MTN(1,i) = x; % MTN(Matrix of Texture Nodes)
MTN(2,i) = y;
i=i+1;
else
MNTN(1,j) = x; % MNTN(Matrix of Non Texture Nodes)
MNTN(2,j) = y;
j=j+1;
end
else
MNTN(1,j) = x; % MNTN(Matrix of Non Texture Nodes)
MNTN(2,j) = y;
j=j+1;
end
end
% Square Texture(ST)
elseif (TOA==2)
% Dividing elements into matrix of Texture and Non Texture
i=1;
j=1;
for I=1:N
% x-coordinates of nodes of an elements
x1=M(1,1,I);
x2=M(2,1,I);
x3=M(3,1,I);
x4=M(4,1,I);
% y-coordinates of nodes of an elements
y1=M(1,2,I);
y2=M(2,2,I);
y3=M(3,2,I);
y4=M(4,2,I);
if (x1>=alpha1)&&(x2<=alpha2)&&(y1>=beta1)&&(y2<=beta2)
V = zeros(1,NTH);
q=1;
for r=1:NTH
if (mod(r,2)~=0)
a=alpha1+(r-1)*LSTH;
b=alpha1+(r)*LSTH;
if (x1>=a)&&(x2<=b)
V(1,q)=1;
q=q+1;
else
V(1,q)=0;
q=q+1;
end
else
V(1,q)=0;
q=q+1;
end
end
V1 = zeros(1,NTV);
q=1;
for r=1:NTV
if (mod(r,2)~=0)
a=beta1+(r-1)*LSTV;
b=beta1+(r)*LSTV;
if (y1>=a)&&(y2<=b)
V1(1,q)=1;
q=q+1;
else
V1(1,q)=0;
q=q+1;
end
else
V1(1,q)=0;
q=q+1;
end
end
if (any(V))&&(any(V1))
% x-coordinates of nodes of an elements
MTE(1,1,i)=x1; % MTE(Matrix of Texture elements)
MTE(2,1,i)=x2;
MTE(3,1,i)=x3;
MTE(4,1,i)=x4;
% y-coordinates of nodes of an elements
MTE(1,2,i)=y1;
MTE(2,2,i)=y2;
MTE(3,2,i)=y3;
MTE(4,2,i)=y4;
i=i+1;
else
% x-coordinates of nodes of an elements
MNTE(1,1,j)=x1; % MNTE(Matrix of Non Texture elements)
MNTE(2,1,j)=x2;
MNTE(3,1,j)=x3;
MNTE(4,1,j)=x4;
% y-coordinates of nodes of an elements
MNTE(1,2,j)=y1;
MNTE(2,2,j)=y2;
MNTE(3,2,j)=y3;
MNTE(4,2,j)=y4;
j=j+1;
end
else
% x-coordinates of nodes of an elements
MNTE(1,1,j)=x1; % MNTE(Matrix of Non Texture elements)
MNTE(2,1,j)=x2;
MNTE(3,1,j)=x3;
MNTE(4,1,j)=x4;
% y-coordinates of nodes of an elements
MNTE(1,2,j)=y1;
MNTE(2,2,j)=y2;
MNTE(3,2,j)=y3;
MNTE(4,2,j)=y4;
j=j+1;
end
end
% Dividing Nodes into matrix of Texture and Non Texture
i=1;
j=1;
for I=1:T
x=Q(1,I);
y=Q(2,I);
if (x>=alpha1)&&(x<=alpha2)&&(y>=beta1)&&(y<=beta2)
q=1;
for r=1:NTH
if (mod(r,2)~=0)
a=alpha1+(r-1)*LSTH;
b=alpha1+(r)*LSTH;
if (x>=a)&&(x<=b)
V(1,q)=1;
q=q+1;
else
V(1,q)=0;
q=q+1;
end
else
V(1,q)=0;
q=q+1;
end
end
q=1;
for r=1:NTV
if (mod(r,2)~=0)
a=beta1+(r-1)*LSTV;
b=beta1+(r)*LSTV;
if (y>=a)&&(y<=b)
V1(1,q)=1;
q=q+1;
else
V1(1,q)=0;
q=q+1;
end
else
V1(1,q)=0;
q=q+1;
end
end
if (any(V))&&(any(V1))
MTN(1,i) = x; % MTN(Matrix of Texture Nodes)
MTN(2,i) = y;
i=i+1;
else
MNTN(1,j) = x; % MNTN(Matrix of Non Texture Nodes)
MNTN(2,j) = y;
j=j+1;
end
else
MNTN(1,j) = x; % MNTN(Matrix of Non Texture Nodes)
MNTN(2,j) = y;
j=j+1;
end
end
% Densely Distributed Square Texture (DDST)
elseif (TOA==3)
% Dividing elements into matrix of Texture and Non Texture
i=1;
j=1;
for I=1:N
% x-coordinates of nodes of an elements
x1=M(1,1,I);
x2=M(2,1,I);
x3=M(3,1,I);
x4=M(4,1,I);
% y-coordinates of nodes of an elements
y1=M(1,2,I);
y2=M(2,2,I);
y3=M(3,2,I);
y4=M(4,2,I);
if (x1>=alpha1)&&(x2<=alpha2)&&(y1>=beta1)&&(y2<=beta2)
V = zeros(1,NTH);
q=1;
for r=1:NTH
if (mod(r,2)~=0)
a=alpha1+(r-1)*LSTH;
b=alpha1+(r)*LSTH;
if (x1>=a)&&(x2<=b)
V(1,q)=1;
q=q+1;
else
V(1,q)=0;
q=q+1;
end
else
V(1,q)=0;
q=q+1;
end
end
V1 = zeros(1,NTV);
q=1;
for r=1:NTV
if (mod(r,2)~=0)
a=beta1+(r-1)*LSTV;
b=beta1+(r)*LSTV;
if (y1>=a)&&(y2<=b)
V1(1,q)=1;
q=q+1;
else
V1(1,q)=0;
q=q+1;
end
else
V1(1,q)=0;
q=q+1;
end
end
V2 = zeros(1,NTH);
q=1;
for r=1:NTH
if (mod(r,2)==0)
a=alpha1+(r-1)*LSTH;
b=alpha1+(r)*LSTH;
if (x1>=a)&&(x2<=b)
V2(1,q)=1;
q=q+1;
else
V2(1,q)=0;
q=q+1;
end
else
V2(1,q)=0;
q=q+1;
end
end
V3 = zeros(1,NTV);
q=1;
for r=1:NTV
if (mod(r,2)==0)
a=beta1+(r-1)*LSTV;
b=beta1+(r)*LSTV;
if (y1>=a)&&(y2<=b)
V3(1,q)=1;
q=q+1;
else
V3(1,q)=0;
q=q+1;
end
else
V3(1,q)=0;
q=q+1;
end
end
if (any(V))&&(any(V1))
% x-coordinates of nodes of an elements
MTE(1,1,i)=x1; % MTE(Matrix of Texture elements)
MTE(2,1,i)=x2;
MTE(3,1,i)=x3;
MTE(4,1,i)=x4;
% y-coordinates of nodes of an elements
MTE(1,2,i)=y1;
MTE(2,2,i)=y2;
MTE(3,2,i)=y3;
MTE(4,2,i)=y4;
i=i+1;
elseif (any(V2))&&(any(V3))
% x-coordinates of nodes of an elements
MTE(1,1,i)=x1; % MTE(Matrix of Texture elements)
MTE(2,1,i)=x2;
MTE(3,1,i)=x3;
MTE(4,1,i)=x4;
% y-coordinates of nodes of an elements
MTE(1,2,i)=y1;
MTE(2,2,i)=y2;
MTE(3,2,i)=y3;
MTE(4,2,i)=y4;
i=i+1;
else
% x-coordinates of nodes of an elements
MNTE(1,1,j)=x1; % MNTE(Matrix of Non Texture elements)
MNTE(2,1,j)=x2;
MNTE(3,1,j)=x3;
MNTE(4,1,j)=x4;
% y-coordinates of nodes of an elements
MNTE(1,2,j)=y1;
MNTE(2,2,j)=y2;
MNTE(3,2,j)=y3;
MNTE(4,2,j)=y4;
j=j+1;
end
else
% x-coordinates of nodes of an elements
MNTE(1,1,j)=x1; % MNTE(Matrix of Non Texture elements)
MNTE(2,1,j)=x2;
MNTE(3,1,j)=x3;
MNTE(4,1,j)=x4;
% y-coordinates of nodes of an elements
MNTE(1,2,j)=y1;
MNTE(2,2,j)=y2;
MNTE(3,2,j)=y3;
MNTE(4,2,j)=y4;
j=j+1;
end
end
% Dividing Nodes into matrix of Texture and Non Texture
i=1;
j=1;
for I=1:T
x=Q(1,I);
y=Q(2,I);
if (x>=alpha1)&&(x<=alpha2)&&(y>=beta1)&&(y<=beta2)
q=1;
for r=1:NTH
if (mod(r,2)~=0)
a=alpha1+(r-1)*LSTH;
b=alpha1+(r)*LSTH;
if (x>=a)&&(x<=b)
V(1,q)=1;
q=q+1;
else
V(1,q)=0;
q=q+1;
end
else
V(1,q)=0;
q=q+1;
end
end
q=1;
for r=1:NTV
if (mod(r,2)~=0)
a=beta1+(r-1)*LSTV;
b=beta1+(r)*LSTV;
if (y>=a)&&(y<=b)
V1(1,q)=1;
q=q+1;
else
V1(1,q)=0;
q=q+1;
end
else
V1(1,q)=0;
q=q+1;
end
end
q=1;
for r=1:NTH
if (mod(r,2)==0)
a=alpha1+(r-1)*LSTH;
b=alpha1+(r)*LSTH;
if (x>=a)&&(x<=b)
V2(1,q)=1;
q=q+1;
else
V2(1,q)=0;
q=q+1;
end
else
V2(1,q)=0;
q=q+1;
end
end
q=1;
for r=1:NTV
if (mod(r,2)==0)
a=beta1+(r-1)*LSTV;
b=beta1+(r)*LSTV;
if (y>=a)&&(y<=b)
V3(1,q)=1;
q=q+1;
else
V3(1,q)=0;
q=q+1;
end
else
V3(1,q)=0;
q=q+1;
end
end
if (any(V))&&(any(V1))
MTN(1,i) = x; % MTN(Matrix of Texture Nodes)
MTN(2,i) = y;
i=i+1;
elseif (any(V2))&&(any(V3))
MTN(1,i) = x; % MTN(Matrix of Texture Nodes)
MTN(2,i) = y;
i=i+1;
else
MNTN(1,j) = x; % MNTN(Matrix of Non Texture Nodes)
MNTN(2,j) = y;
j=j+1;
end
else
MNTN(1,j) = x; % MNTN(Matrix of Non Texture Nodes)
MNTN(2,j) = y;
j=j+1;
end
end
end
%***********************Mapping of element*********************************
Mmap = zeros(NG,2,N);
DANGL = zeros(NG,N);
for I=1:N % I is number of element called
%calling of elements with its nodal coordinates
% x-coordinates of nodes of an element
x1=M(1,1,I);
x2=M(2,1,I);
x3=M(3,1,I);
x4=M(4,1,I);
% y-coordinates of nodes of an element
y1=M(1,2,I);
y2=M(2,2,I);
y3=M(3,2,I);
y4=M(4,2,I);
% Array of x and y coordinates and number of gauss points
x=[x1,x2,x3,x4];
y=[y1,y2,y3,y4];
[MP,DL]=gaussintmap(x,y,n);
for q = 1:NG
Mmap(q,1,I)=MP(1,q);
Mmap(q,2,I)=MP(2,q);
DANGL(q,I)=DL(q);
end
end
%****************Starting Analysis of Various Bearing Surface**************
% Alpha Coordinate Matrix
alphaP = zeros(1,NX+1);
for i=1:NX+1
alphaP(i)= dalpha*(i-1);
end
% Beta Coordinate Matrix
betaP = zeros(1,NY+1);
for j=1:NY+1
betaP(j)=-B+(j-1)*dbeta;
end
for q=1:ITER
Iteration=q
Centre=[XOJ ZOJ];
%************************Calculation of Pressure***************************
pbar = zeros(1,T);
for I=1:T
x=Q(1,I);
y=Q(2,I);
if (TOA==0)||(ismember(y,MNTN(2,x==MNTN(1,:))))
hbar(I)=1.0-(XOJ*cos(x))-(ZOJ*sin(x));
else
hbar(I)=1.0-(XOJ*cos(x))-(ZOJ*sin(x))+deltah;
end
end
% Converting all negative pressure to zeros
pbar(pbar<0)=0;
% Arranging Pressure values in matrix form as like plate meshing
P = zeros(NY+1,NX+1);
for J=1:NX+1
E=J;
for K=1:NY+1
if(mod(K,2)==0)
IE=1+((J-1)*2);
end
if(mod(K,2)~=0)
IE=((2*(NX+1))-1)-((J-1)*2);
end
hbar(K,J) = hbar(E);
P(K,J)= (P(I+1,J) + P(I-1,J) + ((((dalpha)/(dbeta))^2)*(P(I,J+1)+ P(I,J-1)))+ (((3*((XOJ*sin(dalpha))- (ZOJ*cos(dalpha))))*((P(I+1,J)- P(I-1,J)))*(dlpha))/(2*hbar))-((((XOJ*sin(dalpha))-(ZOJ*cos(dalpha)))*((dalpha)^2))/((hbar)^3)))/(2*(1+((dalpha)/(dbata))^2));
E=E+IE;
end
end
P(P<0.000001)=0;
%Starting and end of positive pressure
Pos1 = (alphaP(find(any(P),1,'first')));
Pos2 = (alphaP(find(any(P),1,'last')));
Pos3 = (alphaP(find(all(P==0),1,'first')));
Pos4 = (alphaP(find(all(P==0),1,'last')));
if (Pos1==0)&&(Pos2==(2*PI))
alphaPos1=Pos1;
alphaPos2=Pos3;
alphaPos3=Pos4;
alphaPos4=Pos2;
else
alphaPos1=Pos1;
alphaPos2=Pos2;
alphaPos3=Pos1;
alphaPos4=Pos2;
end
end
Iteration = 1
Index in position 1 exceeds array bounds. Index must not exceed 142.
function [MP,DL] = gaussintmap(x,y,n)
if (n==1)
POSGP(1) = 0.0;
WEIGP(1) = 2.0;
elseif (n==2)
POSGP(1) = -0.577350269189626;
POSGP(2) = -POSGP(1);
WEIGP(1) = 1.0000000000000000;
WEIGP(2) = WEIGP(1);
elseif (n==3)
POSGP(1) = -0.7745966692414833;
POSGP(2) = 0.0;
POSGP(3) = -POSGP(1);
WEIGP(1) = 0.5555555555555556;
WEIGP(2) = 0.8888888888888889;
WEIGP(3) = WEIGP(1);
end
nt=n*n;
MP = zeros(2,nt);
DL = zeros(nt);
for i=1:n
E=i;
for j=1:n
if(mod(j,2)==0)
IE=1+((i-1)*2);
end
if(mod(j,2)~=0)
IE=((2*n)-1)-((i-1)*2);
end
N(1)=0.25*(1-POSGP(i))*(1-POSGP(j));
N(2)=0.25*(1+POSGP(i))*(1-POSGP(j));
N(3)=0.25*(1+POSGP(i))*(1+POSGP(j));
N(4)=0.25*(1-POSGP(i))*(1+POSGP(j));
MP(1,E) = dot(N,x); %xnode
MP(2,E) = dot(N,y); %ynode
J(1,1) = - (1 - POSGP(j))*x(1) + (1 - POSGP(j))*x(2) + (1 + POSGP(j))*x(3) - (1 + POSGP(j))*x(4);
J(1,2) = - (1 - POSGP(j))*y(1) + (1 - POSGP(j))*y(2) + (1 + POSGP(j))*y(3) - (1 + POSGP(j))*y(4);
J(2,1) = - (1 - POSGP(i))*x(1) - (1 + POSGP(i))*x(2) + (1 + POSGP(i))*x(3) + (1 - POSGP(i))*x(4);
J(2,2) = - (1 - POSGP(i))*y(1) - (1 + POSGP(i))*y(2) + (1 + POSGP(i))*y(3) + (1 - POSGP(i))*y(4);
detJ = (J(1,1)*J(2,2) - J(1,2)*J(2,1))/16;
DL(E) = WEIGP(i)*WEIGP(j)*detJ;
E=E+IE;
end
end
end
Note - My error is this, kindly help
Index in position 1 exceeds array bounds (must not exceed 142).
Error in shivamodi (line 680)
P(K,J)= (P(I+1,J) + P(I-1,J) + ((((dalpha)/(dbeta))^2)*(P(I,J+1)+ P(I,J-1)))+ (((3*((XOJ*sin(dalpha))-
(ZOJ*cos(dalpha))))*((P(I+1,J)-
P(I-1,J)))*(dlpha))/(2*hbar))-((((XOJ*sin(dalpha))-(ZOJ*cos(dalpha)))*((dalpha)^2))/((hbar)^3)))/(2*(1+((dalpha)/(dbata))^2));
Pressure is being calculated for hydrodynamic journal bearing with surface texture. With above given data.
  7 Comments
RAHUL KUMAR
RAHUL KUMAR on 22 Jan 2022
function [MP,DL] = gaussintmap(x,y,n)
if (n==1)
POSGP(1) = 0.0;
WEIGP(1) = 2.0;
elseif (n==2)
POSGP(1) = -0.577350269189626;
POSGP(2) = -POSGP(1);
WEIGP(1) = 1.0000000000000000;
WEIGP(2) = WEIGP(1);
elseif (n==3)
POSGP(1) = -0.7745966692414833;
POSGP(2) = 0.0;
POSGP(3) = -POSGP(1);
WEIGP(1) = 0.5555555555555556;
WEIGP(2) = 0.8888888888888889;
WEIGP(3) = WEIGP(1);
end
nt=n*n;
MP = zeros(2,nt);
DL = zeros(nt);
for i=1:n
E=i;
for j=1:n
if(mod(j,2)==0)
IE=1+((i-1)*2);
end
if(mod(j,2)~=0)
IE=((2*n)-1)-((i-1)*2);
end
N(1)=0.25*(1-POSGP(i))*(1-POSGP(j));
N(2)=0.25*(1+POSGP(i))*(1-POSGP(j));
N(3)=0.25*(1+POSGP(i))*(1+POSGP(j));
N(4)=0.25*(1-POSGP(i))*(1+POSGP(j));
MP(1,E) = dot(N,x); %xnode
MP(2,E) = dot(N,y); %ynode
J(1,1) = - (1 - POSGP(j))*x(1) + (1 - POSGP(j))*x(2) + (1 + POSGP(j))*x(3) - (1 + POSGP(j))*x(4);
J(1,2) = - (1 - POSGP(j))*y(1) + (1 - POSGP(j))*y(2) + (1 + POSGP(j))*y(3) - (1 + POSGP(j))*y(4);
J(2,1) = - (1 - POSGP(i))*x(1) - (1 + POSGP(i))*x(2) + (1 + POSGP(i))*x(3) + (1 - POSGP(i))*x(4);
J(2,2) = - (1 - POSGP(i))*y(1) - (1 + POSGP(i))*y(2) + (1 + POSGP(i))*y(3) + (1 - POSGP(i))*y(4);
detJ = (J(1,1)*J(2,2) - J(1,2)*J(2,1))/16;
DL(E) = WEIGP(i)*WEIGP(j)*detJ;
E=E+IE;
end
end
end

Sign in to comment.

Answers (1)

Walter Roberson
Walter Roberson on 22 Jan 2022
T = (NX+1)*(NY+1); % calling is done as Q(X=1 or Y=2, Node number)
for I=1:T
x=Q(1,I);
y=Q(2,I);
if (TOA==0)||(ismember(y,MNTN(2,x==MNTN(1,:))))
hbar(I)=1.0-(XOJ*cos(x))-(ZOJ*sin(x));
else
hbar(I)=1.0-(XOJ*cos(x))-(ZOJ*sin(x))+deltah;
end
end
After that loop, I will be left as the last value it was assigned -- so I will be left the same as T, which is (NX+1)*(NY+1)
P = zeros(NY+1,NX+1);
for J=1:NX+1
E=J;
for K=1:NY+1
if(mod(K,2)==0)
IE=1+((J-1)*2);
end
if(mod(K,2)~=0)
IE=((2*(NX+1))-1)-((J-1)*2);
end
hbar(K,J) = hbar(E);
P(K,J)= (P(I+1,J) + P(I-1,J) + ((((dalpha)/(dbeta))^2)*(P(I,J+1)+ P(I,J-1)))+ (((3*((XOJ*sin(dalpha))- (ZOJ*cos(dalpha))))*((P(I+1,J)- P(I-1,J)))*(dlpha))/(2*hbar))-((((XOJ*sin(dalpha))-(ZOJ*cos(dalpha)))*((dalpha)^2))/((hbar)^3)))/(2*(1+((dalpha)/(dbata))^2));
E=E+IE;
end
end
You refer to P(I+1,J) . But I = T = (NX+1)*(NY+1) and P(I+1,J) would be P((NX+1)*(NY+1)+1,J) which is a problem because P is only size NY+1 in the first dimension.
It is suspicious that you are looping over K but you are not indexing P by values that are computed from K.
  2 Comments
Walter Roberson
Walter Roberson on 23 Jan 2022
My guess would be that everywhere in the line
P(K,J)= (P(I+1,J) + P(I-1,J) + ((((dalpha)/(dbeta))^2)*(P(I,J+1)+ P(I,J-1)))+ (((3*((XOJ*sin(dalpha))- (ZOJ*cos(dalpha))))*((P(I+1,J)- P(I-1,J)))*(dlpha))/(2*hbar))-((((XOJ*sin(dalpha))-(ZOJ*cos(dalpha)))*((dalpha)^2))/((hbar)^3)))/(2*(1+((dalpha)/(dbata))^2));
that you refer to I that you should instead be referring to K

Sign in to comment.

Categories

Find more on MATLAB in Help Center and File Exchange

Products


Release

R2020a

Community Treasure Hunt

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

Start Hunting!