contact_deformation.m

main.m

S.m

SUBAK.m

untitled.jpg

最近又做了個案例，是關于點接觸彈性變形數值計算，具體描述參考《彈性流體動壓潤滑數值計算方法》P23-29，黃平著。以下為一個算例的matlab實現方法(原書為FORTRAN語言)，與大家交流，請大家多提意見。

• 主程序

clear

clc

N=33;

global AK

AK=SUBAK(N);

KL=log(N-1)/log(2)-1.99;

DX=2.4/(N-1);

for I=1:N

    X(I)=-1.2+DX*(I-1);

    A=X(I)^2;

    for J=1:N

        Y(J)=-1.2+DX*(J-1);

        P(I,J)=0;

        H(I,J)=0.5*A+0.5*Y(J)^2;

    end

end

M=0;

for I=1:N

    for J=1:N

        A=1.0-X(I)^2-Y(J)^2;

        if A>=0

            P(I,J)=sqrt(A);

        end

    end

end

V=contact_deformation(N,DX,P,0);

for I=1:N

    for J=1:N

        H(I,J)=H(I,J)+V(I,J);

    end

end

%% painting

figure(1)

subplot(2,2,1)

mesh(X,Y,H)

title('')

subplot(2,2,2)

mesh(X,Y,V)

subplot(2,2,3)

mesh(X,Y,P)

subplot(2,2,4)

mesh(X,Y,(H+V))

• 調用程序1

function S=S(X,Y)

S=X+sqrt(X^2+Y^2);

• 調用程序2

function AK=SUBAK(MM)

AK(1:65,1:65)=0;

for I=1:MM

    XP=I+0.5;

    XM=I-0.5;

    for J=1:MM

        YP=J+0.5;

        YM=J-0.5;

        A1=S(YP,XP)/S(YM,XP);

        A2=S(XM,YM)/S(XP,YM);

        A3=S(YM,XM)/S(YP,XM);

        A4=S(XP,YP)/S(XM,YP);

        AK(I,J)=XP*log(A1)+YM*log(A2)+XM*log(A3)+YP*log(A4);

        AK(J,I)=AK(I,J);

    end

end

• 調用程序3

function V=contact_deformation(N,DX,P,V)

PAI1=0.2026423;

global AK

for I=1:N

    for J=1:N

        H0=0;

        for K=1:N

            IK=abs(I-K);

            for L=1:N

                JL=abs(J-L);

                H0=H0+AK(IK+1,JL+1)*P(K,L);

                V(I,J)=H0*DX*PAI1;

            end

        end

    end

end