1. 서 론

균열선단 부근의 응력 특이성을 지배하고 공학 재료에서 균열의 개시 및 전파를 결정하는 응력확대계수(Stress Intensity Factor, SIF)의 정확한 평가는 파괴역학에서 기본적인 요구사항이다(Anderson, 2017). 응력확대계수를 계산하기 위해 경계요소법(Aliabadi, 2002), 일반/확장 유한요소법(GFEM/XFEM)(Deng et al., 2022; Duarte et al., 2001; Gupta et al., 2015; Kim and Paulino, 2003; Marco et al., 2020), 그리고 무요소법(Meshfree Method)(Belytschko et al., 1994; Chen et al., 2001; Kim et al., 2014) 등 다양한 해석기법이 사용되어 왔다. 다양한 접근법 중에서 전통적인 유한요소법(Finite Element Method, FEM)이 복잡한 기하 형상과 하중 조건을 처리하는 편리함으로 인해 실무에서 널리 사용되고 있다. 그러나 FEM을 이용한 응력확대계수 평가의 정확도는 요소망의 품질과 이로부터 계산되는 응력장의 정확도에 크게 의존한다(Daimon and Okada, 2014; Nagai et al., 2013; Tabaza et al., 2021).

FEM에서는 균열 선단 부근의 응력 특이성을 정확하게 포착하기 위해 일반적으로 구조화된 요소망이 요구되지만(Bouchard et al., 2000; Ingraffea and Manu, 1980; Kim and Paulino, 2004; Lee et al., 2022), 복잡한 균열 형상에서는 이러한 요소망의 생성이 어렵다. 이에 대한 대안으로 비정형 요소망은 높은 유연성을 제공하지만, 균열 선단 부근의 응력 평가 정확도가 저하될 수 있다(Krueger, 2004; Nagai et al., 2013; Paluszny and Zimmerman, 2011; Tabaza et al., 2021). 이러한 요소망에서는 변위장의 수치 미분 과정에서 계산된 응력에 오차가 발생할 수 있으며, 이는 응력확대계수의 정확도에 직접적인 영향을 미친다.

비정형 요소망에서 응력확대계수를 정확하게 계산하는 접근법 중 하나로 가상 격자 응력 복원(Virtual Grid Stress Recovery, VGSR) 기법(Choi et al., 2022; Movahedi et al., 2026b)이 있으며, 이는 기존 유한요소망과 독립적으로 응력장 평가를 향상시키기 위한 효과적인 프레임워크를 제공한다. VGSR 방법은 왜곡되거나 불규칙한 요소와 관련된 오차를 크게 감소시키고, 보다 매끄럽고 신뢰성 있는 응력 분포와 응력확대계수를 제공하는 것으로 확인되었다(Choi et al., 2022; Movahedi et al., 2026a, 2026b). 이러한 장점에도 불구하고, VGSR에서 복원된 응력장의 정확도는 여전히 유한요소 해석으로부터 얻어진 변위장의 품질에 의존한다. 특히 균열 선단과 같이 응력 구배가 큰 영역에서 초기 요소망이 충분히 조밀하지 않은 경우, 보간된 변위장은 국부적인 거동을 정확하게 반영하지 못할 수 있으며, 이는 전체적인 정확도 저하로 이어진다.

변위장의 정확도를 효율적으로 향상시키기 위해 중요한 영역에서 유한요소 해의 해상도를 향상시키기 위한 적응적 요소 세분화(Adaptive Mesh Refinement, AMR)를 적용할 수 있다(Park et al., 2012; Zienkiewicz and Zhu, 1987). 적응적 요소 재구성(remeshing) 기법은 오차 지표를 기반으로 요소망을 국부적으로 세분화하여 계산 비용을 절감할 수 있다. 본 연구에서는 von Mises 응력을 세분화 기준으로 채택하였다(Katragadda and Grosset, 1996; Zienkiewicz and Zhu, 1987). 이 방법은 Abaqus(Simulia, 2014)와 같은 상용 유한요소해석 소프트웨어에서 활용 가능하며, 균열 선단 및 균열 전방과 같이 높은 응력 구배가 발생하는 영역에 요소를 집중시키면서 반복적으로 요소망의 해상도를 개선할 수 있다.

본 연구에서는 VGSR 기법과 적응적 요소 세분화과정을 결합하여 혼합모드 응력확대계수 평가를 단순하고 정확하게 수행하고자 한다. 이 결합된 접근법은 혼합모드 하중 조건을 받는 경사진 penny-shaped crack의 해석에 적용된다. 응력확대계수는 VGSR에서 얻은 응력장을 기반으로 I-적분(Interaction integral)을 사용하여 계산하며, 계산된 혼합모드 응력확대계수의 정확도와 수렴성에 미치는 영향을 분석한다.

2. 해석 기법

2.1 가상격자 응력복원(VGSR) 기법

본 연구에서는 파괴 매개변수 평가에 필요한 보다 정확하고 매끄러운 응력장을 얻기 위해 VGSR 기법을 도입하였다. 특히 사면체 요소로 구성된 비정형 3차원 요소망을 기준으로 수치 미분을 할 경우 발생하는 매끄럽지 않은 응력장을 매끄럽게 만들어 정확한 파괴 매개변수 평가가 가능하다. 이 기법은 기존 유한요소망의 형태에 종속되지 않고 적용될 수 있어, 임의의 균질한 영역에서 수치적분을 할 수 있다는 장점이 있다.

VGSR 방법을 균열 선단 주변에 구조화된 가상 격자(virtual grid, VG)를 생성하고 가상격자 내에서 변위장과 응력장을 평가한다. 가상 격자의 크기(lVG)는 유한요소 크기(lFE)를 기준으로 설정할 수 있다. 가상 격자는 8절점 육면체(HEX8) 요소로 구성되며, 먼저 일정한 크기의 가상 격자가 생성된 후, 실제 균열 선단의 각 절점을 기준으로 정의된 국부 좌표계(local coordinate system)를 이용하여 유한요소 모델에 중첩된다. 가상 격자의 방향은 국부적인 균열 형상에 따라 정의되며, 격자가 균열 선단 및 그 법선 방향과 정렬되도록 설정된다(Fig. 1 참조). 또한, 곡선형 균열 선단의 형상을 보다 정확하게 반영하기 위해 가상 격자의 절점 위치를 균열선단과 일치하도록 국부적으로 조정할 수 있도록 하였다. Fig. 1(a)는 원래의 비정형 사면체 유한요소망을 나타내고, Fig. 1(b)는 균열 선단 주변에 구성된 구조화된 가상 격자를 보여준다. Fig. 1(c)는 가상 격자를 유한요소망 위에 중첩하여 표시한 것으로 가상 격자 절점에서 변위장을 복원하는 과정을 나타낸다.

Fig. 1.

Concept of the VGSR method: (a) unstructured tetrahedral finite element mesh, (b) structured virtual grid composed of hexahedral elements around the crack front, and (c) superposition of the virtual grid and finite element mesh illustrating the interpolation of nodal displacements from the finite element to the virtual grid nodes

가상 격자의 각 절점의 변위는 유한요소 해로부터 보간(interpolation)을 통해 구해진다. 먼저, 각 가상 격자 절점을 포함하는 요소를 결정한 후, 해당 요소의 자연좌표계(natural coordinate system)에서의 위치를 등매개변수 사상(isoparametric mapping)의 역변환을 이용하여 계산한다. 이후, 사면체 요소의 형상함수(shape function)를 이용하여 가상 격자 절점에서의 변위를 보간한다. 이때 보간된 변위장은 다음 식 (1)과 같이 표현될 수 있다.

여기서, ui는 유한요소 해에서의 i번째 절점 변위 벡터를 나타내며, Ni는 HEX8 요소의 Lagrange 형상함수이고, ne,n은 요소 내 절점의 개수를 의미한다. 가상 격자 상에서 변위장이 재구성된 이후, 변형률 및 응력장은 형상함수를 미분하여 구조화된 격자 내에서 계산된다. 가상 격자는 형상이 양호한 요소들로 구성되어 있으므로, 수치 미분은 왜곡된 사면체 요소에서 직접 수행되는 경우에 비해 보다 안정적이고 정확하다.

이 접근법은 변위 해와 응력 계산을 효과적으로 분리하여, 후처리 개념으로 다양한 환경에서 수행될 수 있도록 한다. 그 결과, VGSR 방법은 요소망 왜곡 및 적분영역과 관련된 수치 오차를 감소시키고 응력장의 품질을 향상시키며, 이는 영역 기반(domain-based) 방법을 이용한 응력확대계수와 같은 파괴 파라미터의 정확한 평가에 중요하다.

2.2 적응적 요소 세분화

적응적 요소 세분화(AMR)는 균열 선단 부근에서 응력장과 변위장의 해상도를 향상시키기 위해 적용된다. 세분화는 Abaqus의 자동 요소 재구성 기능을 이용하여 수행되며, 이때 von Mises 응력이 세분화 지표로 사용된다(Simulia, 2014). 응력 구배가 큰 영역은 선택적으로 세분화되고, 나머지 영역은 비교적 조밀도가 낮은 상태로 유지되어 계산 효율성을 확보한다. 본 연구에서는 중요한 영역에서 요소망의 품질을 점진적으로 향상시키기 위해 두 단계의 반복적 세분화를 적용하였다.

이 과정의 예시가 Fig. 2에 도시되어 있다. Fig. 2(a)는 임의 형상의 균열을 포함하는 2차원 해석 영역을 나타낸다. Fig. 2(b)는 초기 유한요소망을 보여주며, 이를 이용하여 1차 해석이 수행된다. 이후 계산된 von Mises 응력 분포를 기반으로 적응적 재요소화가 적용되며, 그 결과 응력 집중이 크게 나타나는 영역, 특히 곡률이 큰 균열 선단 부근에서 요소가 더욱 조밀하게 세분화된 요소망이 Fig. 2(c)에 나타나 있다. 요약하면, 복잡한 균열 형상에 대해 생성된 초기 해석 결과로부터 응력 집중 구간을 파악하고 적응적 재요소화를 통해 요소망을 국부적으로 세분화함으로써 해석의 정밀도와 요소망 생성의 높은 유연성과 범용성을 확보할 수 있다.

Fig. 2.

Illustration of the adaptive mesh refinement procedure: (a) two-dimensional domain with an arbitrarily shaped crack, (b) initial finite element mesh used for the first analysis, and (c) adaptively refined mesh based on the von Mises stress distribution, showing increased mesh density near the crack front

3. I-적분(Interaction integral)

I-적분 방법을 이용하여 응력확대계수를 계산하였으며, I-적분은 J-적분(J-integral)과 보조장(auxiliary field)의 중첩을 이용하여 유도되었다(Shih and Moran, 1986). 균열 선단을 따라 위치 s에서의 J-적분은 다음과 같이 정의된다:

여기서, W는 변형에너지 밀도(W=σijεij/2)이며, δij는 Kronecker delta를 나타낸다. 또한 σij, εij, uj는 각각 응력, 변형률, 변위 성분을 의미한다. 라틴 첨자 i, j는 1부터 3까지의 값을 가지며, 반복되는 첨자는 합산을 의미한다. (∙),i=∂(∙)/(∂ξi)는 (∙)를 ξi 방향에 대해 편미분한 것이다. 일반성을 유지하기 위해 평면 및 비평면 균열 모두에 적용 가능한 곡선 좌표계(ξ1, ξ2, ξ3)가 도입된다(Gosz and Moran, 2002)(Fig. 3 참조). 이 좌표계는 균열 선단 형상을 기준으로 정의되며, ξ1은 균열 진전 방향, ξ2는 균열면의 법선 방향, ξ3는 ξ1과 ξ2의 외적을 통해 정의된다. 적분 경로 𝛤는 국부 곡선 좌표계의 ξ1-ξ2 평면에 위치하며, 균열 하부 면에서 시작하여 상부 면까지 확장된다. ni는 적분 경로 𝛤에 대한 단위 법선 벡터의 성분을 나타낸다.

Fig. 3.

(a) Closed integral on the surface perpendicular to the crack front, (b) tabular domain surrounding a small segment of a curved crack front, and (c) interaction integral domain and boundary surfaces defined within the virtual domain (12 × 12 × 4 elements) used in the VGSR method

Shih와 Moran(1986)은 발산정리(divergence theorem)를 이용하여 식 (2)를 수치적 평가를 위한 영역 기반 I-적분 형태로 재정식화하였다(식 (3)).

연속 스칼라 가중 함수(q)는 외부 경계에서 0이고 균열 선단에서 1의 값을 갖는다. 이 함수는 Fig. 4에 나타난 바와 같이 가상 영역에서 피라미드 형상을 갖도록 설정된다. Fig. 4의 q 함수 정의는 Fig. 3(c)에 제시된 VG 기반 I-적분 영역과 연계되어 정의된다. 체적 V는 폐곡선 적분 영역으로서, Fig. 3(b)에 나타난 바와 같이 상부 및 하부 균열면(S+, S-)과 관형 영역의 외부, 내부, 좌측 및 우측 경계면(S1, St, S2, S3)에 의해 둘러싸인다. 본 연구에서는 이러한 경계면들이 Fig. 3(c)에 나타낸 가상 격자 내부의 적분 영역의 경계와 일치하도록 정의되었다.

Fig. 4.

Weight function q on virtual domain consist of 12 × 12 × 4 elements (a) pyramid shape in ξ1-ξ2 (b) q function along the ξ3 direction

J-적분 식에 실제장과 적절한 보조장(Williams, 1957)을 중첩함으로써 I-적분을 정의할 수 있으며, 이때 보조장은 균열 선단 근방의 2차원 점근 응력장 및 변위장을 제공한다. 따라서 I-적분은 실제 해와 보조 해의 선형 중첩을 기반으로 구성된다. 이에 따라 중첩된 I-적분의 식은 다음 식 (4)와 같이 표현된다:

여기서, J¯(s)와 J¯aux(s)는 각각 실제장과 보조장에 대한 J-적분을 의미한다. 또한 I¯(s)는 다음 식 (5)와 같이 표현된다:

편미분 (∙),i는 직교 곡선 좌표계에서 정의된 그래디언트 연산자를 이용하여 계산된다. 이후 위치 s에서의 I-적분 값은 다음 식 (6)과 같이 계산된다:

여기서, L은 관형 영역의 길이를 의미한다. 최종적으로 에너지 방출률은 혼합모드 응력확대계수와 다음과 같은 관계를 가진다:

여기서, E*=E/(1-ν2)이다. 식 (7)에 보조 응력확대계수(KIaux, KIIaux, KIIIaux)를 중첩하면, 실제장과 보조장 간의 I-적분 관계는 다음 식 (8)과 같이 주어진다:

각 모드의 응력확대계수를 분리하기 위해 보조장은 한 번에 하나의 모드만 활성화되도록 [KIaux,KIIaux,KIIIaux]를 각각 [1,0,0], [0,1,0], [0,0,1]과 같이 선택한다.

4. 원형 균열해석 결과

제안된 방법의 검증을 위해 경사진 원형 균열(penny-shaped crack)을 포함하는 3차원 혼합모드 파괴 문제를 고려하였다. 문제의 기하 형상은 Fig. 5(a)에 나타나 있다. 반지름 a를 갖는 원형 균열은 반지름 R을 갖는 원통형 시편의 중심에 위치하며, 반지름 비는 a/R=0.2이다. 균열면은 하중 방향에 대해 𝛽 = 30°의 각도로 경사져 있으며, 이로 인해 균열 선단을 따라 혼합모드 조건이 발생한다.

Fig. 5.

(a) Geometry representation of inclined penny shaped crack (b) Farfield mesh corresponding to the initial discretization, which remains similarly coarse in the adaptive refinement cases

원통형 시편의 상부 및 하부 면에는 균일한 수직 응력 σ0가 가해지며, 이로 인해 모드 I, 모드 II, 모드 III가 결합된 파괴 모드가 발생한다. 균열의 경사로 인해 각 모드의 상대적 기여도는 균열 선단을 따라 위치에 따라 달라진다.

계산 영역은 비정형 3차원 사면체 요소를 이용하여 이산화된다(Fig. 5(b)). 정확한 응력 평가를 위해 VGSR 방법론에 따라 균열 선단을 따라 구조화된 가상 격자가 구성된다. 본 연구에서는 12 × 12 × 4 가상 격자를 사용하며, 가상 요소의 크기는 유한 요소 크기의 2배로 (lVG≈2lFE) 설정하였다. 참고로, 가상 격자의 요소 크기는 주변 유한요소 크기를 기준으로 정의되므로, 국부적인 격자 세분화가 수행되는 경우 이에 대응하여 가상 격자의 크기 역시 함께 감소하도록 설정하였다.

검증을 위해 계산된 응력확대계수는 σ0πa로 정규화되며, 무한 영역에 대한 해석해와 비교된다. 혼합모드 응력확대계수에 대한 해석식은 다음 식 (9), (10), (11)과 같이 주어진다:

초기 요소망은 균일한 분포를 갖는 비정형 사면체 이산화를 이용하여 생성되었으며, 균열 선단을 따라 약 60개의 요소로 구성된다. 이 구성은 비교적 조밀도가 낮은 요소망을 나타내며, 적응적 세분화의 영향을 평가하기 위한 기준 모델로 사용된다.

von Mises 응력을 기반으로 한 1차 적응적 요소 세분화 이후, 균열 선단 주변 영역에서 요소망이 크게 세분화되며, 멀리 떨어진 영역은 상대적으로 조밀도가 낮은 상태를 유지한다. 이로 인해 균열 선단을 따라 총 172개의 절점이 생성된다. 2차 세분화 단계에서는 동일한 영역에서 요소망의 해상도가 더욱 증가하여 균열 선단을 따라 426개의 절점이 형성되며, 원거리 요소망은 거의 변화하지 않는다. 이러한 점진적 세분화는 관심 영역에서 계산 비용을 크게 증가시키지 않으면서도 균열 선단 부근의 급격한 응력 구배를 효과적으로 포착할 수 있게 한다.

각 요소망 구성에 대한 전체 절점 수와 요소 수는 Table 1에 정리되어 있으며, 이는 각 세분화 수준에 따른 계산 비용을 정량적으로 비교할 수 있도록 한다. Fig. 6는 세 가지 경우에 대한 균열면 이산화를 나타내며, 세분화 단계가 진행됨에 따라 균열 선단 부근에서 요소 밀도가 증가하는 것을 보여준다. 반면, 균열 선단에서 멀리 떨어진 영역이 조밀도가 낮은 요소로 이산화되어 있으며, 적응적 요소 세분화 과정을 거친 이후에도 대부분 세분화되지 않은 상태를 유지한다.

Fig. 6.

Mesh distribution on the crack surface at different refinement stages: (a) initial mesh (60 nodes along the crack front), (b) after first adaptive refinement (172 nodes along the crack front), and (c) after second adaptive refinement (426 nodes along the crack front). The refinement is driven by the von Mises stress, resulting in a progressive increase in element density near the crack front

응력확대계수는 VGSR 방법과 I-적분을 결합하여 평가된다. 대칭성을 고려하여 결과는 균열 선단의 절반에 대해서만 제시하였으며, 각도 좌표 𝜃는 0°에서 180°까지의 범위를 갖는다. 초기 요소망과 1차 세분화, 2차 세분화에 따른 해당 응력확대계수 분포의 결과는 각각 Fig. 7에 제시되어 있다. 요소망의 밀도가 증가함에 따라 결과의 진동이 뚜렷하게 감소하고 균열 선단을 따라 보다 매끄러운 분포를 보이며, 이는 수치적 안정성과 수렴성이 향상되었음을 나타낸다.

Fig. 7.

Variation of normalized mixed-mode stress intensity factors along half of the crack front for (a) initial mesh, (b) first refinement, and (c) second refinement. The angular coordinate 𝜃 ranges from 0° to 180°. The results are plotted using a representative subset of crack front nodes rather than all available data points

적응적 요소 세분화가 계산 결과의 정확도에 미치는 영향을 정량적으로 평가하기 위해 각 요소망의 구성에 대해 균열 선단을 따라 모드 I 응력확대계수의 평균값을 계산하였다. 결과는 Table 2에 정리되어 있으며, 정규화된 평균 응력확대계수와 함께 기준 해에 대한 오차가 제시되어 있다. 또한, 균열 전선을 따라 계산된 점별 응력확대계수의 최소값과 최대값도 표에 포함되어 있다. 결과에 따르면 1차 적응적 세분화 이후 오차가 감소하면서 정확도가 크게 향상됨을 확인할 수 있다. 또한 2차 세분화 이후에는 계산된 평균 응력확대계수가 기준 해와 매우 잘 일치하는 결과를 보인다. 이러한 경향은 적응적 요소 세분화가 응력확대계수 평가의 정확도와 수렴성을 향상시키는 데 효과적임을 명확히 보여준다.

적응적 세분화 기법이 정확도 및 계산 효율성에 미치는 영향을 검토하기 위하여 균일 세분화 모델과 적응적 세분화 모델을 비교하였다. 두 모델의 절반 단면 요소망 형상은 Fig. 8에 제시하였다. 약 87개의 균열 선단 노드에 해당하는 균일 세분화 모델은 약 4,101,400개의 요소를 포함하는 반면, 1차 적응적 세분화 후 172개의 균열 선단 노드를 갖는 적응적 세분화 모델은 총 62,861개의 요소만을 필요로 하였다. 균일 세분화 요소망으로부터 얻어진 SIF 분포는 Fig. 9에 나타내었으며, Fig. 7(b)에 제시된 적응적 세분화 결과와 비교하였다. 적응적 세분화 기법은 훨씬 적은 요소 수를 사용함에도 불구하고, 급격한 응력 구배와 특이 응력장이 존재하는 균열 선단 부근에만 조밀한 요소망을 집중적으로 적용함으로써 더욱 정확한 SIF 분포를 제공하였다.

Fig. 8.

Half-sliced mesh configurations of the fracture models: (a) uniformly refined mesh and (b) adaptively refined mesh after the first refinement step

Fig. 9.

SIF distributions obtained using the uniformly refined mesh

5. 결 론

본 연구에서는 3차원 혼합모드 응력확대계수 평가의 정확도를 향상시키기 위해 VGSR 기법과 적응적 요소 세분화를 결합한 해석 프레임워크를 제안하였다. VGSR 방법을 통해 유한요소망과 독립적인 신뢰성 있는 응력장을 확보하였으며, von Mises 응력을 기반으로 한 적응적 요소 세분화를 적용하여 균열 선단 부근의 변위장 품질을 향상시켰다.

제안된 방법은 혼합모드 하중을 받는 경사진 원형 균열 문제에 적용되었으며, 두 단계의 적응적 요소 세분화 전략을 통해 응력확대계수에 미치는 영향을 분석하였다. 그 결과, 요소 세분화는 균열 선단 부근의 응력장 해상도를 향상시키고, 균열 선단을 따라 보다 일관된 응력확대계수 분포를 제공함을 확인하였다. 모드 I 응력확대계수의 평균값을 기반으로 한 정량적 비교에서는 세분화 수준이 증가함에 따라 오차가 명확히 감소하는 경향이 나타났으며, 이를 통해 제안된 방법의 수렴성을 확인할 수 있었다.

종합적으로, 응력 복원 기법과 적응적 요소 세분화 기법의 결합은 복잡한 3차원 균열 문제의 파괴 해석에 있어 효율적이면서도 정확한 접근법을 제공함을 보여준다. 제안된 방법은 향후 보다 복잡한 기하 형상 및 하중 조건으로 확장될 수 있다.