1. 서론

대부분의 마모시스템들은 정상상태(steady state)에서 요구되는 조건과는 다른 조건으로부터 시작된다. 즉 마모 초기에 접촉압력이 시간에 따라 변화하는 “wearing-in” 상태를 지나 정상상태에 이르게 된다[1,2]. 

이러한 천이 상태의 마모시스템을 분석하기 위해 초기에 접촉역학해석에 기반을 둔 고유치 함수를 전개하는 방법을 개발하였다[3]. 이 방법론은 무한 물체에 대한 접촉압력식을 사용하고 있어서 유한 물체의 마모거동을 파악하는데 있어서 제한적일 수밖에 없다. 이후 많은 연구들이 유한요소방법을 이용하여 천이영역의 마모량을 파악하고자 하였다[4,5]. 특히 이러한 수치해석 프로그램에서는 adaptive meshing기법을 사용하고 있다[6]. 그러나 이러한 방법은 마모가 진행되는 시간에 따라, 설령 세밀하지 않은 시간 간격이라 하더라도, 접촉압력 및 접촉면적, 그리고 구조물의 변위와 응력 등의 모든 정보를 저장하고 활용해야 하므로 수치 해석상의 과도한 계산과 저장공간을 필요로 하게 되어 상당한 부담이 될 수밖에 없다. 이러한 계산을 빠르게 진행하여 수치계산의 용이성을 확보하기 위해 Liu 등[7]은 초기에 제안된 고유치 함수 전개론을 이용하여 간단한 수치해석방법을 제안하였다. 이 연구결과의 장점으로 마모가 진행되기전의 초기접촉압력 정보와 마모계수 만으로 마모천이과정 동안의 접촉압력 및 마모량을 계산할 수 있다. 그러나 이 방법론에 대한 해석적 타당성이 검증되지 않아서 기존 유한 요소 해석 방법에서 활용하는 마모 전과정을 모사하는 해석결과와 비교하고자 한다. 이 비교과정에서 마모고유치 해석 결과를 이용하여 adaptive meshing 기법에서 사용되는 격자 크기에 따른 마모 해석 결과의 민감도 및 격자 개수 등에 대한 지침 등을 제공할 수 있다.

따라서 본 연구에서는 고유치 해석을 이용한 마모 해석 결과와 adaptive meshing 기법에 의한 수치 해석 결과를 비교하여 각 방법의 타당성을 검증하고 마모 해석 모델의 활용을 위한 수치해석적인 기준을 수립하고자 한다.

2. 고유치 해석 및 유한요소 마모 해석

본 연구에서 활용하는 마모해석은 Liu 등[7]에 의해 제시된 고유치 마모해석이며 이는 유한 요소해석 틀을 이용한 해석 방법이다. 이 해석 방법의 결과를 활용하기 위해 다음과 같이 요약한다.

대부분의 마모시스템에서는 마모 초기에 접촉압력이 시간에 따라 변화하는 “wearing-in” 상태로 나타난 후 점차 정상상태(steady state)로 진행된다. 이러한 변화를 반영한 마모시스템의 전체 천이(과도) 해는 다음과 같다. 

\(p(x, y, t)=p_{0}(x, y)+p_{1}(x, y, t)\)      (1)

여기서 p₀(x,y)와 p₁(x,y,t)는 각각 정상상태의 접촉압력과 천이상태의 접촉압력으로서 두번째 항은 시간에 따라 점차 감소하는 경향을 갖는다.

마모시스템의 고유치 특징을 반영하여 다음과 같은 일반 해로 표현할 수 있다.

\(p(x, y, t)=p_{0}(x, y)+\int_{n}^{\infty} C_{n} e^{-\lambda_n t} f_{n}(x, y) \)       (2)

여기서 λ_n은 마모시스템의 무한개의 고유치를 말하며 f_n(x,y)는 상응하는 고유함수이다. 또한 C_n은 초기조건 (예를 들어 t = 0의 접촉압력)으로부터 결정되어지는 상수이다.

문제를 간단히 설명하기 위해 2차원 접촉모델에 대한 평형 방정식은 다음과 같이 구성이 된다.

\(\left\{\begin{array}{l}q \\p\end{array}\right\}=\left[\begin{array}{l}A B^{7} \\B C\end{array}\right]\left\{\begin{array}{l}v \\w\end{array}\right\}+\left\{\begin{array}{l}q^{w} \\p^{w}\end{array}\right\} \)       (3)

여기서 q,p는 각각 접선 방향 및 수직 방향 접촉힘을 나타내며 v,w는 이에 상응하는 변위이다. 또한 q^w,p^w는 접촉에 대한 상대적인 운동이 없을 때 즉 v = w = 0 일 때 발생하는 접촉 힘들이다. 행렬 A,B,C는 전체 강성행렬에 서 얻어지는 부행렬이다. 미끄러짐이 어느 정도 진행 후 초기접촉면 간격 g₀와 수직방향변위 및 마모량의 합에 의해 최종 간격이 0이 된다고 가정하고, 쿨롱의 마찰법칙 q = fp (f는 마찰계수임) 식을 이용하면 다음과 같은 행렬식을 얻을 수 있다.

\(M_{1} p=-M_{2}\left(\delta+g_{0}\right)+p^{w}-B A^{-1} q^{w} \)       (4)

단 여기서

\(M_{1}=I-f B A^{-1}, M_{2}=C-B A^{-1} B^{T}\)       (5)

이제 마모율(rate of wear, \(\dot{\delta}\))이 접촉압력에 비례한다는 Archard 마모 법칙[1]인 \(\dot{\delta}=\alpha L p\)을 적용한다면 다음과 같은 식을 얻는다.

\(M_{1} \dot{p}=-\alpha M_{2} L p\)       (6)

여기서 α는 마모상수이며 L은 노드 i, j에서의 미끄러짐 속도 V_i와 접촉면적 A_j 그리고 크로네커 델타 δ_ij로 이루어진 형태이며 다음과 같이 구성이 된다

\(L_{i j}=\frac{v_{i} \delta_{i j}}{A_{j}}\)     (7)

만일 접촉압력 p가 시간 t에 따라 다음과 같이 구성이 되면

\(p=p_{n} e^{-\lambda_{n} t}\)      (8)

위에서 주어진 식은 다음과 같이 된다.

\(\alpha M_{2} L p_{n}=\lambda_{n} M_{1} p_{n}\)        (9)

이 식은 일반적인 고유치 해석을 위한 식으로서 고유치 값을 알면 다음과 같이 미끄러짐이 일어나는 시간 동안 변할 수 있는 접촉압력을 얻을 수 있다.

\(p(t)=\sum_{n-1}^{N} G_{n} p_{n} e^{-\lambda_{n} t}\)        (10)

상수 G_n을 결정하기 위해 접촉이 시작되는 시간(t = 0)에 마모가 발생하지 않는다는 조건 (δ = 0)을 이용하여 임의의 상수 집합인 G_n은 다음과 같이 구해진다.

\(\sum_{n-1}^{N} G_{n} p_{n}=p_{0}, G=D^{-1} p_{0}\)       (11)

이상의 해석 결과는 식(9)의 고유치 λ_n와 해당 접촉압력 p_n을 결정할 수 있으며 식 (11)을 이용하여 G_n을 결정함으로써 마모가 진행이 되고 있는 시간 동안 접촉 압력을 구할 수 있다. 또한 Archard 마모 법칙을 이용하여 각 시간에 따른 접촉압력으로 마모량을 결정할 수 있다. 

구조해석에서 변형이 심하게 발생하는 곳의 요소를 재구성하여 결과의 정확도를 높이는 기법 중 adaptive meshing 기법이 있다[8]. Adaptive meshing 기법은 우선 접촉압력과 미끄러짐 속도를 도출하여 앞에서 언급한 Archard 마모 법칙을 기반으로 시간당 나타나는 마모량을 계산하고, 마모량을 모델링에 반영한 후 변화된모델에 적합한 형태로 격자 구조를 재배치하는 것을 매해석 증분 시간동안 반복하여 누적 마모량을 도출하는 방법이다. 유한요소해석 프로그램인 ABAQUS에서는 구조해석을 위해Lagrange 해석과 변형된 요소를 재 구성하는 Euler 해석을 이용한다. 이 기법은 두 가지 기본적인 작업을 진행한다. 첫번째 새로운 요소를 만드는 sweeping 과정이며 두번째는 기존 요소로부터 새로운 요소로 변할 때 해석변수를 재구성하는 advection 과정이 있다. 특히 advection 과정에서는 공간과 시간에 대한 변화량의 제곱 차수를 반영하는 Lax-Wendroff 기법을 사용한다. 이러한 기법은 시간에 따라 접촉면적이 변하는 시뮬레이션에서 크게 영향을 끼친다[8]. 

Adaptive meshing 기법을 적용한 유한요소해석과 고유치해석의 비교를 위해 그림 1과 같이 둥근 블록에 대한 격자 모델링을 수행하였다. 이 격자 모델은 둥근 접촉면에서 발생하는 최대접촉압력을 이론식인 식(12)에 근사한 값을 나타내고 있다[9].

Fig. 1. A sliding model with a speed of V.

\(p_{m a x}=\frac{1.80 K_{p}}{\sqrt{d}}=1.073\left(\frac{K_{p}^{2} E^{*}}{R}\right)^{1 / 3}\)       (12)

여기서 \(K_{p}=\frac{P}{\pi \sqrt{2 b}}, \frac{1}{E^{*}}=\frac{1-v_{1}^{2}}{E_{1}}+\frac{1-v_{2}^{2}}{E_{2}}, d=a-b, P R\)은 하중과 둥근 부분의 곡률반경이다. Fig. 1과 같은 모델을 바탕으로 강성행렬과 초기접촉압력 p₀를 도출하여 고유치 해석에 적용하였다.

Fig. 2에서는 고유치 해석과 유한요소해석을 통해 도출한 시간에 따른 접촉압력의 변화 및 누적 마모량을 도시하고 있다. 둥근 블록의 가운데 평면 접촉부에서는 전체적으로 거의 일정한 결과를 보이고 있으며 곡면으로 이루어진 접촉 선단부와 후단부에서는 중간 위치와는 차이가 큰 접촉압력이 나타나고 있다. 이러한 결과를 반영하여 곡면에서 나타나는 최대 접촉 압력과 최대누적마모량을 다른 격자모델과의 비교시 기준으로 삼았다.

Fig. 2. Contact pressure(left) and accumulated wear(right) of eigen analysis and finite element method.

3. 결과 및 고찰

Fig. 1에서 제시된 격자모델의 마모 해석 결과인 Fig. 2의 결과를 기준으로 그림 3과 같이 접촉 곡면 부근 유한요소의 가로 크기를 변경시켜가며 격자모델의 유한요소 민감도해석을 진행하였다. 결과의 타당성 비교를 위해 유한요소 모델에서 고유치 해석에서의 τ가 10, 20, 30, 40일 때의 접촉압력과 누적마모량을 도출하였으며, 최종 τ가 40일 때의 최대접촉압력과 최대누적마모량을 비교하였다. 여기서 τ = αVET/R로 무차원화한 시간 변수이다. 최대접촉압력과 최대누적마모량은 둥근 블록에서의 곡면 부분에서 나타나므로 곡면을 구성하는 유한요소의 가로 크기를 0.2, 0.4, 1.6, 3.2로 증가시켜가며 비교하였다.

Fig. 3. Mesh configurations of optimized(left), coarse (center) and fine(right) mesh, respectively

또한 해석 정확도에는 해석 시간의 수치적 시간 증분의 크기가 영향을 미치는 것으로 알려져 있기 때문에 이를 0.001, 0.004, 0.008, 0.032로 변경하여 해석을 진행하였다. 

위에서 제시된 격자의 가로 크기와 시간 증분을 크기를 동시에 고려하여 다음과 같은 새로운 지표를 정의한다. \(C_{V}=\frac{v d t}{d x}\) 여기에서\( V, dt, dx\)는 각각 슬라이딩 속도, 시간 증분, 접촉면 방향으로의 격자 가로 크기이다. 본연구에서는 시간 증분과 격자 가로 크기의 관계를 분석하기 위해 슬라이딩 속도 V를 300 mm/s에서 500 mm/s의 범위에서 해석을 수행하였다.

C_v에 따른 최대접촉압력과 최대마모량의 상대 오차에 대한 결과는 Fig. 4에 나타나 있다. 여기에서 상대오차란 고유치 해석 결과값을 기준으로 하였을 때, 격자 및 시간 증분의 크기를 변경시켜 만든 유한요소해석 모델에서의 나타나는 결과값이 갖는 차이를 백분율로 나타낸 오차이다.

Fig. 4. Percentage error of contact pressure(upper) and accumulated wear(lower) according to C_V

최대접촉압력과 최적 모델과 가장 근접한 모델은 C_V가 10(속도 500mm/s, 격자 가로 크기 0.4, 시간 증분 0.008)일 때 백분율 오차 6.62%이며 최대마모량의 경우 C_V가5(속도 500 mm/s, 격자 가로 크기 0.2, 시간 증분 0.004)인 모델로 백분율 오차는 0.09%이다. 이때 최적 모델에서 나타나는 최대접촉압력은 729.2054 MPa, 최대마모량은 0.2217 mm이다.

따라서 접촉해석의 정확도를 위해 접촉면을 구성하는 노드의 수를 확보해야 하며, 시간 증분의 크기와 함께 고려된 C_V값 기준 10 이상인 경우 최대마모량 기준 오차 0.15% 이내로 충분한 정확도를 가진 결과를 도출할 수 있다.

4. 결론

본 연구에서는 초기접촉압력 정보와 마모계수 만으로 마모천이과정 동안의 접촉압력 및 마모량을 계산할 수 있는 고유치 해석 결과의 타당성 검증을 위한 유한요소마모해석 기법을 개발하였다. 그 결과, 유한요소 수치해석의 안정성 확보를 위하여 유한요소모델 개발에 지침이 될 수 있는 유한요소의 크기와 수치해석 시간 증분의 크기를 고려한 C_V값을 정의하였다. 다양한 C_V값을 갖는 모델에 대하여 마모해석을 수행하여 C_V값이 10이상인 경우 해석의 정확도를 확보할 수 있음을 보였다.

Acknowledgements 

이 논문은 2020년도 정부(교육부)의 재원으로 한국연구재단의 지원을 받아 수행된 기초연구사업임.(No. 2018R1A2B6008891)