分形表面

分形表面（精选4篇）

分形表面 篇1

0 引言

粗糙表面的接触行为对摩擦学研究具有重要意义。利用分形理论建立的粗糙表面接触模型可以提供多尺度的接触行为预测分析。Majumdar等[1]利用分形几何可以用来描述多尺度粗糙面的优点,将传统的统计参数以分形参数取代,建立了二维粗糙表面的微观接触模型(M-B模型)。Yan等[2]将满足分形性质的Weierstrass-Mandelbrot (W-M) 函数推广至三维的粗糙表面,可以描述更精确的粗糙表面形貌,但是其接触模型(Y-K模型)只探讨了具有弹性-完全塑性(elastic-perfectly plastic)性质的材料,忽略了弹塑性变形与真实的接触情况。盛选禹等[3]利用M-B模型,考虑了微凸体弹性和塑性两状态,对分形表面的静摩擦因数进行了预测。在此基础上,You等[4]根据M-B模型并考虑了微凸体弹性、弹塑性和塑性状态,分析了不同分形参数下的粗糙表面静摩擦特性。上述研究并未针对三维粗糙表面并考虑弹性、弹塑性和塑性的完整状态进行接触分析。本文针对三维分形粗糙表面的接触问题,进一步修正了微凸体的接触模型,模型包含了微凸体接触的弹性-弹塑性-塑性的所有状态,并对接触力进行了理论预估和实验对比分析。

1 三维分形表面的弹塑性接触分析

实际工程领域中,假设工程表面微观轮廓具有连续性时,根据观察可看到当表面微观轮廓被重复放大时,更加精细的结构会不断出现,而且轮廓永远是不光滑的,在任何点均不存在切线,所以轮廓是处处不可微的。另外,当轮廓被放大时,放大后的表面和原始表面的概率分布非常相似,因此大多数轮廓具有自仿射性。而分形几何学中的W-M函数,具有连续性、处处不可微、自仿射性的数学特征,十分适合描述多尺度的粗糙表面。W-M函数表示为[2]

$w(x)={\displaystyle \sum _{n=-\infty }^{\infty }[}\gamma {}^{(D-2)n}(1-\text{e}{}^{\text{i}\gamma {}^{n}x})\text{e}{}^{\text{i}\text{ϕ}{}_n}]$ (1)

式中,x为粗糙表面测量坐标;w为x的复合函数;D为二维粗糙表面分形维数(1<D<2);γ为大于1的常数,γn表示随机轮廓的空间频率,对于服从正态分布的随机表面,γ =1.5;ϕn为一随机项。

取w(x)的实部可得到二维分形表面轮廓的函数:

$\begin{array}{l}z(x)=\mathrm{Re}[w(x)]=\\ {\displaystyle \sum _{n=-\infty }^{\infty }[}\gamma {}^{(D-2)n}(\cos \text{ϕ}{}_n-\cos (\gamma {}^{n}x+\text{ϕ}{}_n))](2)\end{array}$

式中,z(x)为粗糙表面轮廓高度;G为粗糙表面特征长度尺度系数。

Ausloos等[5]基于W-M函数,并引入了几个变量来描述分形表面的多尺度波的叠加特征,从而构造出了各向同性的三维分形表面在极坐标下的表达式:

$\begin{array}{l}z(\rho ,\theta )=(\frac{\ln \gamma }{{\rm M}}){}^{1/2}{\displaystyle \sum _{m=1}^{{\rm M}}A}{}_m{\displaystyle \sum _{n=-\infty }^{\infty }\{}(k\gamma {}^n){}^{D-3}\cdot \\ [\cos \text{ϕ}{}_{m,n}-\cos (k\gamma {}^n\rho \cos (\theta -\alpha {}_m)+\text{ϕ}{}_{m,n})]\}(3)\end{array}$

式中,ρ、θ为表面点的极坐标,与笛卡儿坐标转换关系为$\rho =\sqrt{x{}^2+y{}^2}‚\theta =\arctan (y/x)$;M为分形表面参与叠加的多尺度波的层次;k为同一尺度的波的数量,且k=2π/L,L为采用长度;D为三维表面的分形维数,且2<D<3。

通过引入一个长度参数G,并令Am=2π(2π/G)(2-D),αm=πm/M,然后将Am、αm、ρ、k根据极坐标与笛卡儿坐标转换关系,代入式(3)得到笛卡儿坐标下的三维分形表面表达式[2]:

$\begin{array}{l}z(x,y)=L(\frac{G}{L}){}^{(D-1)}(\frac{\ln \gamma }{{\rm M}}){}^{0.5}{\displaystyle \sum _{m=1}^{{\rm M}}{\displaystyle \sum _{n=0}^{n{}_{\max }}\{}}\gamma {}^{(D-2)n}[\cos \text{ϕ}{}_{m,n}-\\ \cos (\frac{2\text{π}\gamma {}^n\sqrt{x{}^2+y{}^2}}{L}\cos (\arctan (\frac{y}{x})-\frac{\text{π}m}{{\rm M}})+\text{ϕ}{}_{m,n})]\}(4)\end{array}$

$n{}_{\max }=\text{Ι}\text{n}\text{t}(\frac{\ln (L/L{}_{\text{s}})}{\ln \gamma })$

式中,Ls为最小波长。

图1所示即为根据式(4)并通过MATLAB 7.0生成的三维分形表面,具体参数为D=2.4,G=1.36×10-2nm,L=1mm,Ls=5nm,M=10,γ=1.5。

2 三维分形表面的弹塑性接触分析

通常情况下两个粗糙表面接触(图2)可简化为一个刚性平面和一个等效粗糙平面的接触,该等效粗糙平面的弹性模量E*和硬度H*可表示为[6]

式中,下标A,B分别表示两粗糙平面;E、ν、H分别表示弹性模量、泊松比和硬度。

因此,单个接触微凸体的变形量δ和等效半径R可表示为

式中,r′、A′分别为微凸体截断圆半径和面积,A′=π(r′)2。

两粗糙表面的接触磨损可简化为一刚性平面和一粗糙平面的接触磨损。进一步假设粗糙表面是统计各向同性的,那么分析上述刚性平面-粗糙面的接触,可通过分析每个微凸体的接触情况,采用由局部到整体的方法来实现。在Y-K模型中,忽略了微凸体的弹塑性变形,简单将其变形分为两种情况:纯弹性和纯塑性变形,并根据Hertz理论得到开始塑性变形的临界变形量δc和临界接触面积A′c,分别表示为[2]

$\delta {}_{\text{c}}=(\frac{\text{π}{\rm K}{\rm H}}{2E{}^{*}})R$(7)

K=0.454+0.41ν

$A^{\prime }{}_{\text{c}}=[2{}^{11-2D}\text{π}{}^{D-4}G{}^{2D-4}\ln \gamma (\frac{E{}^{*}}{{\rm K}{\rm H}}){}^2]{}^{1/(D-2)}$(8)

多数文献将微凸体变形分为两种情况:纯弹性和纯塑性变形,忽略了弹塑性变形情况。Kogut等[7]应用有限元对球—面接触进行分析,得到包含弹塑性变形的更为完整的接触模型,即当δ≤δc时球形微凸体处于弹性状态;当δc<δ≤6δc时球形微凸体处于第一弹塑性状态;当6δc<δ≤110δc时球形微凸体处于第二弹塑性状态;当110δc<δ时球形微凸体处于完全塑性状态,且在各个状态下接触力和接触面积与微凸体变形量有如下关系:

(1)当δ≤δc时,微凸体处于弹性状态,且接触面积Ae和载荷Fe分别可表示为

(2)当1≤δ≤6δc时,微凸体处于第一类弹塑性状态,且接触面积Aep1和载荷Fep1分别可表示为

(3)当6δc<δ≤110δc时,微凸体处于第二类弹塑性状态,且接触面积Aep2和载荷Fep2分别可表示为

(4)当δ>110δc时,微凸体处于完全塑性状态,且接触面积Ap和载荷Fp分别可表示为

将式(6)代入式(9)并根据A′=π(r′)2得到由截断接触面积表示的微凸体接触载荷,则Fe可重新表示为

$F{}_{\text{e}}=\frac{1}{3}E2{}^{6-D}\text{π}{}^{(D-4)/2}(\ln \gamma ){}^{1/2}G{}^{D-2}(A^{\prime }){}^{(4-D)/2}$(13)

同理,Fep1、Fep2和Fp均可重新表达为

$F{}_{\text{e}\text{p}1}=1.2377E(\frac{{\rm K}{\rm H}}{E}){}^{0.144}2{}^{-0.85D+3.4}\text{π}{}^{0.425D-1.7}\cdot$

(lnγ)0.425G0.85D-1.7(A′)(-0.425D+1.85) (14)

$F{}_{\text{e}\text{p}2}=1.3439E(\frac{{\rm K}{\rm H}}{E}){}^{0.474}2{}^{-0.562D+2.104}\text{π}{}^{0.263D-1.052}\cdot$

(lnγ)0.263G0.526D-1.052(A′)(-0.263D+1.526) (15)

Fp=HA′ (16)

另外,综合式(6)和式(7),并利用关系δ=6δc和δ=110δc得到由完全弹性区向第一弹塑性接触区转变的临界接触面积A′c2,即

A′c2=61/(2-D)A′c (17)

以及由第一弹塑性接触区向第二弹塑性接触区转变的临界接触面积A′c1,即

A′c1=1101/(2-D)A′c (18)

那么,总的真实接触面积Ar包括四个区域的微凸体接触面积之和,即处于完全弹性状态的所有微凸体接触面积Ae,处于第一类弹塑性状态的所有微凸体接触面积Aep1,处于第二类弹塑性状态的所有微凸体接触面积Aep2,处于完全弹性状态的所有微凸体接触面积Ap。具体表达式如下:

Ar=Ae+Aep1+Aep2+Ap=

∫${}_{A^{\prime }{}_{\text{c}}}^{A^{\prime }{}_{\text{L}}}$n(A′)A′dA′+∫${}_{A^{\prime }{}_{\text{c}2}}^{A^{\prime }{}_{\text{c}}}$n(A′)A′dA′+

∫${}_{A^{\prime }{}_{\text{c}1}}^{A^{\prime }{}_{\text{c}2}}$n(A′)A′dA′+∫${}_{A^{\prime }{}_{\text{S}}}^{A^{\prime }{}_{\text{c}1}}$n(A′)A′dA′ (19)

那么总的接触载荷FN包括处于四个区域的所有微凸体接触载荷之和,即

FN=∫${}_{A^{\prime }{}_{\text{c}}}^{A^{\prime }{}_{\text{L}}}$n(A′)FedA′+∫${}_{A^{\prime }{}_{\text{c}2}}^{A^{\prime }{}_{\text{c}}}$n(A′)Fep1dA′+

∫${}_{A^{\prime }{}_{\text{c}1}}^{A^{\prime }{}_{\text{c}2}}$n(A′)Fep2dA′+∫${}_{A^{\prime }{}_{\text{S}}}^{A^{\prime }{}_{\text{c}1}}$n(A′)FpdA′ (20)

式中,A′L、A′S分别为最大和最小接触微凸体截断面积。

3 结果分析

下面将以Y-K模型和修正的Y-K模型进行对比分析。真实接触面积Ar,法向接触载荷FN和位移δ以如下方式进行量纲一化:

式中,Aa为名义接触面积。

为了验证模型的正确性,下面对照前人的实验结果进行分析。Kucharski通过实验研究了粗糙表面接触时载荷和实际接触面积的关系[8]。文献[9]对Kucharski所做实验的材料物理特性进行了总结,即实验采用材料为钢,弹性模量为200GPa,泊松比为0.3,屈服强度为400MPa,硬度(HV)为1.12GPa,粗糙表面分形特征参数为D=2.3,G=1.36×10-11m。本文采用上述参数,基于修正后的Y-K模型,与Y-K模型、有限元分析及实验的结果进行比较,如图3所示。显然,有限元分析结果与实验测试值比较接近,而Y-K模型的结果则与实验值相差甚远,尤其是在载荷较大时。修正Y-K模型的结果,虽然与实验值有一定差距,但与Y-K模型的结果相比已有了相当的改善。

1.Y-K模型 2.修正后模型 3.实验值 4.有限元分析值

另外,Buzzio[10]也针对分形表面的接触作了实验分析。其实验所用材料的物理参数及表面分形参数也参见文献[9],即所用材料的弹性模量为0.8GPa,泊松比为0.3,硬度(HV)为45MPa。表面特征参数为D=2.3,G=1.55×10-10m,名义接触面积为A0=4.8×10-12m2。图4所示为采用上述参数进行分析,并与Y-K模型及实验值作对比的结果。此种情况下,修正Y-K模型与Y-K模型的计算结果均与实验值表现出了一致性。其原因主要是:在该种材料及载荷情况下,临界微凸体接触面积A′c大于最大微凸体接触面积A′L,意味着所有微凸体都处于完全塑性状态,所以此种情况下修正后Y-K模型和未修正时相同,因为修正后Y-K模型只影响弹塑性状态的微凸体。综上,修正后Y-K模型的计算结果较修正前得到了相当的改善。

1.修正Y-K模型 2.实验值 3.Y-K模型

4 结论

(1)采用分形理论研究材料属性、表面形貌及接触载荷的关系,通过考虑微凸体接触的弹性-弹塑性-塑性的所有状态来修正Y-K模型,并将所得结果与实验结果进行对比,发现修正后Y-K模型较Y-K模型的计算精度有了相当的改善,但与实验值相比仍存在一定差距。

(2)针对分形粗糙表面接触行为的预测,可构建更精确的摩擦磨损的理论模型,可在一定程度上减少和替代试验测试,从而为摩擦学设计提供更为精细的计算方法。

参考文献

[1]Majumdar A,Bhushan B.Fractal Model of Elastic-Plastic Contact between Rough Surfaces[J].ASME Journal of Tribology,1991,113:1-11.

[2]Yan W,Komvopoulos K.Contact Analysis of Elas-tic-Plastic Fractal Surfaces[J].Journal of AppliedMechanics,1998,84(7):3617-3624.

[3]盛选禹,雒建斌,温诗铸.基于分形接触的静摩擦系数预测[J].中国机械工程,1998,9(7):15-18.

[4]You J M,Chen T N.A Static Friction Model forthe Contact of Fractal Surfaces[J].Proceedings ofthe Institution of Mechanical Engineers,Part J:Journal of Engineering Tribology,2010,224:513-517.

[5]Ausloos M,Berman D H.A Multivariate Weier-strass-Mandelbrot Function[J].Proceedings of theRoyal Society of London.Series A,Mathematicaland Physical Sciences,1985,400(1819):331-350.

[6]Greenwood J A,Williamson J B P.Contact ofNominally Flat Surfaces[J].Proceedings of theRoyal Society of London,Series A,Mathematicaland Physical Sciences,1966,295(1142):300-319.

[7]Kogut L,Etsion I.Elastic-Plastic Contact Analy-sis of a Sphere and a Rigid Flat[J].ASME Journalof Applied Mechanics,2002,69:657-662.

[8]Kucharski K,Klimczak T,Polijaniuk A,et al.Fi-nite-element Model for the Contact of Rough Sur-faces[J].Wear,1994,177:1-13.

[9]Sahoo P,Ghosh N.Finite Element Contact Analy-sis of Fractal Surfaces[J].J.Phys.D:Appl.Phys.,2007,40(14):4245-4252.

[10]Buzzio R,Boragno C,Valbusa U.Contact Me-chanics and Friction of Fractal Surfaces Probed byAtomic Force Microscopy[J].Wear,2003,254:917-923.

各向异性分形表面建模研究 篇2

随着机械系统进一步向微型化发展，针对机械零件、刀具等表面几何形貌的研究日益受到关注[1,2]。分形几何是粗糙表面表征和建模的重要几何模型，目前已深入开展了基于分形表面的接触、摩擦、磨损等性能的研究。杨红平等[3]推导了基于分形几何的结合面法向接触刚度计算模型；Ji等[4]基于分形几何，建立了粗糙表面的接触热传导预测模型；黄健萌等[5]基于有限元法，建立了在充分考虑粗糙表面接触微凸体间的相互作用及接触界面摩擦热流耦合等影响因素情况下的分形表面与理想刚性平面接触时的热力耦合模型。目前的研究工作中，一般将各向同性的分形表面作为几何模型，然而，实际零件表面存在大量的各向异性特性，研究、完善各向异性分形表面建模算法具有重要的理论和实际意义。文献[6,7,8]进行了各向异性分形表面表征和建模的研究。其中，Wu[8]基于傅里叶变换分析了多种类型的各向异性分形表面，给出了相应的频域表达式以及建模算法。与上述研究工作相比，本文基于二维离散傅里叶变换（discrete Fourier transform,DFT）的投影性质和旋转性质，通过直接修改功率谱图进行各向异性分形表面建模。该表面建模算法在避免了数学复杂性的同时，可精细控制任意期望方向的分形参数。

1 基本理论

1.1 二维离散傅里叶变换

二维DFT建立了函数空域（时域）和频域的一一对应关系，对于采样点数为N2的二维空域信号f(i,k），若令其傅里叶系数为F(u,v），则其离散傅里叶变换可表示为[9]

式中，N为沿一个方向的采样点数，为方便计算，一般取2的整数次幂；j为虚数单位，j 2=-1;i、k为空域采样序号；u、v为频域采样序号。

相应的离散傅里叶逆变换（inverse discrete Fourier transform,IDFT）为[9]

对于实二维信号f(i,k),F(u,v）一般为复数，因此可表示为

F(u,v）共轭对称，即幅度谱

相位谱

并有以下投影性质和旋转性质[9?10]:

(1）投影性质。若函数g(x,y）的傅里叶变换记为G(u,v），且

(2）旋转性质。若函数g(x,y）的傅里叶变换记为G(u,v），在空域中将g(x,y）旋转角度θ，则G(u,v）在频域中也旋转了θ。

1.2 各向同性分形表面建模

基于分数布朗运动分形表面的功率谱，采用二维IDFT的分形表面建模算法具有精度高、速度快及建立的分形表面的高度均方根偏差Sq可控的优点[11]。

分形表面的功率谱满足[12]:

式中，E（*）表示数学期望；D为分形维数；C为尺度系数。

u=v=0时，令F(u,v)=0，此时合成表面各采样点高度的均值为0[13]。若假设分形表面具有平稳性，即合成的各表面采样点高度的均值和均方根偏差不变，则式（4）可写为[11,13]

令表面傅里叶系数的相位谱φ（u,v）为满足均匀分布的随机值，利用式（2）即可完成各向同性分形表面的建模。图1a所示为基于式（5）建立的分形表面，图1b所示为其功率谱。建模时，N=64,D=2.6,C2=4.1457×107。由于D和C是与尺度无关的分形参数，因此本文所建立的分形表面的高度单位为任意长度单位。同时，由文献[14]可知，合成粗糙表面各采样点高度的均方根偏差仅与采样点数N和|F(u,v）|2有关，因此即使在进行各向异性分形表面建模时，也可保证其平稳性。

2 各向异性分形表面建模

2.1 沿空域i轴异性的分形表面

由二维DFT的投影性质可知，功率谱图中，沿u轴方向的功率谱即为表面所有采样点向空域i轴投影后，再进行一维DFT的结果。图2a中，N=64，沿u轴方向的分形参数设置为Da=2.2,Ca2=1×109，其余方向分形参数与图1表面相同，即D=2.6,C2=4.1457×107；图2b所示为经二维IDFT获得的分形表面空域几何形貌，该表面沿i轴方向表现出明显的异性特性。

2.2 沿其他方向异性的分形表面

若要使得合成后分形表面空域的异性方向与i轴夹角为θ，那么只需在功率谱图上沿与u轴夹角为θ的方向，设置相应的各向异性参数即可。该步骤等价于：首先根据二维DFT的投影性质，在表面功率谱图的u轴上设置各向异性参数；然后，根据二维DFT的旋转性质，将具有各向异性特征的谱线旋转θ角。图3a所示的分形表面的异性方向与i轴夹角为45°，图3b为其功率谱；图3c所示的分形表面的异性方向与i轴夹角为135°，图3d为其功率谱。

上面讨论中，分形表面的各向异性性质仅出现在某单一方向，实际上，基于二维DFT的投影性质和旋转性质，通过修改给定方向上的功率谱图，可以方便地合成在多个方向上具有各向异性特性的分形表面。例如，图4a、图4b所示分形表面的异性方向同时平行于i轴和k轴；图4c、图4d所示表面的异性方向轴与i轴的夹角分别为45°和135°，虽然在这两个方向上的功率谱相同，但受DFT随机相位的影响，其空域几何形貌并不相同。

3 影响表面形貌的因素

上节为了清晰地体现表面的各向异性特征，在各向异性方向设定的分形维数较低（沿各向异性方向的分形维数Da=2.2，其余方向的分形维数D=2.6），尺度系数较大（沿各向异性方向的尺度系数Ca2=1×109，其余方向的尺度系数C2=4.1457×107）。这导致了沿异性方向能量占明显优势，合成的表面几何形貌较为平坦[13]。在图5中，令沿u轴的分形维数Da=2.4，尺度系数Ca2=1.658 28×108，获得了较为崎岖不平的各向异性分形表面。

另外，以上表面除了在频点（u=0,v=0）之外，其余各处功率谱都不为0，因此在外观上出现了一个长波长的几何特征。根据二维DFT的定义和有关性质，可将沿各向异性方向上更多的低频成分置为0，从而获得波动更为剧烈的各向异性分形表面。

图6a所示的功率谱中，3个低频成分即（u=0,v=0),(u=1,v=0）和（u=-1,v=0）处的|F(u,v）|2=0，图6b为其空域几何形貌；图6c中，5个低频成分，即（u=0,v=0),(u=1,v=0),(u=-1,v=0),(u=2,v=0）和（u=-2,v=0）处的|F(u,v）|2=0，图6d为其空域几何形貌。

图7a所示的功率谱中，3个低频成分即（u=0,v=0),(u=1,v=1）和（u=-1,v=-1）处的|F(u,v）|2=0，图7b为其空域几何形貌；图7c中，5个低频成分即（u=0,v=0),(u=1,v=1),(u=2,v=2),(u=-1,v=-1）和（u=-2,v=-2）处的|F(u,v）|2=0，图7d为其空域几何形貌。

4 结论

(1）基于二维离散傅里叶变换的投影性质和旋转性质，可灵活、方便地合成在指定的一个或多个方向具有给定分形参数的各向异性分形表面。

(2）功率谱图中，修改沿给定各向异性轴方向的分形参数，可增加或减少沿该方向表面几何形貌与其他方向的差别。

表面粗糙度的若干分形模拟方法 篇3

在各种生产过程中, 任何物体都不可避免的会产生表面, 所谓表面是指该物体材料与周围介质 (通常指空气) 之间的边界。经过加工的表面一般认为有粗糙的表面与光滑的表面之差别, 这可以用视觉或触觉的方法直观地加以区分:光滑表面大多能反射光线, 在表面上看不出有花纹, 手感平滑, 阻力小;粗糙表面不反光, 能看见明显的纹理和沟痕、凹坑, 用手指触摸时感到凹凸不平, 摩擦力大。但是这种判别带有主观因素, 而且只是定性的相对说明比较光滑或粗糙一些, 没有确切的表达表面特征之细节和全貌。

制件在制造过程中产生的表面几何形状总是不够完善的, 如加工后的实际表面形状相对于理想表面形状存在一定的偏差。通过测量会发现, 在表面上有一系列不同间距和高度的峰谷所组成的不规则几何形状叠加在一起的复杂表面结构。

对于这种复杂结构, 现今一般采用三类结构型式——形状误差、表面波纹度和表面粗糙度来描述。这三类表面几何形状偏差在一个表面上并非孤立存在, 大多数加工表面常受其综合影响。但是各类偏差形成的原因和特性以及他们与各种使用功能的因果关系均不相同, 而且受测量手段的约束, 因此目前仍采取分别的评定。

其中关于表面粗糙度完全是由加工方法固有的内在作用所产生, 是制件在加工过程中由实际加工过程中由实际加工介质--切削刀具、磨粒、喷丸等在完工表面上留下的微观不平度。例如切削过程中的残留面积, 刀具对制件表面的摩擦、切屑分裂时材料的塑性变形以及刀瘤和灼伤等因素, 形成各种形式凹凸不平的微细加工痕迹。采用不同的工艺方式便构成特定的表面微观几何结构, 这种微观结构与分形非常类似, 并且有着密切的联系。

随着计算机技术的迅速发展, 数值模拟越来越成为摩擦学研究中的重要手段。由于表面形貌随机性强、形成过程复杂、影响因素多, 使得对表面形貌的定量研究主要依靠测量之后的数据分析, 属于“事后”分析。事实上, 表面形貌在摩擦过程中始终是变化的, 表面形态的演化在时间上具有瞬变性, 在空间上具有离散性, 我们很难记录某一表面形貌在摩擦过程中的演化全过程。随着摩擦学研究的发展, 摩擦磨损行为描述由定性向定量化发展, 表面形貌的研究也大量借助计算机进行仿真模拟和分析, 这样就可以在摩擦过程之前分析表面形貌的演变, 实现“事前”分析。另外, 表面形貌模拟也是虚拟摩擦学的重要技术基础。

大多数固体粗糙表面结构在数埃至数百埃的范围内为分形, 表面的许多功能与其分形结构有关, 其分形表面具有测量的不确定性。根据实际测量发现, 不同的尺度得到的表面粗糙度值是变化的, 这说明传统的表面粗糙度表征参数是基于统计学的, 它们将随着测量尺度的变化而表现出不稳定性, 即具有尺度向关性的特征。正是由于表征参数有尺度相关性, 要求我们从不同的测量分辨率下得到粗糙度的表述, 从而有了下述各种粗糙表面模拟方法。

2 粗糙表面模拟方法

2.1 分形布朗函数模拟法[40]

分形布朗运动 (Fractal Brownian Motion, fBm) 是B.Mandelbro和VanNess于1968年首先提出的, 用于模拟各种具有分形特征的噪声等, 现已成为一个能反映广泛的自然物体性质的数学模型与粗糙表面的许多数学特征相近 (如分形维数、统计自相似、增量平稳等) 。而fBm数学模型的特征可由赫斯特指数H (或分形维数D) , 标准差δ表示。对模拟表面而言, H是粗糙表面复杂情况的一种抽象和概括, 当H越大时, 分形布朗运动变化趋于平缓, 而当H变小时, 变化趋于精细。而δ反应着局部表面总的形状起伏特征。我们只要抽取某个区域的特征参数H和δ, 即可用fBm模拟生成该区域的多分辨率多细节层次的表面。

分形表面的模拟技术实质上就是在数学上合成一个满足fBm特性的空间曲面。分形布朗运动是现代非线性时序分析中的重要随机过程, 它能有效地表达自然界中许多非线性现象, 目前已发展了多种fBm合成技术, 主要有泊松阶跃法、傅里叶滤波法、随机中点位移法、逐次随机增加法等。

2.2 逆傅里叶变化模拟法

逆傅里叶变化模拟法 (Inverse Fourier Transformation, 简称IFT法) 可以模拟表面轮廓和表面形貌。它模拟表面轮廓的步骤如下:生成一组{0, 1}分布的高斯白噪声hk, k=0, 1, …, N-1然后对其做Fourier变换, 得到Fourier系数为

对Fourier系数的实部和虚部进行修正, 使其满足下面的幂律关系, 即

由此得修正后的Fourier系数

最后, 对修正的傅里叶系数做逆Fourier变换, 就得到一个模拟轮廓。

模拟粗糙表面形貌的步骤如下:对二维高斯白噪声作Fourier变换, 使其功率普满足式所示的幂率形式, 对满足条件的S (ω) 做逆Fourier变换, 即可得到fBm序列X (t) 。二维的逆Fourier变换为

其中的变换满足幂型条件

akl是通过对二维白噪声做Fourier变换, 然后除以得到的。IFT法可以通过高频成份的调节来任意接近fBm, 生成的表面具有很好的fBm特性。

2.3 函数模拟法

逆傅里叶变换法实质上是一个频率域中的线性累加, 而随机化的函数则是一个频率域中的几何累加, 它是布朗函数的另一种构造, 其表达式为

上式中An是一个{0, 1}的独立正态分布随机变量;是[0, 2π]区间均匀分布的相位。在一定的尺度变化Nmin~Nmax范围内, W-M法可通过预先设定的系数An和生成fBm曲面。W-M法生成的图像是由预定系数计算的, 因此比较容易控制产生图像的边界, 表面分辨率可以通过累加分量数调节。

2.4 分形插值模拟法

如前文所述, 粗糙表面的轮廓高度变化是一非稳定随机过程, 具有多重尺度特征, 因而用一定分辨率的仪器所测得的表面形貌只反映了该分辨率水平的粗糙度, 不能完全反映表面粗糙度的真实信息, 这就是粗糙度测量的尺度相关性表现。显然, 基于这种情况采集的数据而表征的形貌参数存在一定偏差, 这会影响粗糙表面之间接触和摩擦的认识结果和计算精度。如果能用有限的实验数据模拟出更加精细的粗糙结构, 则可解决仪器分辨率不足所带来的粗糙度测量尺度向关性的问题。另外, 表面形貌模拟也会遇到一种情况, 那就是已知轮廓基线函数时, 如何模拟粗糙轮廓的问题, 分形插值模拟是解决这两个问题的方法之一。

2.4.1 经典插值方法

插值方法是数值模拟与数值计算中常用方法。传统插值方法要求插值函数与被插值函数在插值基点处具有相同的函数值, 甚至到某阶导数。例如, 将原始数据点在坐标中标出, 然后从几何上分析数据点之间的关系, 构造出一个次数尽可能低的多项式或者是一个基本初等函数的复合函数P (x) , 使多项式f (x) 或复合函数P (x) 的曲线的若干坐标点处于原始数据点具有相同值。这里构造出的多项式或复合函数就称为插值函数。

插值函数的选取方法很多, 可以选取代数多项式, 选取三角多项式或有理函数, 也可以选取某区间上的任意光滑函数或者分段光滑函数。最常用、最基本的插值函数是次代数多项式, 即

其中, a0, a1, a2, …an是实数, 把多项式函数的插值问题称为多项式插值。

当数据点的区间范围较大时, 单独采用一个高次的多项式插值, 拟合误差可能较大。为此, 常采用分段插值方法, 即对于不同的数据点范围, 利用不同的插值多项式。分段插值相当于用若干条多项式曲线连接而成, 所以拟合曲线一般是不光滑的, 为此采用样条插值来克服分段插值函数的拟合不光滑问题。

2.4.2 自仿射分形插值函数模拟法

给定闭区间I=[0, 1], 定义映射

讨论当Li (x) , Fi (x, y) 都取线性函数时的情形, 即ωi有如下形式:

根据自仿射分形插值函数的迭代函数系系数的计算公式, 其中参数di为垂直压缩因子,

由上述公式确定的FIF:记作y=f (x) ,

连续函数f (x) 的图像是迭代函数系{ωi}的不变集。

给定三个插值结点, 给定,

可以确定如下两个仿射变换,

当压缩因子取不同的值时FIF构造出的曲线有不同的形状, 当时图像较为光滑, 而当时图像非常粗糙。

变换w1将整个图像变换到直线x=1/2的左侧, 而变换w2将整个图像变换到直线x=1/2的右侧。根据这三个插值结点以及两个变换可以唯一的确定一个分形插值函数, 这样第一次迭代中由3个点迭代得出5个点, 其中仅有两个点为新出现的点, 记作, , 第二次迭代中将上一次的5个点继续用w1w2作用, 得到四个新出现的点 (其他的点没有变化) , 记作

易见每次迭代时所产生的新点个数为上一次两倍, 即第n次有总计2n个新的结点。记第次迭代所生成的所有新出现的点为

n (其中) , 容易得到下列结论:在最终的分形插值曲线的图象

在最终的分形插值曲线y=f (x) 的图象G上.换句话说, 图像G上是所有经过无穷次迭代的不动点集合。

从而可以看出此最大值连续依赖于压缩因子, 并且能够按照此结论类似求得此种分形插值曲线的相关粗糙度参数。

3 当前分形的发展

近年来, 分形的研究受到了非常广泛的重视, 其原因在于分形既有深刻的理论意义, 又有巨大的实用价值。分形所涉及的领域极为广泛, 在数学、物理、化学、材料医学、地理学、天文学、计算机科学以及经济领域等, 甚至艺术、哲学中都不断有大量的学术论文产生。目前, 分形理论已得到广泛应用, 推动了许多学科的发展。但从国内外分形理论发展来看, 仍有许多本质性的问题有待进一步解决, 作为分形理论的核心, 分形数学理论有待于进一步完善。例如如何判断一个对象是分形, 分形维数的本质是什么, 分形生成的原因是什么, 这些问题也有待进一步解决。另外, 分形与其他学科结合过程中的一些问题如分形物理理论, 分形高分子理论等重要研究方向和课题都有待进一步深入。

参考文献

[1]葛世荣朱华.摩擦学的分形[M].北京:机械工业出版社, 2005.

[2]袁长量表面粗糙度及其测量[M]机械工业出版社1989.

[3]毛起广.表面粗糙度的评定和测量[M].北京:机械工业出版社, 1991.2.

[4]高永生李柱表面统计粗糙度理论与评定方法计量学报198910 (3) 176-181.

[5]俞汉清表面粗糙度标准及应用[M].北京:中国计量出版社, 1996.12.

[6]孙霞吴自勤黄分形原理及其应用[M]合肥:中国科学技术大学出版社, 2003.10.

分形表面 篇4

表观上看似光滑平整的机械加工表面经放大后呈现出大小各异、复杂排列的凸峰和凹谷,说明表面形貌实际上是粗糙和不规则的。机械工程中的表面几何形貌对构件接触处的摩擦、磨损、润滑、密封以及接触性能有很大的影响。研究表明[1],机械加工表面和摩擦磨损表面等粗糙表面轮廓具有非平稳性、自相似性和多尺度特性,分形理论是描述这些特征的一种有效途径;然而,分形理论用于描述机械加工表面时需要分形维数与特征长度尺度参数两个表征参数。准确地计算出其分形维数尤为关键,因为分形维数反映了粗糙表面轮廓结构的复杂程度,定性地刻画了高频成分所占比重。分形维数是机械结合面接触刚度[2,3,4]、阻尼能耗[5,6]等模型中一个主要的参数,接触刚度与表面分形维数间存在非线性关系(分形维数在1.1~1.4之间)或近似线性关系(分形维数在1.4~1.9之间)[7];在实际应用时,对轮廓维数的计算是必不可少的。

国内外已有许多学者进行了分形维数计算的研究,目前,常用的分形维数计算方法有盒维数法、尺码法、均方根法、协方差法、功率谱密度(PSD)法、结构函数法。李成贵等[8,9]对已知分形维数为1.6的W-M函数所生成的轮廓,分别采用功率谱密度法和结构函数法计算其分形维数,计算结果分别为1.65、1.63。王建军等[10]对尺码法、盒维数法、R/S法、结构函数法以及功率谱密度法5种方法的原理进行了阐述并比较了后3种方法对理论维数分别为1.2、1.5、1.8的W-M函数的计算结果,结果表明,结构函数法的精度最高,功率谱密度法次之,R/S法精度最低且计算量大,其精度随着理论分形维数的减小而大幅降低。葛世荣等[11]比较了尺码法、盒维数法、方差法和结构函数法等方法计算分形维数分别为1.2、1.5、1.8的W-M函数的计算结果,结果表明,尺码法的误差最大,方差法和盒维数法的误差在10%以上,结构函数法的误差较小,计算结果分别为1.164、1.455、1.761。另外,他们还提出了均方根法[12,13],并将均方根法与结构函数法进行了比较,结果表明,两种方法均有较好的效果;对于较小的分形维数,均方根法的计算精度较高,而对于较大的分形维数,结构函数法的计算精度较高。此外,还可以通过表面轮廓的高度标准差σ和斜率标准差σ′组成方程组解出分形参数D和G[3],本文称之为方程组法。

上述各方法中,由Majumdar等[4]提出的结构函数法对分形维数的计算精度最高,已得到了较多的应用,然而,寻求更高精度的计算方法来提高结合面建模精度是非常必要的。小波具有多尺度分析的能力,可用于处理多尺度自相似问题。王安良等[14,15]提出了用小波变换计算粗糙表面分形维数的方法,根据各分解层小波变换系数的一范数与分解层数k间的关系通过最小二乘拟合求出分形维数,对于由W-M函数及M-B函数生成的轮廓,计算误差在4%以内,在分形维数为1.3~1.6时,计算误差在1%以内。杨红平等[1]通过分形轮廓小波分解系数方差与分形维数间的关系求出分形维数,但没有明确给出通用性的小波分解层数。

鉴于此,本文将小波的多尺度分析能力与分形的多尺度自相似特征相结合,先对分形轮廓进行多尺度小波分解,基于对数小波谱,识别出轮廓具有分形特征的小波分解尺度(从对数小波谱上看在一条直线上或接近同一直线),进而评价轮廓的分形特征并计算其分形维数,提出了机械加工表面轮廓分形维数对数小波谱计算方法,并与功率谱密度法、结构函数法、均方根法以及方程组法的计算精度进行对比,最后将其应用于实际机械加工表面轮廓分形维数的计算,说明了其实用性。

1 机械加工表面轮廓分形维数的对数小波谱计算方法

小波被称为“数学显微镜”,在时域和频域都具有良好的局部化性质,有多尺度/多分辨的特点。小波分析是把待分析的信号分解为一系列由母小波函数ψ(t)经尺度伸缩和平移所得到的小波函数族ψkm(t)((m,k)∈Z2)的线性叠加,类似于傅里叶变换将信号分解为一系列不同频率的正弦函数的叠加。小波分析是将信号在一组正交基{ψkm(t),(m,k)∈Z2}中展开,当尺度按二进制变化时,母小波的尺度伸缩和平移可以表示为

ψkm(t)=2-m/2ψ(2-mt-k)(1)

式中,2-m为尺度参数;k为平移参数。

任意信号f(t)∈L2(R)(平方可积函数空间)的离散小波级数表示为

小波分解系数可以表示为

式中,符号“〈〉”表示内积运算;m为尺度序号;dmk为m层分解的第k个分解系数。

自相似过程的功率谱密度满足[16],更确切的功率谱的表达式为[17]

S(ω)=ωσ2xα(4)

其中,α=2 H+1,H为Holder指数,也称Hurst指数;σx2为自相似过程的方差;ω为自相似过程的频率(时间或者空间)。

自相似过程具有分形特征,根据文献[17]对自相似的定义f(t)=γ-Hf(γt),分形维数

序列f(t)在小波基{ψkm(t),(m,k)∈Z2}上以及序列f(τ)在小波基{ψk′m′(τ),(m′,k′)∈Z2}上的小波分解系数分别为

定义1将序列f(t)的自相关函数定义为

基于多尺度分析的方法,对具有自相似性的粗糙轮廓高度序列进行小波分解,在不同的尺度上,则其小波分解系数满足[18]

其中,E[]表示求均值,“*”表示卷积运算。将ψ(t)的傅里叶变换记为Ψ(ω),则ψkm(t)的傅里叶变换为

将式(8)右端利用Parseval定理以及时域卷积定理转化到频域后,式(8)可改写为

特别地,m=m′,k=k′时,有

从而有

式中,var[]为求方差运算。

对式(11)两边取以2为底对数,有

将式(13)左端记为Ym,得

对给定的具有分形特征的实际机械加工表面轮廓进行M层的小波分解,各分解尺度(层)对应的m(m=1,2,…,M)、Ym便是已知的,相应地,对数据点(m,Ym)序列进行直线拟合即可求出直线斜率α,再由前述关系α=2 H+1及D=2-H便求出分形维数D,即

在此称式(13)或式(14)为对数小波谱。以上分析推理求解机械加工表面轮廓分形维数的方法,本文称之为对数小波谱方法。

2 机械加工表面轮廓分形模拟与分形维数计算

M-B函数被广泛用于模拟实际机械加工表面轮廓,M-B函数是一个确定性自相似过程,满足Z(x)=γ-(2-D)Z(γx)。其具体表达式为

式中,z(x)为表面轮廓高度;G为分形特征长度尺度系数;nl为最低频率指数;γnl=1/L;L为采样区间长度。

在实际应用中,最高频率的取值是有限的,即

应用M-B函数对机械加工表面分形轮廓进行数值仿真,根据文献[19]中对实际机械加工表面的测量,轮廓参数为:采样间隔l=0.204μm,分形特征长度尺度系数G=2.86×10-10m,采样区间长度L=213l=1671μm,分形维数D取1.1、1.2、…、1.8、1.9。为满足采样定理,取最高频率指数nmax=log1γ2l。在相同分形维数和相同分形特征长度尺度系数的情况下,应用了3种不同采样区间轮廓,采样区间随机地分别取(0,L)、(3L,4L)和(10L,11L)。由此获得的模拟轮廓如图1所示(本文仅仅给出了D=1.5时对应各采样区间的M-B函数生成的轮廓)。

在对机械加工表面轮廓进行小波分解时,首先要考虑的是小波基函数的选择和分解尺度(层数)的确定。如何选择小波基函数,目前还没有一个理论标准,大多凭经验。工程中较多采用的小波基函数有db小波、sym小波和haar小波。对于给定长度的序列,最大分解尺度(层数)也是有限制的。从小波理论分析角度来看,最大分解尺度仅由数据长度和小波阶数便可确定,理论上小波分析最大分解尺度M为[20]

式中,floor(·)为向下取整函数;lx为信号长度;lw为滤波器长度,与小波类型有关,对于dbN小波,lw=2 N;N为小波函数的消失矩。

本文分别选用db2、sym2、db4、sym4以及haar小波基函数,由式(18)可确定db2与sym2小波对应的最大分解尺度为11,db4与sym4小波对应的最大分解尺度为10,采用haar小波时则为13。根据自相似特征确定具有分形特征的尺度,对这些尺度进行拟合就可以计算出分形维数。由M-B函数生成理论分形维数Dt为1.1、1.2、…、1.9的模拟轮廊,对上述轮廓在采样区间(3L,4L)上分别进行小波分解,应用本文提出的对数小波谱方法,对具有分形特征的分解尺度进行拟合获得直线斜率α,从而由式(15)计算出相应的分形维数D,结果见表1。

由表1可见,采用sym4小波时计算结果精度最高。图2给出了理论分形维数Dt分别为1.2、1.5、1.8时的对数小波谱及其拟合直线。根据图2,分解尺度m为2、3、4、5、6对应的对数小波谱几乎在一条直线上,其余理论维数类似,轮廓在这几个尺度上(对应轮廓高频部分)具有精确的自相似性,说明轮廓在这几个尺度上具有良好的分形特征(从对数小波谱上看在一条直线上,称之为有效分解尺度),因此分解尺度2、3、4、5、6为有效分解尺度;在尺度大于8(对应轮廓低频成分)时,与小于8的尺度没有严格的相似性,这是因为在生成模拟轮廓时采用的是式(17),而式(17)是有限项级数,存在误差,因此不具有严格的自相似性。对有效分解尺度2~6进行拟合计算的精度很高,误差都在0.15%以内。

应用对数小波谱方法,采用sym4小波计算获得的分形维数与应用功率谱密度法、结构函数法、均方根法以及解方程组法获得的结果的比较如图3所示。

由图3可见,功率谱密度法计算误差最大,在理论分形维数Dt=1.3时,识别误差最小,但Dt偏离1.3时误差急剧变大,Dt=1.1时,误差高达31.6332%;方程组法次之,最大误差也超过了10%;均方根法误差也较大;结构函数法计算精度很高,在Dt<1.3时计算结果大于理论分形维数,Dt>1.3时计算结果小于理论分形维数,Dt在1.3~1.6时计算结果基本接近理论分形维数,误差绝对值都在1%以内;本文提出的对数小波谱法采用db2小波与sym2小波时计算结果一样,误差也都在1%以内且比文献[1]方法的误差小,采用sym4小波时误差最小,小于0.15%。由于功率谱密度法、解方程组法和均方根法都是建立在分形轮廓的功率谱密度函数近似的连续表达式基础上,M-B函数的功率谱密度是其自相关函数的离散傅里叶变换,实际上是一系列冲击函数的和式[21],因此导致了利用这些方法计算分形维数时产生较大的误差;轮廓经离散傅里叶变换后的功率谱数据点多,振荡幅度大,不易于直线拟合,功率谱密度法对轮廓的采样点进行离散傅里叶变换得到的频谱和功率谱也是近似的,这也导致功率谱密度法的误差很大。由Majumdar等[4,13]引入的结构函数法的计算误差较小,这同文献[10,11]的结论一致;结构函数法直接通过轮廓高度信息计算维数,过程简单,避免了近似计算带来的误差,计算精度较高,其误差的绝对值在4%以内。由图2可见,用sym4小波计算M-B函数模拟轮廓分形维数时,得到的对数小波谱数据点直线度很好,且受最大分解层数限制,数据点也少,作直线拟合时易于实现,表1和图3结果表明其计算误差都在0.15%以内。需要指出的是,对于前面提到的采样区间(0,L)和(10L,11L),所得结论与上面一致,可见对数小波谱法能很好地识别出轮廓的分形维数,要优于其他方法。

1.对数小波谱方法(sym4小波)2.均方根法3.结构函数法4.文献[1]方法5.方程组法6.功率谱密度法7.分形维数理论值

将对数小波谱方法应用到文献[14]中的M-B函数算例,得到的计算结果与文献中其他方法计算结果的比较见表2。由表2可见,对数小波谱方法相对于其他方法有更高的精度。

3 计算实例

图4为由Talysurf5-120轮廓测试仪对车削、铣削和磨削加工表面进行测试获得的表面轮廓曲线,试件材料为45钢,测试的采样间隔为1.25μm,采样长度为3.75mm,共有3000个离散采样点。对该轮廓进行db2小波分解,得到图5所示的各加工表面轮廓的对数小波谱,经直线拟合求出其分形维数。分别采用功率谱密度法与结构函数法对上述各轮廓分形维数进行计算,结果如表3所示。

4 结论

(1)提出了一种机械加工表面轮廓分形维数的对数小波谱计算方法以及有效分解尺度的概念。

(2)计算轮廓分形维数的关键是选择小波基函数和分解尺度,然而其最大分解尺度可由表面轮廓序列长度和小波阶数确定,可先按最大分解尺度分解。

(3)对数小波谱法能很好地处理分形的多尺度特征,对于由M-B函数模拟的分形轮廓,特别是采用sym4小波时,对数小波谱方法计算误差在0.15%以内,对文献[14]中的算例计算也具有高的准确度。