Numerical Wear Analysis of a Three-dimensional Rough Surface

수치적 방법을 이용한 3차원 거친 표면의 마모 해석

Abstract

It is essential to predict the amount of wear and surface parameters for a surface where relative motion occurs. In the asperity-based model for wear prediction, only the average contact pressure can be obtained. Hence, the accuracy of wear analysis is poor. In this study, DC-FFT is used to obtain the pressure of each node, and wear analysis is performed by considering the effect of the pressure gradient. The numerical surface generation method is used to create Gaussian, negatively skewed, and positively skewed surfaces for wear analysis. The spatial and height distributions of each surface are analyzed to confirm the effectiveness of the generated surface. Furthermore, wear analysis is performed using DC-FFT and Archard's wear formula. After analysis, it is confirmed that all peaks are removed and only valleys remain on the surface. The RMS roughness and Sk continue to decrease and Ku increases as the cycle progresses. It is observed that the surface parameters are significantly affected by the radius of curvature of the asperity. This analysis method is more accurate than the existing average wear and truncation models because the change in asperity shape during the wear process is reflected in detail.

1. 서론

마모는 접촉하는 두 표면 간의 상대 운동으로 인해 표면의 일부가 닳아 재료가 제거되는 현상으로, 복잡한 물리적, 화학적 현상의 조합에 의해 일어난다. 마모는 재료의 수명 단축 뿐 만 아니라 표면 거칠기를 변화시키며 이외에도 진동, 소음, 마찰 등 기계 시스템의 많은 측면에 영향을 주기 때문에, 마모에 대한 연구는 공학적, 경제적 측면에서 아주 중요하다. 하지만 마모의 원리에 대한 연구가 수 세기 동안 진행되어 왔음에도 불구하고, 아직 마모에 대한 이해는 완전하지 않다. 

마모가 일어나는 원리는 아주 복잡하기 때문에 실험적 방법으로 마모 계수를 도출할 수 있다. 하지만, 마모 실험은 결과의 산포가 매우 크고, 온도, 습도 등 환경적인 조건에도 민감하게 바뀌어 마모 실험에서 표면 거칠기까지 다루는 것은 상당히 어려운 작업이다. 이러한 이유로 마모 실험은 많은 시간이 소요될 뿐 아니라 진행과정 동안의 정보를 알 수 없기 때문에 마모를 유발하는 주된 파라미터를 규명하기에는 무리가 있다. 따라서, 수치적으로 표면의 마모를 해석하는 것은 초기 표면이 마모에 미치는 영향과 마모 진행에 따른 표면의 변화를 예측하기에 적절한 방법이라 할 수 있다.

표면의 마모에 대한 연구에서, 마모가 진행되는 동안 표면 파라미터들이 크게 변화하는 것이 관찰되었다. Hutt[1]는 mixed-EHL에서 작동하는 기어의 길들임(running-in)에 대해 연구하기 위하여 두 디스크의 마모 실험을 하였고, 길들임 과정에서 주요 돌기들이 대부분 평평해짐을 확인하였다. Clarke[2]또한 초기의 몇 사이클에서는 돌기가 소성 변형한다는 사실을 실험들을 통해 확인하였다. 이러한 길들임 과정 동안의 표면의 마모 현상에 대한 해석적 연구는 King과 Spedding[3]의 truncation model이 대표적이다. 이들은 비 정규 표면의 프로파일이 특정한 면에 의해 잘리는 것으로 마모를 해석하는truncation 모델을 제안하였다. Ghosh와 Sadeghi[4]가 마모 해석에 사용한 평균 마모 해석은 King과 Spedding의 truncation 모델과 아주 비슷하다. 이들은 넓은 범위의 표준 편차, 첨도(Skewness), 왜도(Kurtosis)를 갖는 다양한 거친 표면에서 일어나는 마모를 해석하였다. 2D 표면 프로파일에 대해 Hertzian 접촉 이론 및 Archard의 마모식을 사용하여 마모량 및 표면 파라미터를 분석하였으며, 마모 실험 결과와 해석 결과를 비교하여 모델을 검증하였다. 하지만 표면의 특성이 2D 프로파일에서는 완전히 표현되지 않기 때문에[5] 3D 표면 데이터에 대한 마모 해석이 필요하다. Prajapati와 Tiwari[6]는 수치적인 3D 마모 모델을 개발하여 마모 과정 중에서 일어나는 표면 파라미터들의 변화 경향을 살펴보았다. 그들은 사이 클 수에 따른 표면 파라미터들의 변화를 관찰하였고 시편과 동일한 특성의 표면을 수치적 생성, 마모 시뮬레이션 한 결과 실험 결과의 경향과 잘 일치하였다.

하지만 지금까지의 거친 표면의 마모에 대한 연구는, Greenwood-Williamson(GW)의 모델을 기반으로 각 돌기의 곡률을 상당하여 단순한 방법으로 마모를 계산하는 방법이다. 하지만, 접촉 압력은 돌기의 형상에 따라 크게 다르며 접촉 압력이 마모량을 지배하는 주 요소이기 때문에, 마모 진행에 따른 접촉 면적의 변화, 최종 돌기 형상 및 마모 후의 표면 파라미터 등을 면밀히 관찰하기 위해서는 압력 계산 과정에서 돌기 형상 자체를 그대로 반영할 수 있어야 한다. 따라서 본 논문에서는 각 격자점의 압력을 계산하여 마모 해석을 수행하였으며, 시간이 많이 소요되는 접촉 압력 해석은 Discrete Convolution Fast Fourier Transform(DC-FFT)를 사용하여 해석의 속도를 향상하였다. 마모가 일어나는 동안의 돌기의 형상 및 표면 파라미터들의 변화를 살펴보았다.

2. 연구방법 및 내용

2-1. 거친 표면의 수치적 생성

고차원적인 거친 표면의 마모 해석을 수행하기 위해서는 먼저 원하는 거친 표면을 정확히 만들 수 있어야 한다. 거친 표면의 정보는 표면 위의 높이 방향의 분포를 의미하는 높이 분포 함수(Height Distribution Function, HDF)과 거친 표면 높이의 표면 위에서 연관성을 나타내는 자기 상관 함수(Auto correlation function, ACF)로 표현할 수 있다. 이에 따라 거친 표면의 수치적 생성은 목표로 하는 HDF와 ACF를 가진 표면을 구현하는 것이다. 하지만, 수치적으로 생성된 표면의 두 함수는 목표로 하였던 함수에 완전히 일치할 수 없기 때문에, 수치적 거친 표면의 생성의 주안점은 적합한 HDF와 ACF를가지는 표면을 만들어 내는 것이다.

거친 표면 생성에서 중요한 두 함수에 대해서 살펴보자. ACF는 x − y 평면의 표면 높이 프로파일 z(x, y)의 공간적 관계와 의존성에 대한 정보를 제공하며 다음과 같이 표현된다. 

\(R_{zz}\left(\Delta x_{\mathrm{n}} \Delta y\right)=E\{Z(x, y) z(x+\Delta x, y+\Delta y)\}\)       (1)

여기서, E는 기대값 연산자(expectancy operator)이며, 식 (1)의 ACF를 N × M개의 격자점으로 이산화(discrete) 하여 표현하면 다음과 같다.

\(​​​​R_{z=}(p, q)=\frac{\sum_{i=0}^{N-k} \sum_{j-0}^{M-1} z(i, j) z(i+p, j+q)}{s_{q}^{2}(N-p)(M-q)}\)       (2)

여기서, i, j은 임의의 위치, p, q는 i, j로부터 떨어진 임의의 거리를 나타낸다. 자기 상관 함수를 대표하는 파라미터로는 자기상관길이(Auto correlation length)가 있는데, 이 값은 표면을 표현하는 유용한 정보로 ACF가 임의의 값(일반적으로 0.1)을 가지는 지점의 길이이다. 참고로, 선삭, 밀링 등에 의해 가공된 표면과 같은 공학적 표면의 ACF는 대부분 지수 함수 형태이다[7]. HDF는 명칭에서 의미하는 바와 같이 높이에 대한 정보를 나타내며, 이 함수로부터 거친 표면을 대표하는 다양한 통계적 파라미터가 계산된다. 표면의 높이가 이산적인 데이터 \(\left\{z_{i, j} \mid i=0,1,2, \ldots, N-1, j=0,1,2, \ldots, M-1\right\}\)로 주어졌을 때 4가지 중심 적률(central moments)은 각각 평균 (average) s_a, 표준편차(root mean square) s_q, 왜도(skewness)s_sk, 첨도(kurtosis) s_ku라고 하며 다음 식과 같다.

\(s_{a}=\frac{\sum_{i-1}^{N} \sum_{j^{-1}}^{M} z_{i j}}{N M}\)       (3)

\(s_{q}=\frac{\sum_{i-1}^{N} \sum_{j^{-1}}^{M} z_{i j}^{2}}{N M}\)       (4)

\(s_{s k}=\frac{\sum_{i-1}^{N} \sum_{j}^{M} z_{i j}^{3}}{s_{q}^{3} N M}\)       (5)

\(s_{k u}=\frac{\sum_{i-1}^{N} \sum_{j^{-1}}^{M} z_{i, j}^{4}}{s_{q}^{4} N M}\)       (6)

높이 분포는 s_sk = 0이고 s_ku = 3인 정규(Gaussian) 분포와 그 외의 비정규(Non-Gaussian) 분포로 구분한다.

높이 분포 함수와 자기상관함수가 결정되고 나면 표면을 생성할 수 있다. 거친 표면을 생성하기 위한 모델에는 Moving Average(MA) 모델[8], Auto-Regressive(AR) 모델[9], Auto Regressive Moving Average(ARMA) 모델 등이 있다[10]. AR, ARMA모델의 경우 낮은 차수를 가진 시스템이기 때문에 큰 자기상관길이에서 불안정하다. MA모델은 정밀도가 높지만 자기상관길이가 길어질수록 비선형 방정식 행렬의 커져 계산 정확도 및 효율이 떨어진다. 본 연구에서는 MA모델과 수학적으로 일치하는 유한 임펄스 응답 필터(Finite impulse responsefilter, FIR)를 사용하고 FFT(Fast Fourier Transform)를 활용하여 계산 속도를 향상 시킨 Hu와 Tonder[11]의 기법을 활용하였다. 지금부터는 이들의 방법을 활용하여 거친 표면을 생성하는 방법을 소개하고자 한다. 먼저 정규 분포 형태의 거친 표면 생성에 대해 설명한 뒤, 비 정규 분포 형태의 거친 표면 생성에 대해서도 다룰 것이다. N × M개의 점을 갖는 정규 분포 표면의 높이는 다음 식으로 구할 수 있다. 

\(\begin{array}{l} Z_{i, j}=\sum_{k=0}^{n-1} \sum_{l=0}^{m-1} h_{k, l} \eta_{i+k, j+l} \\ i=0,1,2, \ldots, N-1, j=0,1,2, \ldots M-1 \end{array}\)       (7)

여기서, h는 FIR 필터 계수이며, η는 정규 분포를 갖는 독립적인 임의의 수의 행렬이다. ACF는 [n × m]의 크기를 가지는데, n, m은 절단 길이(truncation length)로 정의되며 Hu와 Tonder의 기법에 따라 n = N/2, m = M/2로 설정된다. 

식 (7)은 합성 곱 형태이기 때문에, FFT를 사용하면 합성 곱에 비해 계산 속도를 증가시킬 수 있다. FFT 계산을 위해서는 R_zz를 L_x × L_y의 가상의 도메인으로 확장하여 R_f를 새롭게 정의해야 한다. 주파수 영역에서 FFT를 수행하기 위한 L_x, L_y의 크기의 제약 조건은 각각 N + n −1, M + m −1보다 커야 한다. R_zz은 식 (8)에 의해 R_f로 확장된다.

\(\begin{aligned} &R_{j}\left(L_{x}-i, L_{y}\right)=R_{f}\left(L_{x}, L_{y}-j\right)=R_{f}\left(L_{x}-i, L_{y}-j\right)=R_{j}(i, j)\\ &i=0,1, \ldots, \frac{n}{2}, \quad j=0,1, \ldots, \frac{m}{2} \end{aligned}\)      (8)

R_zz와 마찬가지로, η 역시 가상의 도메인의 크기 L_x × L_y만큼 확장하여 η_f를 정의해야 하며, 그 과정은 식 (9)와 같다.

\(\begin{array}{l} \eta_{f}(i, j)=\eta(i, j) \\ i=0,1,2, \ldots, N-1, j=0,1,2, \ldots, M-1 \\ \eta_f(i, j)=0 \\ i=N, N+1, \ldots, L_{x}, \quad j=M, M+1, \ldots, L_{y} \end{array}\)         (9)

식 (7)을 L_x × L_y영역으로 확장한 FT는 다음과 같다.

\(\hat{\hat{z}}_{f}=\hat{\hat{h}}_{f} \cdot \hat{\hat{\eta}}_{f} \)       (10)

z_f는 L_x × L_y의 크기를 갖는 확장된 z 행렬로, N × M 영역에서 의미 있는 높이 데이터를 갖는다. \(\hat{\hat{z}}_{f}, \hat{\hat{\eta}}_{f}\), 는 z_f, η_f 의 푸리에 변환(Fourier Transform, FT)이다. 식 (10)에 의해 다음과 같은 관계식이 성립한다.

\(\hat{\hat{z}}_{f}^{2}=\left|\hat{h}_{f}\right|^{2} \hat{\hat{\eta}}_{f}^{2}\)       (11)

여기서 ,\(\hat{\hat{h}}_{f} \)는 h_f 의 FT이며, η_f 는 독립적인 임의의 수의 행렬이기때문에 η_f 의파워스펙트럼밀도(Power Spectrum Density, PSD)인 는 상수이다. 또한, = FFT(R_f)이므로, 식 (11)은 식 (12)와 같이 나타낼 수 있다.

\(\hat{\hat{h}}_{f}=\left(\frac{\left|F F T\left(R_{f}\right)\right|}{\hat{\hat{\eta}}_{f}}\right)^{\frac{1}{2}} \)       (12)

식 (12)를 통해 \(\hat{\hat{h}}_{f}\)를 구한 뒤, 식(10)을 역 푸리에 변환(Inverse Fast Fourier Transform, IFFT)하면 z_f 를 구할 수 있다. z_f 의 N × M영역에 존재하는 값이 설정한 ACF를 갖는 거친 표면의 높이 z이다.

여기까지 정규 분포를 갖는 거친 표면을 생성하는 방법에 대하여 설명하였다. 하지만 대부분의 공학적 표면은 비정규 분포이기 때문에[12], 비 정규 분포를 가지는 거친 표면 또한 생성할 수 있어야 한다[13,14]. 비 정규 분포를 가진 표면을 만들기 위해서는 높이 분포에 대한 정보인 η를 변환시켜야 하며, 이 과정에는 Johnson tran-slation system이 사용된다. Johnson translation system은 일반적으로 비 정규 분포의 확률 과정(stochastic process)를 정규 분포의 확률 과정으로 변환하기 위해 쓰인다. 원하는s_sk,z, s_ku,z를 갖는 거친 표면을 생성하기 위해서는, 아래의 관계식을 통해 s_sk,η, s_ku,η을 계산하여 비정규 분포의 임의의 수의 행렬인 η^'를 생성해야 한다.

\(s_{s k, z}=\frac{\sum_{i-0}^{q} \theta_{i}^{3}}{\left(\sum_{-0}^{q} \theta_{i}^{2}\right)^{3 / 2}} s_{s k, \eta}\)       (13)

\(s_{k w, z}=\frac{\sum_{-0}^{q} \theta_{i}^{4} s_{k u, \eta}+6 \sum_{r=0}^{q-1} \sum_{i^{\prime} i+1} \theta_{i}^{2} \theta_{j}^{2}}{\left(\sum_{-0}^{q} \theta_{i}^{2}\right)^{2}}\)       (14)

식 (13)과 (14)에서 \(\theta_{i}=h_{k, l}, i=k m+l, \quad q=n m-1\)이다. Johnson translation system 변환은 3가지 함수 곡선인 Unbounded, Log-normal(η' > ξ), Bounded(η^' > ξ + λ)로 구성되어 있다.

\(\text { Unbounded: } \eta=\gamma_{J}+\delta_{J} \sinh ^{-1} \frac{\eta^{\prime}-\xi_{d}}{\lambda_{J}}\)       (15)

\(\text { Log-normal: } \eta=\gamma_{J}+\delta_{J} \log \frac{\eta^{\prime}-\xi_{J}}{\lambda_{J}}\)       (16)

 \(\text { Bounded: } \quad \eta=\gamma_{J}+\delta_{J} \log \frac{\eta^{\prime}-\xi_{J}}{\xi_{J}+\lambda_{J}-\eta^{\prime}}\)       (17)

여기서, η는 정규 분포 데이터이며, η^'는 비 정규 분포 데이터이다. 주어진 조건에 따른 4가지 중심 적률에 대한 γJ, δJ, λJ, ξJ는 Hill의 알고리즘[15]을 사용하여 구할 수 있다. \(s_{k u}-s_{s k}^{2}-1>0 \)을 만족해야 Johnson translation system을 사용 할 수 있기 때문에 특정한 조건에서는 변환이 불가능하다는 제약이 있다[16]. 결과적으로 η^'을 구
한 뒤 정규 분포와 같은 방식으로 비 정규 분포 높이를 얻을 수 있다. 표면을 생성하는 방법을 요약하면 다음과 같다.

A1) 생성하려는 표면의 특성에 따라 적절한 ACF를 선 택한다. 이때, ACF 행렬 Rzz의 크기는 [n × m]이다.

A2) 조건에 따른 η를 생성한다. 정규 분포의 경우, 독립적인 임의의 수의 행렬을 생성하며, 비 정규 분포의 경우 식 및 Johnson translation system을 사용하여 요구되는 s_sk,η, s_ku,η를 가지는 η^'를 생성한다. 행렬 η의 크기는 [N × M]이다.

A3) FFT를 하기 위해 식 (8), (9)를 사용하여 R_f, η_f를생성한다.

A4) 식 (12)를 통해 필터 계수의 FFT인 \(\hat{\hat{h}}\)를 구한다.

A5) 식 (10)에 대한 IFFT를 계산하여 특정한 ACF 및 s_sk,z, s_ku,z를 갖는 표면 z를 생성한다.

2-2. 거친 표면의 접촉 압력 해석

거친 표면의 접촉 압력을 해석하는 방법은 크게 두 가지로 나뉘는데, 전자는 거친 표면의 돌기들의 곡률을 특정 등가 곡률로 상당한 후 각 돌기들을 Hertzian 접촉으로 가정하여 계산하는 방법이며, 후자는 대상 영역의 격자계 위에 거친 표면의 높이 분포를 있는 그대로 반영하여 접촉 압력을 해석하는 방법이다. 전자의 가장 대표적인 방법은 GW 모델[17]로 해석의 편리함에 의해 현재도 널리 활용되고 있다. GW 이론에 따르면 한 돌기의 접촉 면적에 대한 평균 압력은 계산이 가능하지만, 실제 접촉 압력은 접촉 면적에 대해 포물선 형태를 가지기 때문에 한 돌기 내의 압력 구배에 대한 정보는 얻을 수 없다. 압력 구배는 마모 형상 및 깊이에 영향을 미치기 때문에 더욱 정밀한 압력 분포 데이터가 필요하다. 본 연구에서는 거친 표면이 마모되는 과정을 보다 정확히 살펴보기 위해 후자의 방법을 채택하였다. 거친 표면을 모사하기 위해서는 많은 노드 수가 필요하기 때문에, 거친 표면의 접촉 해석은 많은 시간이 소요된다. 거친 표면의접촉을 수치적으로 해결하기 위해서 Conjugate Gradient Method[18], Multi-Level Multi-Integration technique[19], FFT[20] 등 다양한 효율적인 방법들이 연구되었다. 그 중에서Chen등[21,22]에 의해 제안된 Discrete Convolution (DC)-FFT방법이 단순하면서도 빠른 해석 속도를 가지는 장점으로 현재 가장 널리 활용되고 있다. 그러므로 본 연구 에서도 DC-FFT를 이용하여 접촉 해석을 수행하였으며, 그 방법을 본 절에서 간단히 살펴보고자 한다. 이 해석에 대한 보다 자세한 설명은 Wang 등 연구[21,22]에매우 잘 기술되어 있다.

균질한 반무한 표면의 수직 하중에 관련된Boussinesq 문제에 대해 수직 변형에 관한 Green함수 G는 식 (18)과 같고, 이에 따라 변형량 u는 식 (19)와 같이 정의된다.

\(G(x, y)=\frac{1}{\pi E^{\prime} \sqrt{x^{2}+y^{2}}}\)       (18)

\(u(x, y)=\left(\int_{-\infty}^{\infty} \int_{-\infty}^{\infty}\right) G\left(x-x^{\prime}, y-y^{\prime}\right) \cdot p\left(x^{\prime}, y^{\prime}\right) d x^{\prime} d y^{\prime}\)       (19)

여기서, \(E^{\prime}=\left(1-v_{1}^{2}\right) / E_{1}+\left(1-v_{2}^{2}\right) / E_{2}\) 상당 탄성 계수이고, E₁, E₂, v₁, v₂는 각각 상대 접촉하는 두 물체의 탄성 계수와 푸아송 비이다. Green함수는 접촉 문제에서 영향 함수(Influence coefficient)로 표현된다. 접촉 표면을 N_x × N_y의 유한개의 격자로 구성되어 있다고 가정하고, 각 요소의 크기는 2Δ_x × 2Δ_y로 두자. 영향 함수(Influence coefficient, IC) D_{α, β} 는 원점에서 단위 압력이 작용할 때 (p_0,0 = 1), (2αΔ_x, 2βΔ_y) 위치에서의 변형을 의미하며, 식(20)과 같이 나타낼 수 있다.

\(D_{\alpha, \beta}=\int_{-\Delta_{x}}^{\Delta_{x}} \int_{-\Delta_{y}}^{\Delta_{y}} G\left(2 \alpha \Delta_{x}-x^{\prime}, 2 \beta \Delta_{y}-y^{\prime}\right) d x^{\prime} d y^{\prime}\)       (20)

여기서, Green 함수의 적분은 다음과 같다.

\(\begin{aligned} \iint G(x, y) d x d y &=\frac{1}{\pi E^{\prime}}\left[x \ln \left(y+\sqrt{x^{2}+y^{2}}\right)\right.\\ &\left.+y \ln \left(x+\sqrt{x^{2}+y^{2}}\right)\right] \end{aligned}\)       (21)

식 (21)를 IC를 이용하여 DC 형태로 다시 쓰면 다음과 같다.

\(u_{\alpha, \beta}=\sum_{\xi=0}^{N_{x}-1} \sum_{\psi=0}^{N_{y}-1} D_{\alpha-\xi, \beta-\psi p \xi, \psi}\)       (22)

여기서 \(0 \leq \alpha, \xi이다. FT의 연속 합성 곱 이론에 따라, 공간 영역의 합성 곱은 주파수 도메인의 직접 곱이 된다. 식 (22)에 2차원 FT를 적용하여 주파수 영역에서의 변형량을 계산하면 식 (23)과 같다.

\(\widehat{\hat{u}}_{\alpha, \beta}=\widehat{\bar{D}}_{\alpha, \beta} \hat{\hat{p}}_{\alpha, \beta}\)       (23)

여기까지는 전체 표면에서 압력에 따른 변형량을 설명 하였다. 접촉 문제를 풀기 위해서는 다음의 구속 조건에 의한 반복 연산이 필요하다. 접촉 영역 Ω 안에서는

접촉 영역 Ω 안에서는

\(p(x, y)>0, \quad g(x, y)=0\)       (24)

접촉 영역 Ω 밖에서는

\(p(x, y)=0, \quad g(x, y)>0\)       (25)

이다. 여기서, g(x, y)는 변형 후 표면 형상의 높이 분포로 다음과 같이 나타낼 수 있다.

\(g(x, y)=z(x, y)+u(x, y)-\delta\)       (26)

여기서, δ은 유효 강체 변위(effective rigid body displacement)이며, z(x, y)는 변형되기 전의 거친 표면의 높이 분포이다. DC-FFT 방법에서는 Conjugate gradient(CG) 알고리즘을 이용하여 압력을 반복 계산한다[22]. 접촉 표면의 압력 분포는 거친 표면의 높이 분포 z와 유효 강체변위 δ에 의해 결정된다. 이 압력 분포는 특정 변위에 대한 값으로 접촉 표면에 가해지는 수직 하중 FN에 대해 식 (27)과 같이 힘의 평형을 만족하는 δ를 반복 계산하여 전체 표면의 압력 분포를 획득한다.

\(F_{N}=\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} p(x, y) d x d y\)       (27)

2-3. 마모량 계산

2-2절의 방법을 이용하여 이산화 된 거친 표면의 각 격자점 위치에서의 접촉 압력을 계산할 수 있다. 이 때, 각 격자점에서는 접촉 압력이 직사각형 미소 표면(dxdy) 위에 균일하게 분포되어 있다고 가정하였다. 각 미소 표면에서 일정 압력 p가 작용했을 때, s만큼 미끄러 진 후의 마모량은 Archard 마모식[23]으로 계산할 수 있다.

\(h_{w}=\frac{k p s}{H}\)       (28)

여기서, k는Archard 마모 계수이고, H는 거친 표면의 경도이다.

2-4. 거친 표면의 미끄럼 마모 해석 과정

거친 표면의 미끄럼에 따른 마모를 예측하는 것은 매우 어려운 작업이며, 수치 해석적으로 미끄럼 마모를 완벽히 모사하는 것은 불가능하다. 그러므로, 미끄럼 마모해석을 수행하기 위해 본 연구에서는 다음과 같은 가정으로 마모 모델을 단순화하였다.

B1) 거친 표면에는 오직 탄성 변형만 일어나며, 표면의 접촉 압력 역시 탄성 역학에 기초하여 계산되었다.

B2) 미끄러 지는 순간마다 두 거친 표면은 완벽한 대칭 형태인 것으로 가정하였다. 즉, 두 거친 표면이 미끄러지면서 상대적 접촉 위치가 변화하는 것은 고려하지 않았다. 이 가정을 바탕으로 두 거친 표면 사이의 접촉을 하나의 등가 된 거친 표면과 완벽히 매끄러운 강체 사이의 접촉으로 치환할 수 있다.

B3) Archard마모 계수 k는 미끄럼 상황에 영향을 받지 않는다. 그러므로 전체 해석 구간에 대해서 k는 일정하다고 가정한다.

마모 해석은 접촉 평형 상태에 도달할 때까지 진행한다. 여기서, 접촉 평형은 일반적으로 겉보기 면적에 대한 실 접촉 면적의 비 m이며, 일반적으로 m = 0.75를 접촉 평형의 기준으로 삼는다[4]. Fig. 1에 마모 해석의 흐름도를 나타내었으며, 상세한 내용은 다음과 같다.

Fig. 1. Flow chart for wear simulation.

C1) 표면의 간극 분포와 작용 하중으로부터 DC-FFT 방법을 활용하여 압력 분포를 계산한다.

C2) 각 표면 위의 격자점에 대하여, 마모량을 계산하여 새로운 마모 표면으로 교체한다.

C3) 수렴 조건을 만족하였는지 체크하고, 아닐 경우 마모된 새로운 표면을 가지고 C1)의 과정으로 되돌아간다.

3. 결과 및 고찰

3-1. 표면 생성

거친 표면의 마모 해석을 수행하기 위하여, 서로 다른 표면 파라미터를 가지는 3가지 거친 표면을 2.1절의 방법을 활용하여 수치적으로 생성하였다. 생성된 세 표면의ACF는 지수 함수 형태로 식 (29)와 같으며, 각 표면의 수치 생성에 관련된 변수들은 Table 1에 정리하였다.

Table 1. Input parameters used in generating surface

\(R_{z=}(p, q)=\exp \left[-2.3\left\{\left(\frac{\beta_{x}}{\beta_{x}^{*}}\right)^{2}+\left(\frac{\beta_{y}}{\beta_{y}^{*}}\right)^{2}\right\}\right]\)       (29)

여기서, β ^*은 자기상관길이 이며 ACF가 0.1이 되는 길이로 정의하였고, \(\beta_{x}^{\prime}=\beta_{y}^{\prime}\)로 등방성 표면의 자기 상관함수이다. 정규 분포를 갖는 표면(Gaussian surface)과 2개의 비 정규 표면을 갖는 표면(Negatively skewed surface, Positively skewed surface)을 생성하였으며, 생성된 세 표면의 형상을 Fig. 2에 나타냈다. 생성된 표면이 목표하는 ACF 및 높이 분포를 만족하는 지 확인하기 위해 이론적 ACF와 표면의 x방향 프로파일의 ACF를 비교하였으며, 이를 Fig. 3에 나타내었다. 생성된 표면의 ACF는 채택된 이론적 ACF 함수와 유사하지만, 오차가 발생하여 정확성이 떨어진다. 시스템의 오차는 잘림 길이와 자기 상관 길이와의 관계의 영향을 받는데, 잘림 길이는 자기 상관 길이의 약 6배 이상이어야 안정적이며 정확한 표면 파라미터를 얻을 수 있다[25]. 본 연구에서는 두 파라미터가 4배의 비율을 가져 ACF의 오차가 존재한다. Gaussian surface와 Negatively skewed surface에서는 β^*= 16에서 ACF가 0.1에 가까운 값을 가지며 Positively skewed surface의 경우 약간의 오차가 있음을 확인할 수 있다. Fig. 4의 표면 높이의 분포를 나타내는 확률 밀도 함수(Probability Density Function, PDF)에서 세 표면의 s_sk,z, s_ku,z에서 기인한 차이를 확인할 수 있다. Fig. 4(a)의 Gaussian surface의 PDF는 굵은 실선으로 그려진 정규 분포 함수와 일치한다. Negatively skewed surface의 경우 Fig. 4(b)에서 보이듯 분포가 왼쪽으로 치우쳐져 있으며, Positively skewed surface는 오른쪽으로 치우쳐져 있고 다른 두 표면에 비해 sku,z가 커 분포의 피크가 뾰족하다.

Fig. 2. Generated surface. (a) Gaussian, (b) Negatively skewed, (c) Positively skewed.

Fig. 3. Autocorrelation function with \(\boldsymbol{\beta}^{*}=\boldsymbol{\beta}_{x}^{*}=\boldsymbol{\beta}_{y}^{*}=\mathbf{1 6}\). (a) Gaussian surface, (b) Negatively skewed surface, (c) Positively skewed surface.

Fig. 4. Probability density function of surface height. (a) Gaussian surface, (b) Negatively skewed surface, (c) Positively skewed surface.

3-2. 마모 진행 과정 해석

해석에 사용된 마모되는 두 물체는 Prajapati의 연구[6]에서 제시된 것과 같은 탄소강으로 가정하였으며, 재료 물성은 Table 2에서 확인할 수 있다. 수치적으로 생성된 거친 표면은 128 μm × 128 μm의 작은 면적을 가지며, 여기에 50 MPa(수직 하중 F_N은 0.82N)의 일정한 접촉 압력 p*가 가해지는 상태에 대한 마모 해석을 수행하였다. 이 때, 마모 계수 k역시 Prajapati의 연구[6]에서 사용된 값이다. 지속해서 미끄러지며 마모되는 두 표면은 일정 거리 동안 같은 표면을 가진다고 가정하였는데, 그 미끄럼 간격 s를 10 μm로 가정하였다. 이 가정의 의미는 s만큼 미끄러지는 동안에는 일정한 거친 표면을 가진다는 것이다. 작은 면적에 대해 해석을 진행한 이유는 돌기를 좀 더 자세히 모사하기 위해서이며, 거친 표면이 주기적으로 반복된다고 보면, 면압이 일정할 때는 해석 면적에 배수를 곱하여 보다 넓은 면적에 대한 마모 해석으로 확장해서 볼 수 있다.

Table 2. Operating condition

3-2-1. 평균 접촉 압력과 마모 깊이

돌기 기반 모델에서는 평균 압력으로 마모 깊이를 계산하기 때문에 접촉한 모든 격자의 높이가 같다. 따라서 실제 현상보다 돌기 곡률 반경이 급격히 증가하며 이는 해석에서 일어나는 오차의 원인이 된다. 엄밀한 해석을 하기 위해 본 논문에서는 DC-FFT로 계산한 압력을 사용하여 각 격자의 마모 깊이를 계산하였다. 두 방법에 따른 차이를 돌기 형상 변화를 통해 확인해보자. Fig. 5는 정규 분포를 갖는 표면의 x축 프로파일 중 일부분으로,굵은 실선은 초기 형상을, 굵은 파선은 마모 과정이 끝나고 접촉 평형을 이루었을 때의 형상이다. 기존 형상인 굵은 실선을 기준으로 40 마모 사이클회마다 같은 위치의 프로파일을 그려 돌기의 변화를 나타냈다. Fig. 5(a)는 접촉 영역 내 압력 구배를 사용하여 각 격자점에 대한 마모 깊이를 계산하였을 때의 마모 과정을 나타낸다.돌기의 형상이 점차 평평해지며 최종 돌기 형상이 약간의 곡선을 그린다. Fig. 5(c)는 평균 압력으로 계산된 마모 깊이에 따른 마모 형상이다. DC-FFT로 획득한 압력의 평균을 사용하여 돌기 기반 모델과 유사하게 나타나도록 하였다. 해석 초기부터 돌기 끝이 평평하며, 압력의 평균값을 사용했기 때문에 마모 깊이가 실제보다 얕게 도출되었다. 또한, 돌기 곡률 반경이 증가하여 접촉 면적이 커지기 때문에 다음 사이클에서 접촉 압력이 작게 계산되고, 마모 깊이가 작아진다.

Fig. 5. (a) Evolution of surface with DC-FFT method and (b) pressure at the last cycle, (c) evolution of surface with Truncation model and (d) pressure at the last cycle.

평균 접촉 압력(P^*= F_N/A_real)은 작용 힘을 실 접촉 면적으로 나눈 것이다. 접촉하고 있는 돌기 수가 많아질수록, 돌기가 마모되어 돌기 곡률이 커질수록 돌기의 접촉 압력이 작아진다. 따라서 접촉 압력은 표면의 형상에 따라 크게 달라지는데, 이는 Fig. 6에서도 나타난다. Gaussian surface와 비교하였을 때, 접촉하는 방향에 끝이 뾰족한 산이 많은 형상을 갖는 Positively skewed surface는 초기 평균 압력이 2.8GPa에 달하는데, 이는 s_ku.z가 3보다 커 돌기가 뾰족하여 실 접촉 면적이 작기 때문이다. 압력에 따라 돌기는 소성, 탄소성, 탄성 변형을 하게 되는데, 마모 초기에는 압력이 매우 크기 때문에 소성변형을 한다고 알려져 있다[24]. 하지만 이번 연구에서는 가정에 따라 소성에 의한 효과는 고려되지 않았다. Negatively skewed surface에서는 상대적으로 낮은 1 PGa를 초기 압력으로 가지며 사이클이 진행되는 동안 세 표면 중 가장 낮은 압력을 보이는데, 이는 s_ku,z의 영향으로 골이 많은표면 형상에 의해 접촉 면적이 크기 때문이다.

마모 깊이는 마모되기 전 표면의 최고점에서 마모 과정 중 표면의 최고점까지의 거리를 나타낸 것이다. Fig. 7의 마모 깊이 곡선은 크게 두가지 구간으로 나눌 수 있다. 표면의 종류에 따라 사이클 시작 지점부터 약 50~80번 까지의 구간에서는 마모 깊이가 급격한 증가세를 보이는데 이는 돌기의 평탄화로 인한 효과이다. 두 구간의 경계점은 Fig. 6의 평균 압력 곡선의 기울기가 급격하게 바뀌는 지점과 유사하며 이 지점에서 대부분의 돌기 산이 없어져 접촉이 안정해졌다고 볼 수 있다. 표면이 평평해진 이후부터는 사이클 당 마모 깊이가 일정한 값에 도달하며, 따라서 마모 깊이 그래프가 선형의 형태에 가까워지는 것을 알 수 있다.

Fig. 6. Mean pressure (P^* = FN /A_real).

Fig. 7. Wear depth.

3-2-2. 표면 거칠기

평균 거칠기 s_a와 RMS 거칠기 s_q는 마모가 진행될수록 감소한다. 돌기가 갈리면서 피크 점(산)의 높이가 감소하여 산과 골의 거리가 짧아지기 때문이다. 식 (28)을 면적에 대해 적분했을 때, 사이클 당 마모되는 부피는 같으나 마모가 진행될수록 접촉 면적이 점점 늘어남에 따라 마모 깊이는 줄어든다. 높이 감소량이 점차 줄어들기 때문에 s_a의 감소폭은 줄어들게 됨을 Fig. 8(b)에서 확인 할 수 있다. 세 표면의 RMS 거칠기 초기값은 모두 0.1 μm로 동일하며 해석이 종료되고 난 후에는(m > 0.75) 정규 분포 표면은 s_q = 0.35 μm, Negatively skewed surface는 s_q = 0.43 μm, Positively skewed surface는 s_q = 0.030 μm의 RMS 거칠기 값을 가져 감소량에 차이를 보인다. Positively skewed surface는 s_ku,z가 다른 표면에 비해 높고, s_ku,z가 양수이기 때문에 돌기가 위쪽으로 많이 분포하여 있다. 마모 초기에는 접촉 돌기 수가 적고 돌기 모양으로 인해 접촉 면적이 크지 않아 초기 접촉 압력이 약 2.8 GPa으로 매우 높으며, 이후에도 다른 표면에 비해 높은 값을 갖는다. 따라서 다른 표면에 비해 접촉 압력이 높아 마모 깊이가 깊기 때문에, 이것이 RMS 거칠기로 나타난 것이다. 반대로, Negatively skewed surface는 감소 폭이 가장 작은데, 이 경우에도 표면의 형상으로 설명 가능하다. s_ku,z가 음수이기 때문에 산이 낮고 뭉툭하여 접촉 면적이 크다. 접촉 압력이 Positively skewed surface에 비해 작은 편이며 초기 압력 1 GPa을 시작으로 점점 감소한다. 따라서 돌기의 높이 변화가 작아 RMS 거칠기 감소량이 가장 작다.

Fig. 8. Statistical parameter during wear simulation.(a) Average roughness, (b) RMS roughness, (c) Skewness, (d) Kurtosis.

s_sk는 확률 분포의 비대칭성을 나타내는 지표로, 양수일 경우 확률 분포의 오른쪽에 긴 꼬리가 존재하여 골에 비해 산이 많은 표면을 의미하며 음수일 경우 확률 분포의 왼쪽에 긴 꼬리가 존재하여 골이 많은 표면을 의미한다. Fig. 8(c)에 나타나 있듯이, 세 표면 모두 s_sk,z 사이클이 진행될수록 점점 감소하는 경향을 보인다. 돌기가 평평해지면서 산의 개수보다 골의 개수가 더 많아 지기 때문이다. Positively skewed surface는 s_sk,z = 0.5를 가져 산이 비교적 많은 표면이었으나, 마모가 진행됨에 따라 산이 모두 없어져 s_sk,z가 음수가 되었다. 세 표면 모두 s_sk의 변화 경향은 비슷하지만, Positively skewed surface의 경우 초기에 그래프의 기울기가 급격한 것을 확인할 수 있다. s_sk,z가 양수이기 때문에 뾰족한 산이 많아 접촉 면적이 작기 때문에 접촉 압력이 높다. 따라서 다른 평면에 비해 초기 마모 깊이가 깊어 산이 빠르게 평평해지기 때문에 s_sk,z가 더 크게 감소한다. 이후로는 넓어진 접촉 면적에 의한 영향이 더 커져 다른 평면과 비슷한 경향을 띈다. Negatively skewed surface의 경우 돌기가 낮아 위와 같은 효과가 작게 나타난 것으로 보인다. 결과적으로, Positively skewed surface에 비해 Negatively skewed surface에서 골이 더 많이 남아있기 때문에 마모시뮬레이션 이후에도 s_sk,z가 가장 작으며, 이는 Fig. 9에서도 확인할 수 있다.

Fig. 9. Surface after wear simulation. (a) Gaussian, (b) Negatively skewed, (c) Positively skewed​​​​​​​.

s_ku는 확률 분포의 모양을 나타내는 지표로, 정규분포는 3을 값을 가지며 3보다 크면 분포의 피크가 정규분포보다높아지고, 3보다 작으면 분포의 피크가 낮아진다. 따라서 s_ku,z가 높을수록 돌기의 모양이 뾰족하며, 낮을수록 돌기가 뭉툭해진다. Fig. 8(d)에서 Gaussian surface와 Positively skewed surface는 초기에 감소하는 경향을 보이다가 후에는 8~11까지 증가하였다. 마모가 시작되면서 뾰족한 돌기의 끝이 둥글어 지기 때문에 뾰족함을 나타내는 지표인 s_ku,z가 감소한다. 이후 돌기가 완전히 평평해지면서 s_ku,z가 증가하게 된다. Negatively skewed surface에서는 돌기의 접촉 면적이 넓기 때문에 돌기의 제거에 의한 효과가 적다. 따라서 s_ku,z가 감소하는 구간이 거의 없이 바로 증가하는 양상을 보인다.

4. 결론

본 연구에서는 3차원 거친 표면을 수치적 생성하고 이에 대한 마모 과정을 해석하여, 마모가 진행되는 동안의 거친 표면의 변화를 분석하였다. 마모 해석의 대상이 되는 거친 표면의 생성에는 FIR 필터를 사용한 방법이 채택되었으며 계산 속도를 높이기 위해 합성 곱을 주파수 영역에서 계산하는 FFT를 사용하였다. 비 정규 분포를 갖는 높이를 생성할 때에는 Johnson translation system을 통해 정규 분포 데이터를 임의의 원하는 s_sk,η, s_ku,η를 갖는 비 정규 분포 데이터로 변환하였다. 수치적으로 생성된 표면을 가지고 마모 해석을 수행하였는데, 세 표면 모두RMS 거칠기와 s_sk,z는 감소하였으며 정규 표면 및 Positively skewed surface에서 s_ku,z는 짧은 감소 구간 후 지속적인 증가 추세를 보였다. 이는 기존의 돌기 형상 기반의 마모 모델에서의 거친 표면 파라미터 변화와 같은 경향이다. 표면 파라미터의 변화를 살펴본 결과 형상으로 인한 접촉 압력, 접촉 면적의 차이가 마모에 지배적인 영향을 미치는 것으로 확인되었다.

향후에는 본 연구에서 고려하지 못한 돌기의 접촉 압력 측면에서는 소성 변형을 마모 측면에서는 미끄럼에 따른 접촉 온도 상승과 그에 상응하는 Archard 마모 계수 변화 효과를 고려할 것이며, 더 나아가 돌기가 미끄러지면서 접촉 위치가 변하는 상황까지 고려할 수 있는 모델을 개발하여 거친 표면 접촉 문제를 보다 더 자세히 연구할 것이다.