表面纹理方向性研究 - 图文

for n=1:50 Ry=0; for i=1:100 for j=1:100

Ry=Ry+I1(i,j)*I1(i,j+n)/(100*(100-n)); end end

y=abs(Ry-Rz/2.0) if y

r=a1/b1,

对于给定的自相关函数，用线性转化的方法产生随机表面：首先生成一个均值为0，方差为1的随机矩阵Z；表面粗糙度高度用下面的线性转化方法：

B(i,j)???Z(i?k,j?l)k?1l?1nm

这种线性转化会导致x和y方向的自相关函数为：

Rx(?)?(n??*)m,0??*?m

0

*， n??

(m??*)*n,0??*?m

Ry(?)? 0，m??

*n?l?n??0.5x?2??0.5ym?lm2

产生的表面B将是一个高斯分布，均值为0，方差为mn，方向比为γ。 对应的MATLAB程序为：

%产生随机表面B，表面的均值为0，方差为m*n，方向比n/m。 n=1,m=9;

Z=normrnd(0,1,100,100); for i=1:50 for j=1:50 sum=0; for k=1:n for l=1:m

sum=sum+Z(i+k,j+l);

- 27 -

end end B(i,j)=sum; end end mesh(B)；

2、压力的数字化计算方法及其程序

根据上面的雷若方程解出各点的压力。采用有限差分法将上面的微分方程离散化；

有限差分法

Фi-1，j ?y

Фi，j-1 ?x Фi，j ?x Фi，j+1

?y Фi+1，j

非边界处：

?i?1,j??i?1,j???i,j?1??i,j?1??()?()i,j?i,j?y2?y?x2?x ，

?i,j?1??i,j?1?2?i,j?i?1,j??i?1,j?2?i,j?2??2?(2)i,j?(2)i,j?2?x(?x)2?y(?y) ，

边界处：

?i,j?1??i,j?i,j??i,j?1????()i,j?()i,j??x?x?x ， ?x

?i?1,j??i,j?i,j??i?1,j????()i,j?()i,j??y?y?y ， ?y

根据以上的有限差分法离散下面的方程：

?H3?p?H3?pH3?2p3H2?H?pH3?2p3H2?H?p()?()?????0?x12??x?y12??y12??x212??x?x12??y212??y?y ?2p?H?p?2p?H?pH2?3?H2?3?0?x?x?y?y?x?y化简得：

将上面的离散公式代入得：

- 28 -

p(i,j?1)?p(i,j?1)?2p(i,j)H(i,j?1)?H(i,j?1)p(i,j?1)?p(i,j?1)?3?22?x2?x(?x)p(i?1,j)?p(i?1,j)?2p(i,j)H(i?1,j)?H(i?1,j)p(i?1,j)?p(i?1,j)H(i,j)?3?022?y2?y(?y)H(i,j)另?x=?y，化简得：

3?4H(i,j)?p(i,j)?[H(i,j)?(H(i?1,j)?H(i?1,j))]?p(i?1,j)?433[H(i,j)?(H(i,j?1)?H(i,j?1))]?p(i,j?1)?[H(i,j)?(H(i,j?1)?H(i,j?1))]?p(i,j?1)443?[H(i,j)?(H(i?1.j)?H(i?1,j))]?p(i?1,j)?04 根据此公式用MATLAB编程计算每一点的压力，假设共有50*50个点，每点的压力为一个未知数，共2500个未知数，也需要2500个线性方程，非边界点处按上面的公式确定对应的系数， 边界点处按边界条件确定系数。这样产生一个大小为2500维的系数矩阵A，具体的程序如下：

%产生系数矩阵A，a,b分别为线性方程组的个数和每个方程对应的系数；假设有50*50个点组成数组，则共有2500个数据，i,j为对应点所在的行和列。 A=zeros(2500,2500); for i=1:50

for j=1:50 a=(i-1)*50+j;

if i==1 %第一行点的压力系数，边界条件导数为0； for b=1:2500 if b==a A(a,b)=1; end if b==a+50 A(a,b)=-1; end end end

if i==50 %最后一行点的压力系数，边界条件导数为0 for b=1:2500 if b==a A(a,b)=1; end if b==a-50 A(a,b)=-1; end end end

if i>=2&i<=49

- 29 -

if j==1|j==50 %第一列和最后一列点系数，边界条件为常数； for b=1:2500 if b==a A(a,b)=1; end end end

if j>=2&j<=49 %中间点所按照有限差粉法确定系数 for b=1:2500 if b==a-50

A(a,b)=H(i,j)-3/4*(H(i+1,j)-H(i-1,j)); end if b==a-1

A(a,b)=H(i,j)-3/4*(H(i,j+1)-H(i,j-1)); end if b==a

A(a,b)=-4*H(i,j); end if b==a+1

A(a,b)=H(i,j)-3/4*(H(i,j+1)-H(i,j-1)); end if b==a+50

A(a,b)=H(i,j)-3/4*(H(i+1,j)-H(i-1,j)); end end end end end end

方程组的常数矩阵F按照边界点对应的方程为常数，其它都为0；这样得到的常数项F的程序为：

%方程组常数项F的确定，常数项为一个2500行，一列的数组；压力流量方程组中除了第一列和最后一列的数据分别为PA，PB，其他都为零 PA=10;PB=2; for i=1:2500 for m=1:50

if i==(m-1)*50+1 F(i,1)=PA; end if i==m*50

F(i,1)=PB; end end end

- 30 -

这样根据MATLAB中线性方程组的解法，得到2500个点对应的压力X=A\\F；并将X转化为对应的两维矩阵，转化过程为： %方程组的解X，将其转化为二维矩阵P； X=A\\F; for i=1:50

for j=1:50

k=(i-1)*50+j; P(i,j)=X(k,1); end end

- 31 -