aokengdizai.rar_scaleom7_弹流润滑_接触_润滑_表面凹坑

global vis0 dens0 dX U0 ph b R E w X K z arf G W D S1 S2 vis0=1.25; % Pa.s dens0=870; % kg/m^3 z=0.68; arf=z*5.1e-9*(log(vis0)+9.67); % m^2/N dX=1/80; % Mesh Density E1=2e11;E2=2e11; pois1=0.3;pois2=0.3; E=2/((1-pois1^2)/E1+(1-pois2^2)/E2); % 当量弹性模量 R1=19.5e-3;R2=5e-3;R=R1*R2/(R1+R2); % 当量曲率半径 u0=0.2; U0=u0*vis0/E/R % 速度参数 L=0.018; w=1000/L; b=sqrt(8*w*R/pi/E); ph=sqrt(w*E/2/pi/R); G=arf*E W=w/E/R % 材料参数、载荷参数 K=0.75*pi^2*U0/W^2; % 变形矩阵 X=1:521;X=dX*(X-401); for j=2:2:520 a1(j)=1/dX^2; a2(j)=X(j)/dX^2-(X(j)+X(j+1))/2/dX^2; a3(j)=-2/dX^2;a4(j)=-2*X(j)/dX^2+(X(j-1)+X(j+1))/dX^2; a5(j)=1/dX^2; a6(j)=X(j)/dX^2-(X(j)+X(j-1))/2/dX^2; end dD=zeros(521); D=zeros(521); for i=1:521 for j=2:2:520 Zmin(i,j)=X(i)-X(j-1); Zmax(i,j)=X(i)-X(j+1); if Zmin(i,j)==0 M0(i,j)=-Zmax(i,j)^2*(log(Zmax(i,j)^2)-3); N0(i,j)=Zmax(i,j)^3*(log(abs(Zmax(i,j))^3)-4); else if Zmax(i,j)==0 M0(i,j)=Zmin(i,j)^2*(log(Zmin(i,j)^2)-3); N0(i,j)=-Zmin(i,j)^3*(log(abs(Zmin(i,j))^3)-4); else M0(i,j)=Zmin(i,j)^2*(log(Zmin(i,j)^2)-3)-Zmax(i,j)^2*(log(Zmax(i,j)^2)-3); N0(i,j)=Zmax(i,j)^3*(log(abs(Zmax(i,j))^3)-4)-Zmin(i,j)^3*(log(abs(Zmin(i,j))^3)-4); end end L0(i,j)=M0(i,j)*(X(i)-X(j))/2+2/9*N0(i,j); dD(i,j-1)=dD(i,j-1)-1/2/pi*(a1(j)*L0(i,j)+a2(j)*M0(i,j)/2); dD(i,j)=-1/2/pi*(a3(j)*L0(i,j)+a4(j)*M0(i,j)/2); dD(i,j+1)=dD(i,j+1)-1/2/pi*(a5(j)*L0(i,j)+a6(j)*M0(i,j)/2); end end for i=1:521 for j=1:521 D(i,j)=dD(i,j); end end % 初始压力 P=zeros(1,521); for i=321:481 P(i)=sqrt(1-X(i)^2); end clear i hc=3.07*arf^0.57*R^0.4*(vis0*u0)^0.71/E^0.03/w^0.11;Hc=hc*R/b^2 Xe=1; TP=zeros(521,1); TP(2:520)=P(2:520); TP(1)=Hc; TP(521)=Xe; plot(X(2:520),TP(2:520));hold on %%%%%%%%%%%%%%%%%% N=floor((Xe+5)/dX)+1; Defl=D*P'-0.25*log(R^2*8*W/pi); S1=zeros(1,521);S2=zeros(1,521); for i=1:521 S1(i)=R/b^2*0.04e-6*sin(2*pi/1*X(i)); % 波纹度 end for i=457:473 S2(i)=1.1*sqrt(1-((X(i)-X(465))/0.1)^2); % 凹坑 end for i=1:521 H(i)=Hc+X(i)^2/2+Defl(i)-Defl(401)+S2(i); end for i=1:521 Vis(i)=exp((log(vis0)+9.67)*((1+5.1e-9*ph*P(i))^z-1)); Dens(i)=1+0.6e-9*ph*P(i)/(1+1.7e-9*ph*P(i)); end He=H(N)+(Xe-X(N))/dX*(H(N+1)-H(N)); plot(X,H);hold on dPdX=zeros(1,521); for i=1:N dPdX(i)=(P(i+1)-P(i))/dX; end ddPX=zeros(520,520); for i=1:N for j=i:i+1 ddPX(i,j)=(eq(i+1,j)-eq(i,j))/dX; end end ddPX(N,N+1)=0; clear f g f=H.^3.*dPdX+K*Vis.*(He./Dens-H); clear PD1 PD2 for i=1:521 PD1(i,i)=Vis(i)*5.1e-9*ph*z*(log(vis0)+9.67)*(1+5.1e-9*ph*P(i))^(z-1); PD2(i,i)=-0.6e-9*ph/(1+2.3e-9*ph*P(i))^2; end clear i clear dfdP for i=1:520 for j=1:520 dHdP(i,j)=D(i,j)-D(401,j); end end for i=1:520 for j=1:520 dfdP(i,j)=dPdX(i)*3*H(i)^2*dHdP(i,j)+H(i)^3*ddPX(i,j)-K*(H(i)-He/Dens(i))*PD1(i,j)-K*Vis(i)*(dHdP(i,j)-He*PD2(i,j)); end end clear i j clear dfdHc dfdHc=3*H.^2.*dPdX-K*Vis.*(1-1./Dens); clear dfdXe dfdXe=K*Vis./Dens*Xe; dgdP=zeros(1,521); if mod(N,2) for i=2:2:N-1 dgdP(i)=dX*4/3; end clear i for i=3:2:N-2 dgdP(i)=dX*2/3; end clear i dgdP(N)=dX/3+(Xe-X(N))/3; else for i=2:2:N-2 dgdP(i)=dX*4/3; end clear i for i=3:2:N-3 dgdP(i)=dX*2/3; end clear i dgdP(N-1)=dX/3+dX/2; dgdP(N)=dX/2+(Xe-X(N))/3; end g=pi/2-dgdP*P'; clear PDF F DF PDF=zeros(521,521); PDF(1:520,1)=dfdHc(1:520)'; PDF(1:520,2:520)=dfdP(:,2:520); PDF(1:520,521)=dfdXe(1:520)'; PDF(521,1)=0; PDF(521,521)=0; PDF(521,2:520)=dgdP(2:520); F(1:520)=-f(1:520); F(521)=g; F=F'; DF=eliminationsolve(PDF,F); TP=DF+TP; Xe=TP(521); N=floor((Xe+5)/dX)+1; TP(N+1:520)=0; TP(TP<0)=0; plot(X(2:520),TP(2:520)); %%%%%%%%%%%%%%%%%% err=1; iter=0; while err>1e-3 TP0=TP; TP=system80three(TP); iter=iter+1; err0=zeros(521,1); for i=1:521 if TP(i)~=0&&TP0(i)~=0 err0(i)=abs(TP(i)-TP0(i))./TP(i); end end err=max(err0); end Ps=0.648*W^0.185*U0^0.275*G^0.391*E/ph; Xs=1.111*W^0.606*U0^-0.021*G^0.077*R/b; P=[0;TP(2:520);0]; p=ph*P; Hc=TP(1); hc=Hc*b^2/R; Defl=D*P-0.25*log(R^2*8*W/pi); for i=1:521 H(i)=Hc+X(i)^2/2+Defl(i)-Defl(401)+S2(i); end !del C:\Users\Administrator\Desktop\aokengdizai\P54.txt fid=fopen('C:\Users\Administrator\Desktop\aokengdizai\P54.txt','wt');%写入文件路径 [m1,n1]=size(P); %获取矩阵的大小，p为要输出的矩阵 format long; for i=1:m1 for j=1:n1 if j==n1 %如果一行的个数达到n1个则换行，否则空格 fprintf(fid,'%f\n',P(i,j)); else fprintf(fid,'%f\t',P(i,j)); end end end fclose(fid); %关闭文件 !del C:\Users\Administrator\Desktop\aokengdizai\H54.txt fid=fopen('C:\Users\Administrator\Desktop\aokengdizai\H54.txt','wt');%写入文件路径 [m1,n1]=size(H); %获取矩阵的大小，p为要输出的矩阵 format long; for i=1:m1 for j=1:n1 if j==n1 %如果一行的个数达到n1个则换行，否则空格 fprintf(fid,'%f\n',H(i,j)); else fprintf(fid,'%f\t',H(i,j)); end end end fclose(fid); %关闭文件