Dislocation in Semi-infinite Half Plane Subject to Adhesive Complete Contact with Square Wedge: Part II - Approximation and Application of Corrective Functions

직각 쐐기와 응착접촉 하는 반무한 평판 내 전위: 제2부 - 보정 함수의 근사 및 응용

Abstract

In Part I, developed was a method to obtain the stress field due to an edge dislocation that locates in an elastic half plane beneath the contact edge of an elastically similar square wedge. Essential result was the corrective functions which incorporate a traction free condition of the free surfaces. In the sequel to Part I, features of the corrective functions, F_kij,(k = x, y;i,j = x,y) are investigated in this Part II at first. It is found that F_xxx(ŷ) = F_xyx(ŷ) where ŷ = y/η and η being the location of an edge dislocation on the y axis. When compared with the corrective functions derived for the case of an edge dislocation at x = ξ, analogy is found when the indices of y and x are exchanged with each other as can be readily expected. The corrective functions are curve fitted by using the scatter data generated using a numerical technique. The algebraic form for the curve fitting is designed as F_kij(ŷ) = $\frac{1}{\hat{y}^{1-{\lambda}}I+yp}$$\sum_{q=0}^{m}{\left}$$\left[A_q\left(\frac{\hat{y}}{1+\hat{y}} \right)^q \right]$ where λ_I=0.5445, the eigenvalue of the adhesive complete contact problem introduced in Part I. To investigate the exponent of F_kij, i.e.(1 - λ_I) and p, Log|F_kij|(ŷ)-Log|(ŷ)| is plotted and investigated. All the coefficients and powers in the algebraic form of the corrective functions are obtained using Mathematica. Method of analyzing a surface perpendicular crack emanated from the complete contact edge is explained as an application of the curve-fitted corrective functions.

Nomenclature

a, b : Length of free surfaces along y and x axes, respectively (y및 x축 상의 자유면 길이)

b_x, b_y : Burgers vector components(버거스 벡터 성분)

c : Crack length (균열 길이)

fⁱ_ij(θ), f^II_ij(θ): Angular function (각도 함수)

p, A_q, q : Exponents and coefficients of the algebraic form of the corrective functions(보정 함수의 대수적 표현 식의 지수 및 계수)

G_kij, (k = x,y;i,j = xy) : Influence functions by Burgers vector component, bk (버거스 벡터 성분bk에 대한 영향 함수)

F_kij, (k = x,y;i,j = xy) : Corrective functions by Burgers vector component, bk(버거스 벡터 성분bk에 대한 보정 함수)

K^contact_I, K^contact_I : Stress intensity factors of complete contact(완전접촉 경계의 응력확대계수)

K^crack_I, K^crack_I : Stress intensity factors of crack(균열 선단의 응력확대계수)

κ : Kolosov constant (Kolosov 상수)

λ_I, λ_II :Eigenvalues (고유값)

σ_rθ, σ_θθ :Stress components (응력 성분)

서론

접촉 경계 부근에서 발생하여 물체 내부로 진전하는 균열은 접촉 응력에 의해 가속될 수 있으며 이때 일반적으로 접촉 응력이 존재하지 않는 경우에 실험적으로 구한 물체의 피로수명이 현저히 감소하게 된다. 이러한 피로현상 중 대표적인 것이 프레팅 피로(fretting fatigue)이다. 여러 연구자들은 프레팅 피로에서의 수명 감소량이 접촉 응력이 없는 경우의 일반 피로수명에 비해 약 5-10%[1], 50%[2], 심지어 90%[3]까지 발생한다고 보고하였다. 따라서 접촉 응력에 의한 균열의 거동을 해석하는 것은 공학 구조물의 안전한 설계를 위해 매우 중요하다.

접촉 응력에 의한 균열 거동을 해석하기 위해서는 우선 접촉 응력을 엄밀하게 해석하는 것이 중요하다. 이때 접촉면의 크기와 접촉 하중의 관계에 따라 완전 접촉 (complete contact)과 불완전 접촉(incomplete contact)으로 구분할 수 있다. 완전접촉은 접촉 하중이 증가하더라도 접촉면 크기의 변화가 없는 경우를 말하며 하중의 증가에 따라 접촉면이 증가하는 경우를 불완전 접촉이라 한다.

불완전 접촉은 접촉 물체의 접촉 경계부에서의 형상이 부드러운 곡선의 형태로 되어 접촉부 경계에서의 접촉 응력이 0이 되는 것을 말한다(대표적인 예는 접촉부의 형상을 2차 함수로 정의할 수 있는 헤르츠 응력[4]이다). 이에 반해 완전 접촉은 접촉하는 두 물체 중 하나의 모서리가 날카로운(ρ → 0인 경우, 여기서 ρ는 모서리 반지름) 경우에 발생한다. 이때 날카로운 모서리로 접근함에 따라 접촉 응력이 무한대가 되는 특징이 있으며 이는 균열 선단에서의 응력 집적(stress intensification)과 유사한 현상이다.

접촉 응력을 해석한 후 접촉부 경계에서 발생하는 균열을 해석하기 위해서는 먼저 균열을 갖는 물체 내부의 응력장을 계산하여야 한다. 이때 불완전 접촉 상태에서의 응력해석 방법은 Flamant 응력 포텐셜을 이용하여 잘 체계화 되어 있다[5]. 그러나 접촉 경계에서 접촉 응력이 무한대가 되는 완전 접촉의 경우는 Williams[6]가 날카로운 모서리를 갖는 노치문제 해석을 위해 개발한 응력 포텐셜과 점근해석(asymptotic analysis)을 이용한 방법이 비교적 근래에 정리되었다[7-9].

물체 내부의 응력장을 구한 후 균열을 해석하기 위해서는 전위밀도함수법(dislocation density function method)[10]을 사용할 수 있는데, 이때 균열위치에서의 응력 성분을 구하는 것이 필요하다. 전위밀도함수법은 균열을 가상 전위(pseudo-dislocation)가 연속적으로 분포되어 있는 것으로 가정하는데, 이에 의한 응력은 전위에 의한 응력포텐셜을 이용하여 계산할 수 있다. 그런데 현재 개발되어 있는 전위에 의한 응력 포텐셜 중 응착 완전접촉 문제에서의 균열해석에 적용가능한 것으로는 Dundurs가 개발한 무한 평편에 하나의 전위가 있는 경우이다[11].

그러나 두 개의 물체가 접촉하는 경우의 응력 해석영역은 두 물체가 접촉하지 않는 자유면이 존재하게 되므로 무한 평판으로 가정할 수 없다. 따라서 자유면에서의 무응력(traction free) 조건을 만족할 수 있도록 하는 방법이 필요하다. 이를 위해 본 논문에서는 보정 함수 (corrective function) 개념을 도입하여 무한 평판에서 전위에 의해 계산되는 응력장을 완전접촉 문제에서 사용할 수 있도록 수정하는 보정 함수를 유도하는 것을 목적으로 하였다.

이에 따라 제1부에서는 보정 함수의 유도방법 및 결과에 초점을 맞추어 기술하였으며 제2부에서는 유도된 보정 함수의 특징을 관찰하고 본 연구의 궁극적인 목적인 균열문제의 해석을 위한 수식화 작업과 이를 위해 필요한 보정 함수의 수치적 근사 방법 및 그 결과에 대해 기술한다. 본 논문에서 균열문제의 최종 해는 수록하지 않았으며 이는 이미 알려진 해가 있어 검증가능한 균열 형상을 찾아 향후 수행하고자 한다. 한편 제2부에서 필요한 수식 중 제1부에서 이미 기술된 것은 반복하여 수록하지 않고 제1부의 수식번호를 언급하여 지정하였다.

2. 보정 함수 분석

2-1. 전위의 위치에 따른 보정 함수의 유사성

제1부에서는 Fig. 1과 같이 모서리가 날카로운 사각 쐐기가 반무한 평판에 접촉하고 이때 접촉 경계로부터 수직 하방으로 성장하는 균열문제를 해석하기 위해 적용할 수 있는 전위에 의한 보정 함수를 유도하였다. 여기서는 Fig. 2와 같이 균열이 접촉면에 존재하는 경우를 고려하여 전위가 x축 상에 있는 경우의 보정 함수를 계산하여 제1부에서의 결과와 비교해 본다.

Fig. 1. Geopetrical description of the present problem: a single edge dislocation beneath the complete contact contact edge of a square wedge.

Fig. 2. Geopetrical description of an edge dislocation and a crack locating at the interface of a complete contact between a square wedge and a half plane of the same material.

Fig. 2와 같이 접촉면에 균열이 있는 문제로서 두 물체의 접합 계면을 생각할 수 있다. 즉 모재에 코팅이 되어 있는 경우 또는 삼차원 프린팅으로 물체를 적층하여 제조하는 경우 계면 균열의 전파는 박리에 의한 파손을 생각할 수 있으며 이때 Fig. 2에 대한 보정 함수를 활용할 수 있을 것이다.

보정 함수에 영향을 미치는 것은 문제의 기하학적 형상과 균열의 형상을 고려한 전위의 위치이다. 따라서 Fig. 2와 같이 전체 해석 영역은 Fig. 1과 동일하지만 균열이 접촉면 상에 존재하고 이의 해석을 위해 x = + ξ에 있는 전위를 고려한 보정 함수는 제1부에서 유도한 보정 함수에 단지 좌표축만이(즉 x와 y에 대해서만) 변경된 유사한 형태로 구할 수 있을 것임을 쉽게 예상할 수 있다. 즉 Fig. 2는 Fig. 1을 시계방향으로 90° 회전한 후 좌우를 바꾸고, x 및 y좌표를 부호를 바꾸면서 상호 교환하면 얻을 수 있다.

이를 확인하기 위해 Fig. 2인 경우의 보정 함수를 제 1부에서와 동일한 방법으로 구하여 보았다. 이때 영향 함수, G_kij(x, y; x_d, y_d), (k, i, j = x, y) (제1부의 식 (12))를 포함하여 제1부에서의 식 (14)-(17)은 정확히 동일하다. 단지 (x,y) = (ξ, 0)에 있는 전위에 의한 자유면에서의 응력을 표현하기 위해 제1부에서의 식 (18)-(21)은 다음과 같이 바뀌게 된다.

i) −b ≤ x ≤ 0 (y = 0) 구간

\(\begin{aligned}\begin{array}{l}\sigma_{y y}^{D H}= \\ \frac{2 G}{\pi(\kappa+1)}\left[G_{x y y}(x, 0 ; \xi, 0) b_{x}(\xi, 0)+G_{y y y}(x, 0 ; \xi, 0) b_{y}(\xi, 0)\right]\end{array}\end{aligned}\)       (1)

\(\begin{aligned}\begin{array}{l}\sigma_{x y}^{D H}=\frac{2 G}{\pi(\kappa+1)}\left[G_{x x y}(x, 0 ; \xi, 0) b_{x}(\xi, 0)+\right. \\ \left.G_{y x y}(x, 0 ; \xi, 0) b_{y}(\xi, 0)\right],\end{array}\end{aligned}\)       (2)

ii) 0 ≤ y ≤ + a(x = 0) 구간

\(\begin{aligned}\begin{array}{l}\sigma_{x x}^{D V}= \\ \frac{2 G}{\pi(\kappa+1)}\left[G_{x x x}(0, y ; \xi, 0) b_{x}(\xi, 0)+G_{y x x}(0, y ; \xi, 0) b_{y}(\xi, 0)\right]\end{array}\end{aligned}\)       (3)

\(\begin{aligned}\begin{array}{l}\sigma_{x y}^{D V}= \\ \frac{2 G}{\pi(\kappa+1)}\left[G_{x x y}(0, y ; \xi, 0) b_{x}(\xi, 0)+G_{y x y}(0, y ; \xi, 0) b_{y}(\xi, 0)\right] .\end{array}\end{aligned}\)       (4)

이제 자유면에서의 무응력의 조건을 만족하기 위한 조건식은 제1부의 식 (22)와 동일한 다음의 연립방정식이 된다.

σ_yy^FH = -σ_yy^DH , σ_xy^FH = -σ_xy^DH,

σ_xx^FV = -σ_xx^DV , σ_xy^FV = -σ_xy^DV. (5)

식 (5)에서 상하첨자의 의미는 제1부에서 이미 설명하였으므로 생략한다.

이후의 수치해석 방법은 제1부의 식 (23)-(31)을 동일하게 이용하면 된다. 이에 따라 식 (5)를 수치계산을 위해 기술한 제1부의 식 (32) 및 그 외 3 개의 식에서 b_x 및 b_y항에 있는 η는 ξ로 적절히 바뀌게 된다. 그러나 중요한 것은 계산에 사용된 수치해석 방법은 완벽히 동일하다는 것이다.

Fig. 2에서 전위는 x축 상에 있으므로 보정 함수는 x축 상의 σ_yy 및 σ_xy에 대해 구한다. 따라서 제1부의 식 (33)과 (34)가 각각 다음의 식 (5)와 (6)과 같이 변하게 되며, 이때 F_xyy(\(\begin{aligned}\hat{x}\end{aligned}\)), F_yyy(\(\begin{aligned}\hat{x}\end{aligned}\)), F_xxy(\(\begin{aligned}\hat{x}\end{aligned}\)) 및 F_yxy(\(\begin{aligned}\hat{x}\end{aligned}\))가 보정 함수이다(여기서 \(\begin{aligned}\hat{x}=x / \xi\end{aligned}\)이다).

\(\begin{aligned}\begin{array}{l} \bullet \; \sigma_{y y}(x, 0 ; \xi, 0) \\ \frac{2 G}{\pi(\kappa+1)}\left[\int_{0}^{a}\left\{G_{x y y}\left(x, 0 ; 0, y_{d}\right) B_{x}\left(y_{d}\right)+\right.\right. \\ \left.G_{y y y}\left(x, 0 ; 0, y_{d}\right) B_{y}\left(y_{d}\right)\right\} d y_{d}+ \\ \int_{-b}^{0}\left\{G_{x y y}\left(x, 0 ; x_{d}, 0\right) B_{x}\left(x_{d}\right)+\right. \\ \left.G_{y y y}\left(x, 0 ; x_{d}, 0\right) B_{y}\left(x_{d}\right)\right\} d x_{d} \\ \left.+G_{x y y}(x, 0 ; \xi, 0) b_{x}(\xi, 0)+G_{y y y}(x, 0 ; \xi, 0) b_{y}(\xi, 0)\right]= \\ \frac{2 G}{\pi(\kappa+1)}\left[\left\{G_{x y y}(x, 0 ; \xi, 0)+F_{x y y}(\hat{x})\right\} b_{x}(\xi, 0)+\right. \\ \left.\left\{G_{y y y}(x, 0 ; \xi, 0) b_{y}(\xi, 0)+F_{y y y}(\hat{x})\right\} b_{y}(\xi, 0)\right]\end{array}\\\end{aligned}\)       (6)

\(\begin{aligned}\begin{array}{l} \bullet\; \sigma_{x y}(x, 0 ; \xi, 0) \\ \frac{2 G}{\pi(\kappa+1)}\left[\int_{0}^{a}\left\{G_{x x y}\left(x, 0 ; 0, y_{d}\right) B_{x}\left(y_{d}\right)+\right.\right. \\ \left.G_{y x y}\left(x, 0 ; 0, y_{d}\right) B_{y}\left(y_{d}\right)\right\} d y_{d}+ \\ \int_{-b}^{0}\left\{G_{x x y}\left(x, 0 ; x_{d}, 0\right) B_{x}\left(x_{d}\right)+\right. \\ \left.G_{y x y}\left(x, 0 ; x_{d}, 0\right) B_{y}\left(x_{d}\right)\right\} d x_{d} \\ \left.+G_{x x y}(x, 0 ; \xi, 0) b_{x}(\xi, 0)+G_{y x y}(x, 0 ; \xi, 0) b_{y}(\xi, 0)\right]= \\ \frac{2 G}{\pi(\kappa+1)}\left[\left\{G_{x x y}(x, 0 ; \xi, 0)+F_{x x y}(\hat{x})\right\} b_{x}(\xi, 0)+\right. \\ \left.\left\{G_{y x y}(x, 0 ; \xi, 0) b_{y}(\xi, 0)+F_{y x y}(\hat{x})\right\} b_{y}(\xi, 0)\right]\end{array}\\\end{aligned}\)       (7)

식 (6) 및 (7)에서 각각의 등호 전후에 있는 영향 함수, G_kij(x, 0; ξ, 0), (k, i, j = x, y)는 소거되므로 결국 보정 함수는 자유면에 연속적으로 분포한 전위에 의한 응력 항이 된다. 제1부의 수치적 방법을 이용하여 전위가 x = + ξ에 있는 경우의 보정 함수를 유도하면 다음과 같다.

\(\begin{aligned}\begin{array}{l}F_{x y y}(x, \xi) \equiv F_{x y y}(\hat{x}) \\ =\sum_{i=1}^{n}\left[\left(-\frac{a}{2}\right)\left\{G_{x y y}\left(x, 0 ; 0,\left(-\frac{a\left(u_{y i}-1\right)}{2}\right)\right)\right.\right. \\ \quad \times W\left(u_{x i}\right) \phi_{x}\left(u_{y i}\right)+G_{y y y}\left(x, 0 ; 0,\left(-\frac{a\left(u_{y i}-1\right)}{2}\right)\right) \\ \left.\times W\left(u_{y i}\right) \phi_{y}\left(u_{y i}\right)\right\} \\ \quad+\left(\frac{b}{2}\right)\left\{G_{x y y}\left(x, 0 ;\left(\frac{b\left(u_{x i}-1\right)}{2}\right), 0\right)\right. \\ \quad \times W\left(u_{x i}\right) \phi_{x}\left(u_{x i}\right) \\ \left.\left.\quad+G_{y y y}\left(x, 0 ;\left(\frac{b\left(u_{x i}-1\right)}{2}\right), 0\right) \times W\left(u_{y i}\right) \phi_{y}\left(u_{x i}\right)\right\}\right]\end{array}\end{aligned}\)       (8)

\(\begin{aligned}\begin{array}{l}F_{y y y}(x, \xi) \equiv F_{y y y}(\hat{x}) \\ =\sum_{i=1}^{n}\left[\left(-\frac{a}{2}\right)\left\{G_{x y y}\left(x, 0 ; 0,\left(-\frac{a\left(u_{y i}-1\right)}{2}\right)\right)\right.\right. \\ \quad \times W\left(u_{x i}\right) \phi_{x}\left(u_{y i}\right)+G_{y y y}\left(x, 0 ; 0,\left(-\frac{a\left(u_{y i}-1\right)}{2}\right)\right) \\ \left.\times W\left(u_{y i}\right) \phi_{y}\left(u_{y i}\right)\right\} \\ \quad+\left(\frac{b}{2}\right)\left\{G_{x y y}\left(x, 0 ;\left(\frac{b\left(u_{x i}-1\right)}{2}\right), 0\right)\right. \\ \quad \times W\left(u_{x i}\right) \phi_{x}\left(u_{x i}\right) \\ \left.\left.\quad+G_{y y y}\left(x, 0 ;\left(\frac{b\left(u_{x i}-1\right)}{2}\right), 0\right) \times W\left(u_{y i}\right) \phi_{y}\left(u_{x i}\right)\right\}\right]\end{array}\end{aligned}\)       (9)

\(\begin{aligned}\begin{array}{l}F_{x x y}(x, \xi) \equiv F_{x x y}(\hat{x}) \\ =\sum_{i=1}^{n}\left[\left(-\frac{a}{2}\right)\left\{G_{x x y}\left(x, 0 ; 0,\left(-\frac{a\left(u_{y i}-1\right)}{2}\right)\right)\right.\right. \\ \quad \times W\left(u_{x i}\right) \phi_{x}\left(u_{y i}\right)+G_{y x y}\left(x, 0 ; 0,\left(-\frac{a\left(u_{y i}-1\right)}{2}\right)\right) \\ \left.\times W\left(u_{y i}\right) \phi_{y}\left(u_{y i}\right)\right\} \\ \quad+\left(\frac{b}{2}\right)\left\{G_{x x y}\left(x, 0 ;\left(\frac{b\left(u_{x i}-1\right)}{2}\right), 0\right)\right. \\ \quad \times W\left(u_{x i}\right) \phi_{x}\left(u_{x i}\right) \\ \left.\left.\quad+G_{y x y}\left(x, 0 ;\left(\frac{b\left(u_{x i}-1\right)}{2}\right), 0\right) \times W\left(u_{y i}\right) \phi_{y}\left(u_{x i}\right)\right\}\right]\end{array}\end{aligned}\)       (10)

\(\begin{aligned}\begin{array}{l}F_{y x y}(x, \xi) \equiv F_{y x y}(\hat{x}) \\ =\sum_{i=1}^{n}\left[\left(-\frac{a}{2}\right)\left\{G_{x x y}\left(x, 0 ; 0,\left(-\frac{a\left(u_{y i}-1\right)}{2}\right)\right)\right.\right. \\ \quad \times W\left(u_{x i}\right) \phi_{x}\left(u_{y i}\right)+G_{y x y}\left(x, 0 ; 0,\left(-\frac{a\left(u_{y i}-1\right)}{2}\right)\right) \\ \left.\quad \times W\left(u_{y i}\right) \phi_{y}\left(u_{y i}\right)\right\} \\ \quad+\left(\frac{b}{2}\right)\left\{G_{x x y}\left(x, 0 ;\left(\frac{b\left(u_{x i}-1\right)}{2}\right), 0\right)\right. \\ \quad \times W\left(u_{x i}\right) \phi_{x}\left(u_{x i}\right)\end{array}\end{aligned}\)       (11)

식 (8)-(11)의 결과를 Fig. 3에 보여 준다. 제1부에서와 같이 계산은 상용 프로그램 Mathematica [12]를 사용하였고 이때, 제1부에서 수렴성이 확인된 n = 500, N(= a/η, b/η) = 80으로 동일하게 하였다.

Fig. 3. Corrective functions (dislocation at x = ξ).

한편 용이한 비교를 위해 Fig. 4에 제1부의 Fig. 5를 여기에 다시 도시한다.

Fig. 4. Corrective functions (dislocation at y = -η). Duplication of Fig. 4 of Part I.

Fig. 5. Log plot of the corrective functions (dislocation at y = -η).

Fig. 3과 Fig. 4를 비교할 때 서로 동일한 형태의 보정 함수는 Table 1과 같이 정리할 수 있다(각각의 변수인 \(\begin{aligned}\hat{x}=(x / \xi)\end{aligned}\)와 ŷ = y / η 표기를 생략함). 앞서 예상한 대로 Fig. 3과 Fig. 4의 보정 함수는 x와 y를 서로 교환하면 동일하다는 것을 확인할 수 있다.

Table 1. Analogy of the corrective functions

2-2. 보정 함수의 근사

식 (8)-(11) (또는 제1부의 식 (35)-(38))과 같이 구한 보정 함수는 식 (5)의 연립방정식의 해인 𝜙_x(u_x), 𝜙_x(u_y), 𝜙_y(u_x), 및 𝜙_y(u_y)로부터 얻었으므로 대수 함수(algebraic function)의 형태를 갖지 않는다. 따라서 향후 전위의 배열로 모사하게 될 균열 문제를 해석하기 위해서는 Fig. 3 또는 Fig. 4의 보정 함수를 대수 함수의 꼴로 바꾸어야한다. 이를 위해 본 논문에서는 Mathematica[12]에서 제공하는 기능으로서 “ListPlot” 및 “Table”을 사용하여 Fig. 3(또는 제1부의 Fig. 4)의 그래프 데이터를 취득한 후 “FindFit”을 사용하여 curve fitting함으로써 대수 함수 형태의 근사식을 생산하였다.

이때 근사식은 미리 요구하는 형태의 함수를 도출하여 이를 적용하는 것이 논리적이다. 이러한(요구하는 형태의 함수 형태를 도출하는) 작업은 그 함수가 현재의 문제에 맞는 논리적 타당성을 가져야 하는 동시에 향후 균열 문제 해석에서 등장하게 되는 다른 함수를 고려하여 큰 어려움 없이 계산할 수 있을 것이 요구된다. 이를 위해 다음과 같은 점을 고려하였다.

i) 전위의 위치

향후 해석하게 될 균열이 길이가 짧은 경우와 긴 경우를 구분할 필요가 있다. 그 이유는 균열을 진전하게 하는 응력이 접촉 경계에서 가까운 경우와 먼 경우에 있어 많은 차이가 발생하기 때문이다. 즉 본 연구의 주제인 완전접촉 문제에서 접촉경계의 응력장은 제1부의 식 (5)로 표현되고, 이때 접촉 경계로 접근할수록(r → 0) 아래 식(12)와 같이 특이 항(singular terms)만이 유효하게 되고 그 이외의 항은 무시될 수 있기 때문이다.

σ_ij = K_I^contactγ^λ_I-1f_ij^I(θ) + K_II^contactγ^λ_II-1f_ij^II(θ).       (12)

한편, 식 (6)과 (7)에서 알 수 있듯이 영향 함수와 보정 함수(각각 G_Kij 및 F_Kij, k, i, j = x, y)는 응력의 성격을 갖게 되므로 접촉 경계로 접근할수록 접촉 응력장의 특이성 지수(r^{λI - 1}))¹를 가져야 할 필요가 있다. 현재의 문제에서 λ_I = 0.5445 (제1부 식 (8) 참조)이다. 또한 전위가 접촉 경계로부터 먼 경우에는 제1부의 식 (5)와 같이 특이(singular) 항에 비특이(bounded) 항이 더해진 응력장이 형성된다. 따라서 접촉 경계로부터의 거리(r)에 임의의 지수가 사용될 필요가 있다. 이러한 조건으로부터 curve fitting을 위한 보정 함수의 대수 함수꼴이 다음의 표현이 포함되도록 한다.

\(\begin{aligned}\frac{1}{\hat{y}^{1-\lambda_{I+\hat{y}^{p}}}},(k=x, y ; i, j=x, y)\end{aligned}\)       (13)

여기서 ŷ = y/η이다.

ii) 다항식의 사용

일반적으로 curve fitting에 의해 대수 함수를 생산할 때 동일한 변수의 멱함수(power function)를 이용하면 문제의 해석 영역(여기서는 x-y 평면 또는 r-θ 평면) 전체를 포함할 수 있으므로 편리하다. 이에 따라 보정 함수는 식 (13)에 멱함수가 곱해진 형태가 되도록 하였다. 이때 접촉 경계로 접근할 때(r → 0) 응력의 지수가 (λ_I –1)이 되어야 하므로 식 (13)과 곱해질 멱급수가 r → 0일때 상수가 될 필요가 있다. 이를 고려한 멱급수의 형태를 다음과 같이 제안할 수 있다.

\(\begin{aligned}\sum_{q=0}^{m}\left[A_{q}\left(\frac{\hat{y}}{1+\hat{y}}\right)^{q}\right]\end{aligned}\)       (14)

여기서A_q는 상수로서 curve fitting에 의해 구해질 것이며 \(\begin{aligned}\left(\frac{\hat{y}}{1+\hat{y}}\right)^{q}\end{aligned}\) 의 계수가 된다. 이로부터 보정 함수의 최종적인 형태를 다음과 같이 결정하였다.

\(\begin{aligned}F_{k i j}=\frac{1}{\hat{y}^{1-\lambda_{I+\hat{y}}}} \sum_{q=0}^{m}\left[A_{q}\left(\frac{\hat{y}}{1+\hat{y}}\right)^{q}\right],(k=x, y ; i, j=x, y) \end{aligned}\).       (15)

식 (15)를 접촉 경계에 접근하는 경우와 경계로부터 멀리 떨어진 경우에 대해 표현하면 다음과 같이 된다.

접촉 경계에 접근하는 경우(ŷ « 1 )

\(\begin{aligned}F_{k i j}=\frac{A_{0}}{\hat{y}^{1-\lambda_{I}}}, \because(\hat{y} /(1+\hat{y})) \rightarrow 0\end{aligned}\)       (16)

접촉 경계에서 먼 경우(ŷ » 1 )

\(\begin{aligned}F_{k i j}=\frac{1}{\hat{y}^{p}} \sum_{q=0}^{m}\left[A_{q}\right], \because(\hat{y} /(1+\hat{y})) \rightarrow 1\end{aligned}\)       (17)

식 (14)-(17)에서 m이 증가함에 따라 더욱 근사한 함수를 얻을 수 있을 것이다. 본 논문에서는 m = 3으로 하여 근사한 정도를 검토하였다. 한편, Fig. 3과 같이 전위가 x = ξ에 있고 접촉면 상에 균열이 있는 문제의 경우에는 위의 식 (13)-(17)에서의 ŷ 을 \(\begin{aligned}\hat{x} \quad(=x / \xi)\end{aligned}\)로 치환하면 된다.

식 (13)에서 ŷ의 지수는 Fig. 4에서 구한 보정 함수를 Log-Log로 변환한 후 그 기울기를 찾으면 된다. 이때 음(−)의 영역에 있는 F_kij를 표현하기 위해 절대값을 취하는 점에 주의한다. Fig. 5에 그 결과를 보여 준다.

Fig. 4에서 예상할 수 있듯이, Fig. 5에서 −ŷ 가 증가하는 동안(y축을 따라 아래로 내려갈 때)보정 함수가 음(−)에서 양(+)으로 바뀌는 F_xxy와 F_yxx는 중간에0으로 감소했다 증가하는 거동이 존재하고 계속 음(−)의 영역에 있는 F_xxx와 F_yxy는 그러한 거동을 보이지 않는 것을 확인할 수 있다.

Fig. 5의 curve fitting에서 Log(−ŷ) → 0인 부분의 기울기를 (λ_I − 1) (= −0.4555)로 하고 Log(−ŷ) → ∞(Fig. 5에서 Log(−ŷ) ≥ ~100)에서의 기울기를 식 (13)에서의 지수 p로 추출한다. 예로서, Log|F_yxxŷ|에 대해 각각의 기울기를 Fig. 6에 도시한다.

Fig. 6. Slopes of Log|F_yxx(ŷ)|.

실제 Mathematica를 이용하여 식 (15)의 형태로 각각의 보정 함수를 근사할 때에는 λ_I = 0.5445와 m = 3을 지정하면 p, A_q, (q = 0,1,2,3)가 계산된다. Fig. 7은 Mathematica가 계산한 하나의 예로서, F_yxx(ŷ)의 데이터(점으로 표시됨)와 이를 근사한 대수 함수(선으로 표시됨)를 보여준다.

Fig. 7. Generated data (scatter points) and its curve fit (line) of F_yxx(ŷ) carried out by Mathematica.

이와 같은 방법으로 각각의 보정 함수의 대수적 표현식의 각 변수 값을 Table 2에 정리하였다. 이때 식 (15) 에서 λ_I = 0.5445으로 고정하였으므로 별도로 표기하지 않았다.

Table 2. Power and coefficient values of the curve fitted corrective functions

3. 균열문제 해석 방법

균열 문제를 해석하기 위한 방정식은 제1부의 식 (3)에 기술한 것과 같이 균열면을 열리게 하는 수직 응력(σ_θθ)과 균열면을 따라 미끄럼이 발생토록 하는 전단 응력(σ_rθ)이 균열면 상에서 0이 되어야 하는 무응력 조건 식으로 구성한다(제1부 식 (3) 참조).

본 논문에서는 Fig. 1과 같이 접촉 경계에서 연직 하방으로 진전하는 길이, c의 균열을 고려하고 있다. 따라서, 데카르트(Cartesian) 좌표로 나타낼 때 균열면에 대한 수직 응력은 σ_xx이고 전단 응력은 σ_xy이다. 한편 반무한 평판의 크기에 비해 균열의 길이가 짧다고 생각할 때, 균열면에서의 수직 및 전단 응력은 다음과 같이 나타낼 수 있게 된다.

\(\begin{aligned}\begin{array}{l}\sigma_{\theta \theta}^{C}=\sigma_{x x}^{C}= \\ K_{I}^{\text {contact }} y^{\lambda_{I}-1} f_{\theta \theta}^{I}\left(-\frac{\pi}{2}\right)+K_{I I}^{\text {contact }} y^{\lambda_{I I}-1} f_{\theta \theta}^{I I}\left(-\frac{\pi}{2}\right)\end{array}\end{aligned}\)       (18)

\(\begin{aligned}\begin{array}{l}\sigma_{r \theta}^{C}=\sigma_{x y}^{C}= \\ K_{I}^{\text {contact }} y^{\lambda_{I}-1} f_{r \theta}^{I}\left(-\frac{\pi}{2}\right)+K_{I I}^{\text {contact }} y^{\lambda_{I I}-1} f_{r \theta}^{I I}\left(-\frac{\pi}{2}\right)\end{array}\end{aligned}\)        (19)

여기서 응력 성분의 상첨자 C는 제1부에서와 같이 완전 접촉에 의한 것을 의미한다.

한편 균열 길이 c의 구간에 전위가 연속적으로 분포한다고 생각하여, 균열에 의한 수직 및 전단 응력을 다음과 같이 쓸 수 있다(식 (6) 및 (7)을 참조하면 쉽게 유도할 수 있다).

\(\begin{aligned}\begin{array}{l}\sigma_{\theta \theta}^{D}=\sigma_{x x}^{D}=\frac{2 G}{\pi(\kappa+1)}\left[\int_{0}^{c} B_{x}(\eta)\left\{G_{x x x}\left(0, y ; 0, y_{d}\right)+\right.\right. \\ \left.F_{x x x}(y, \eta)\right\} d \eta+\int_{0}^{c} B_{y}(\eta)\left\{G_{y x x}\left(0, y ; 0, y_{d}\right)+\right. \\ \left.\left.F_{y x x}(y, \eta)\right\} d \eta\right]\end{array}\end{aligned}\)       (20)

\(\begin{aligned}\begin{array}{l}\sigma_{r \theta}^{D}=\sigma_{x y}^{D}=\frac{2 G}{\pi(\kappa+1)}\left[\int_{0}^{c} B_{x}(\eta)\left\{G_{x x y}\left(0, y ; 0, y_{d}\right)+\right.\right. \\ \left.F_{x x y}(y, \eta)\right\} d \eta+\int_{0}^{c} B_{y}(\eta)\left\{G_{y x y}\left(0, y ; 0, y_{d}\right)+\right. \\ \left.\left.F_{y x y}(y, \eta)\right\} d \eta\right]\end{array}\end{aligned}\)       (21)

여기서 B_x(η), B_y(η)는 균열면에 분포시킨 전위의 Burgers vector 성분, b_x, b_y 각각에 관련한 전위밀도함수이며, 응력 성분의 상첨자 D는 전위에 의한 것임을 나타낸다. 한편 식 (20) 및 (21)에서 영향 함수, G_kij(0, y; 0, y_d), (k = x, y; i, j = x, y)는 제1부의 식 (12)를 이용하여 다음과 같이 계산된다.

\(\begin{aligned}G_{x x x}\left(0, y ; 0, y_{d}\right)=G_{y x y}\left(0, y ; 0, y_{d}\right)=-\frac{1}{y-\eta}\end{aligned}\)       (22)

G_yxx(0, y;0, y_d) = G_xxy(0, y;0, y_d) = 0       (23)

따라서 균열 해석을 위한 최종 식은 식 (18)과 (20), 그리고 식 (19)와 (21)을 각각 더한 것이 0이 되어야 한다는 무응력 조건으로부터 다음 식 (24) 및 (25)와 같이 된다.

\(\begin{aligned}\begin{array}{l}K_{I}^{\text {contact }} r^{\lambda_{I}-1} f_{\theta \theta}^{I}\left(-\frac{\pi}{2}\right)+K_{I I}^{\text {contact }} r^{\lambda_{I I}-1} f_{\theta \theta}^{I I}\left(-\frac{\pi}{2}\right)+ \\ \frac{2 G}{\pi(\kappa+1)}\left[\int_{0}^{c} B_{x}(\eta)\left\{-\frac{1}{y-\eta}+F_{x x x}(y, \eta)\right\} d \eta+\right. \\ \left.\int_{0}^{c} B_{y}(\eta) F_{y x x}(y, \eta) d \eta\right]=0\end{array}\end{aligned}\)       (24)

\(\begin{aligned}\begin{array}{l}K_{I}^{\text {contact }} r^{\lambda_{I}-1} f_{r \theta}^{I}\left(-\frac{\pi}{2}\right)+K_{I I}^{\text {contact }} r^{\lambda_{I I}-1} f_{r \theta}^{I I}\left(-\frac{\pi}{2}\right)+ \\ \frac{2 G}{\pi(\kappa+1)}\left[\int_{0}^{c} B_{x}(\eta) F_{x x y}(y, \eta) d \eta+\int_{0}^{c} B_{y}(\eta)\left\{-\frac{1}{y-\eta}+\right.\right. \\ \left.\left.F_{y x y}(y, \eta)\right\} d \eta\right]=0\end{array}\end{aligned}\)       (25)

여기서 네 개의 보정 함수, F_xxx(y, η), F_yxx(y, η), F_xxy(y,η), F_yxy(y, η)는 Table 2의 값이 적용된 식 (15)가 된다(괄호 (y, η)는 (ŷ)와 같은 의미이다). 예를 들어, F_yxx(y,η)는 다음과 같다(m = 3인 경우).

\(\begin{aligned}\begin{array}{l}F_{y x x}(\hat{y})=\frac{1}{\hat{y}^{0.4555}+\hat{y}^{1.187}}\left[-0.902+3.832\left(\frac{\hat{y}}{1+\hat{y}}\right)-\right. \\ \left.3.186\left(\frac{\hat{y}}{1+\hat{y}}\right)^{2}+0.265\left(\frac{\hat{y}}{1+\hat{y}}\right)^{3}\right]\end{array}\end{aligned}\)       (26)

식 (24)와 (25)를 풀어 균열의 응력확대계수를 구하기 위해서는 제1부에서 보정 함수를 구하기 위해 사용하였던 수치해석 방법을 적용할 수 있다. 즉, Gauss-Jacobi 적분 공식을 이용하여[13], 각각 η와 y를 c에 대해 정규화 한 후 균열 구간 내 적분점과 collocation point를 생성하고(제1부의 식 (23), (24), (28), (29) 참조), 또 전위밀도함수 B_x(η) 및 B_y(η)를 fundamental function, w(u) (제1부의 식 (27) 및 (30), (31) 참조)와 구하고자 하는 미지 함수, 𝜙_x(u_i)와 𝜙_y(u_i)의 곱으로 한 후(제1부의 식(26) 참조) 연립방정식을 구성하여 𝜙_x(u_i)와 𝜙_y(u_i)를 구한다. 이때 영향 함수 및 보정 함수의 변수(η 및 y)도 같은 방법으로 치환한다.

최종적으로 균열의 응력확대계수는 제1부의 식 (1) 또는 (2)의 정의에 의해 다음의 공식으로 구할 수 있다[10].

\(\begin{aligned}K_{I}^{\text {crack }}=\frac{2 \sqrt{2} G}{\kappa+1} \sqrt{\pi c} \phi_{x}(1), K_{I I}^{\text {crack }}=\frac{2 \sqrt{2} G}{\kappa+1} \sqrt{\pi c} \phi_{y}(1)\end{aligned}\)       (27)

여기서 𝜙_k(1), (k = x, y)는 다음의 보간법을 사용하여 구할 수 있다[14].

\(\begin{aligned}\begin{array}{l}\phi_{k}(1)= \\ \frac{2}{2 n+1} \sum_{i=1}^{n}\left[\cot \left\{\left(\frac{2 i-1}{2 n+1}\right) \frac{\pi}{2}\right\} \sin \left\{\left(\frac{2 i-1}{2 n+1}\right) n \pi\right\} \phi_{k}\left(u_{i}\right)\right] \\ (k=x, y)\end{array}\end{aligned}\)       (28)

한편 식 (24) 및 (25)로부터 그 해로 구하게 되는 𝜙_x(u_i)와 𝜙_y(u_i)는 완전접촉 경계에서의 응력확대계수인 K_I^comntact과 K_II^comntact을 포함하게 된다. 따라서 본 논문에서의 균열 응력확대계수는 다음과 같은 형태가 될 것임을 예측할 수 있다.

\(\begin{aligned}\begin{array}{l}K_{I}^{\text {crack }} / \sqrt{\pi c}=N_{1} K_{I}^{\text {contact }} c^{\lambda_{I}-1}+N_{2} K_{I I}^{\text {contact }} c^{\lambda_{I I}-1} \\ K_{I I}^{\text {crack }} / \sqrt{\pi c}=S_{1} K_{I}^{\text {contact }} c^{\lambda_{I}-1}+S_{2} K_{I I}^{\text {contact }} c^{\lambda_{I I}-1}\end{array}\end{aligned}\)       (29)

여기서, N₁, N₂, S₁, S₂는 상수이고, 제1부의 식 (8)에 따라 λ_I = 0.5445, λ_II = 0.9085이다.

본 논문은 응착 완전접촉 상태에서의 균열 문제를 해석할 수 있는 전위에 대해 기술하는 것이 목적이므로 실제의 균열 문제를 해석하여 식 (29)와 같은 균열의 응력 확대계수를 생산하지는 않았다. 이에 대해서는 향후 엄밀해가 있는 Fig. 1 또는 2와 같은 형상의 균열 문제를 찾아 본 논문에서의 방법을 사용하여 검증할 계획이다.

4. 결론

본 논문의 제1부에서 유도한 보정 함수에 대해 제2부에서는 각 보정 함수의 특성과 균열문제에 적용하기 위한 대수적 함수를 구하고 균열의 응력확대계수를 구할수 있는 방법에 대해 기술하였다.

제1부에서 구한, 완전접촉 경계의 연직 하방으로 성장하는 균열의 해석에 필요한 보정 함수는 좌표축(x, y)을 교환하는 것만으로 접촉경계에서 접촉면을 따라 성장하는 균열 해석의 보정 함수가 된다는 것을 확인하였다. 이것은 접촉경계에서 시작된 균열이 어떤 방향으로 진전하든 그에 대한 보정 함수는 그 균열에 직각 방향으로 진전하는 균열의 보정 함수와 동일한 대수적 함수의 형태를 갖게 된다는 것을 의미하므로 해석에 매우 유용하게 사용할 수 있다.

한편, Mathematica에 의한 수치적 curve fitting을 수행하여 보정 함수의 대수적 형태를 결정하고 균열문제의 해석방법에 적용하였다. 따라서 본 제2부 논문은 제1부와 함께 응착 완전접촉 상태에 있는 기계구조물의 트라이볼로지적 균열 손상문제를 해석하는 데에 기본적인 도구를 제공한다.

향후의 연구로서, 본 논문에서 개발한 방법을 이용하여 실제 완전접촉 상태에서 발생하는 균열의 응력확대 계수를 구하고자 한다. 이를 위해 본 논문에서 해석하고 있는 기하학적 형상(Fig. 1 및 2)과 동일한 균열 문제의 엄밀해를 찾아 본 논문에서의 방법론을 검증할 필요가 있으며 이때 제2부의 제3절을 함께 상술할 것이다.