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)

Question

RAHUL KUMAR 2022-1-20

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1632945-i-have-written-a-code-for-pressure-calculation-in-bearing-i-have-found-this-error-index-in-positi

评论： Walter Roberson 2022-1-23

gaussintmap.m

%*****************************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 个评论
显示 5更早的评论隐藏 5更早的评论

Rik 2022-1-22

You have edited your question, but you haven't attached the gaussintmap function.

RAHUL KUMAR 2022-1-22

在 MATLAB Online 中打开

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

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Walter Roberson 2022-1-22

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1632945-i-have-written-a-code-for-pressure-calculation-in-bearing-i-have-found-this-error-index-in-positi#answer_879725

在 MATLAB Online 中打开

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 个评论
显示无隐藏无

RAHUL KUMAR 2022-1-23

So what should I do next?will yo please suggest me.

Walter Roberson 2022-1-23

在 MATLAB Online 中打开

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

请先登录，再进行评论。