Numerical Buckling Analysis of Cracked Graphene-Reinforced Cylindrical Structures using a Phase-Field Crack Formulation

Jin-Rae Cho

doi:10.7734/COSEIK.2025.38.6.355

Preview

Research Paper

Journal of the Computational Structural Engineering Institute of Korea. 31 December 2025. 355-362
https://doi.org/10.7734/COSEIK.2025.38.6.355

Numerical Buckling Analysis of Cracked Graphene-Reinforced Cylindrical Structures using a Phase-Field Crack Formulation

위상필드 균열 정식화를 이용한 그래핀이 강화된 원통형 구조물의 좌굴 수치해석

Jin-Rae Cho¹^*

조 진래¹^*

¹Professor, Department of Naval Architecture and Ocean Engineering, Hongik University, Sejong, 30016, Korea

¹홍익대학교 조선해양공학과 교수

^{*Corresponding Author}

License (open-access, https://creativecommons.org/licenses/by-nc/4.0/):

This is an Open-Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (https://creativecommons.org/licenses/by-nc/4.0) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.

ABSTRACT

This study presents a numerical investigation of the buckling load of a cracked cylindrical panel reinforced with graphene platelets (GPLs). An effective phase-field crack model was incorporated into the framework of the 2-D natural element method (NEM) to accurately evaluate the buckling loads of graphene-reinforced panels. Graphene platelets were introduced into a cracked cylindrical panel according to a prescribed functional distribution pattern throughout the thickness. Conventional crack node separation methods for representing cracks are not only difficult but also prone to numerical inaccuracy and instability. To overcome this limitation, the proposed phase-field model represents crack initiation using a phase-field value without grid manipulation. The numerical stability of the developed crack model was verified through a convergence test with respect to the grid density and comparison with a reference solution. A detailed parametric study was then conducted to examine the influence of key parameters on the buckling load characteristics.

Keywords

graphene-reinforcement

cylindrical composite panel

internal crack

critical buckling load

phase-field formulation

본 논문은 그래핀 혈소판(GPL)으로 보강되고 내부 균열을 가진 원통 복합 구조물에 대한 임계 좌굴하중의 수치적 고찰을 다루고 있다. 임계 좌굴하중을 정확하게 평가하기 위해 이차원 자연요소법(NEM) 기반으로 개발한 효과적인 위상필드 균열모델을 제시하였다. 그래핀 혈소판은 두께방향으로 특정한 기능적 분포패턴으로 원통형 구조물에 삽입되어 있다. 수치적으로 균열의 존재를 표현하는 전통적인 절점분리 기법은 모델링의 번거로움은 물론 수치적 불안정성을 야기할 수 있다. 이러한 문제점을 극복하기 위해 본 논문에서의 위상필드 정식화에서는 수치 그리드의 복잡한 작업 없이 위상 필드를 도입하여 균열을 표현하였다. 개발된 수치기법의 안정성과 신뢰성은 그리드 밀도에 따른 수렴성과 참고문헌과의 비교를 통해 입증하였으며, 그래핀이 보강된 원통 복합재의 좌굴특성을 관련된 주요 인자들에 따른 파라메트릭 수치실험을 통해 고찰하였다.

키워드

그래핀 강화

원통형 복합재

내부 균열

임계 좌굴하중

위상필드 정식화

MAIN

1. 서 론
2. 그래핀 혈소판이 보강된 원통형 구조물
2.1. 그래핀 혈소판의 기능적 보강
2.2. 위상필드 정식화(PFF)
3. 위상필드 균열모델을 이용한 좌굴하중 근사화
3.1. 향상된 전단 변형 쉘이론
3.2. 록킹-프리(locking-free) 자연 요소 근사화
3.3. 가상회전을 이용한 가상필드 근사화
4. 수치실험 및 해석결과
5. 결 론

1. 서 론

최근 들어 그래핀 혈소판(GPL)은 탁월한 물리적 특성으로 전통적인 폴리머 및 금속재료 기반 복합재를 위한 첨단 보강재로서 각광 받고 있다(Shi et al., 2018; Young et al., 2012). 연구에 따르면 균질재료에 그래핀 혈소판을 소량만 첨가하여도 열기계적 특성을 급격히 증가시킬 수 있는 것으로 알려져 있다(Ramanathan et al., 2008). 하지만, 그래핀 혈소판은 현재 상대적으로 높은 생산 비용으로 실제로 보강할 수 함량은 제한적이다. 이러한 가격 문제를 효과적으로 해결하기 위해 소량의 그래핀 혈소판이 기능경사 복합재(FGM)(Cho and Oden, 2000; Shen, 2009) 개념에 따라 구조물의 두께 방향으로 특정한 기능적 분포로 삽입되고 있다. 이에 따라 그래핀 혈소판이 보강된 구조물의 열기계적 특성은 이러한 기능적 분포 패턴에 상당히 좌우되는 것으로 보고되고 있다(Liew et al., 2015).

이에 따라 그래핀 혈소판이 보강된 구조물의 성공적인 실용화를 위해 분포 패턴과 같은 주요 인자들에 대한 파라메트릭 분석이 활발하게 수행되고 있다(Lee and Ghayesh, 2023). 그렇지만 거의 대부분의 연구는 구조물이 균열되지 않은 완전한 상태라는 가정에 기초하고 있다. 하지만, 구조물은 균열에 의해 손상이 되면 구조적 강성이 현저히 감소하게 된다(Seifi et al., 2015). 따라서, 균열이 발생한 구조물의 기계적 거동을 분석하는 것은 구조 안전성 측면에서 대단히 중요하다. 하지만 그리드 기반의 수치해석 기법으로 좌굴과 같은 균열된 구조물의 거동을 평가하는 것은 쉬운 일이 아니다.

그리드 기반의 수치해석에서 균열 모델링은 번잡한 절점분리와 그리드 재구성이 필요할뿐더러, 균열된 비균등 그리드는 수치해석 정확도와 안정성을 저하시킬 수 있다(Cho, 2019). 그리드 기반 수치해석 기법의 이러한 어려움을 해결하기 위한 방안으로 위상필드 정식화(phase-field formulation) 기법이 적용되고 있다(Ambati et al., 2015). 이 기법에서는 절점을 분리하여 균열을 모델링하는 대신 위상 필드(phase field)를 추가적으로 도입하여 구조물의 손상 정도를 표현하여 균열을 모델링하게 된다(Natarajan and Annabattula, 2019). 따라서, 기존의 모델링의 번잡함과 수치해석 정확도와 안정성 저하를 효과적으로 방지할 수 있다.

본 연구는 위상필드 정식화를 활용한 균열 모델링으로 그래핀 혈소판이 보강된 원통형 구조물의 좌굴해석을 내용으로 한다. 위상필드 정식화와 임계 좌굴하중 평가기법은 2차원 무요소기법의 일종인 자연요소법(natural element method)(Cho and Lee, 2006; Sukumar et al., 1998)을 기반으로 개발하였다. 아울러 경사진 균열라인과 그리드 라인의 불일치를 해결하기 위해 경사 균열선의 가상 회전 개념(virtual rotation concept)을 도입하였다. 개발한 수치해석 기법의 안정성과 신뢰성은 그리드 밀도에 따른 수렴성과 참고문헌과의 비교를 통해 입증하였으며, 균열된 원통형 구조물의 임계 좌굴하중을 주요 인자들에 대해 파라메트릭하게 분석하였다.

2. 그래핀 혈소판이 보강된 원통형 구조물

2.1. 그래핀 혈소판의 기능적 보강

Fig. 1은 내부에 경사진 균열이 있는 균일한 두께 $h$ 를 가진 원통형 구조물을 나타내고 있다. 구조물 전체의 영역 𝛺는 중립면(neutral surface) $ϖ$ 와 두께 $h$ 로 $Ω = ϖ \times [- h / 2, h / 2]$ 로 표현된다. 구조물 내부에는 그래핀 혈소판이 그림에서 예시한 바와 같이 네가지 서로 다른 기능적 분포로 삽입되어 있다. 두께방향으로 균일한 분포인 FG-U, 중립면에 밀집된 FG-O, 상하면에 집중된 FG-X, 그리고 아랫면에 집중된 FG-Λ로 구분된다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2025-038-06/N0040380601/images/Figure_jcoseik_38_06_01_F1.jpg

Fig. 1.

A cracked cylindrical panel functionally reinforced with graphene platelets

모재에 삽입되는 보강재의 함량은 일반적으로 단위 체적당 상대체적 비율로 정의되는 체적분율(volume fraction)로 나타낸다. 그래핀 혈소판의 체적분율을 $f_{G P L} (z)$ 로 표시하면, 위 네 가지 기능적 분포는 식 (1)과 같이 표현된다.

(1)

f_{G P L} (z) = \{\begin{cases} V_{G P L}^{*}, FG - U \\ (2 - 4 | z | / h) V_{G P L}^{*}, FG - O \\ 4 | z | V_{G P L}^{*}, FG - X \\ (1 - 2 z / h) V_{G P L}^{*}, FG - Λ \end{cases}

한편, 물리적 관계식으로부터 모재의 체적분율 $f_{m} (z)$ 은 자동적으로 $f_{m} = (1 - f_{G P L})$ 로 계산된다. 여기서, 그래핀 혈소판의 총 체적분율 $V_{G P L}^{*}$ 은 총 질량분율 $g_{G P L}^{*}$ 로부터 결정된다(Yang et al., 2017).

그래핀이 보강된 복합 구조물은 일반적으로 등방성 재료로 모델링되며, 기계적 물성치는 체적분율 기반의 균질화기법(homogenization method)으로 계산된다. 등가 탄성계수 $E_{C}$ 는 잘 알려진 Halphin-Tsai 마이크로 모델(Halphin and Kardos, 1976)을 이용하여 다음과 같이 계산된다.

(2)

E_{C} = \frac{3}{8} \frac{1 + ξ_{L} η_{L} f_{G P L}}{1 - η_{L} f_{G P L}} E_{m} + \frac{5}{8} \frac{1 + ξ_{T} η_{T} f_{G P L}}{1 - η_{T} f_{G P L}} E_{m}

위 식에서 $η_{L}$ 과 $η_{T}$ 는 그래핀 혈소판의 기하학적 치수와 연관되며 다음과 같이 정의된다.

(3)

η_{L} = \frac{E_{G P L} - E_{m}}{E_{G P L} + ξ_{L} E_{m}}, η_{T} = \frac{E_{G P L} - E_{m}}{E_{G P L} + ξ_{T} E_{m}}

여기서, $ξ_{L}$ 과 $ξ_{T}$ 는 그래핀 혈소판의 길이 $ℓ_{G P L}$ , 폭 $w_{G P L}$ , 두께 $t_{G P L}$ 에 의해 $ξ_{L} = 2 ℓ_{G P L} / t_{G P L}$ 그리고 $ξ_{T} = 2 w_{G P L} / t_{G P L}$ 로 계산된다. 한편, 등가 프와송비와 등가 밀도는 단순 선형혼합법칙(Cho and Oden, 2000)으로 계산된다.

2.2. 위상필드 정식화(PFF)

Fig. 2는 경사각 𝛼 그리고 상대 균열길이 $a / S$ 인 내부 경사 균열을 보여주고 있다. 원통형 구조물은 원주방향으로 이동이 가능하고 압력 $p$ 를 받고 있으며, $r_{11}$ 과 $r_{22}$ 는 0과 1사이의 값을 나타낸다. 균열 선단의 응력 특이성(singularity) 표현을 위해 극좌표계( $r, θ$ )를 각 선단에 도입하였다. 그리드 기반의 수치해석에서 절점을 분리하여 균열을 모델링하는 것은 모델링의 번잡함은 물론 수치해석의 불안정 더 나아가 불가능을 초해할 수 있다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2025-038-06/N0040380601/images/Figure_jcoseik_38_06_01_F2.jpg

Fig. 2.

Phase field $ϕ (x)$ within a cracked structure

이러한 문제점은 위상필드 개념(Tobari and Ansari, 2020)을 도입하면 효과적으로 해결할 수 있으며, 위상필드 $ϕ (x)$ 는 구조물의 손상 상태를 나타내는 변수로서 𝜙=1은 완전한 손상을 반면 𝜙=0는 완전히 온전한 상태를 나타낸다. 따라서, 균열 라인은 𝜙=1을 가지는 절점들을 연결함으로서 표한할 수 있고, 균열 라인과 인접한 천이영역(transition region)은 0과 1의 위상값을 가지게 된다. 이 개념에서 구조물의 변형에너지 $U$ 는 다음과 같이 표현된다.

(4)

U = \int_{Ω} (1 - ϕ)^{2} \hat{U} d V

여기서, $\hat{U}$ 은 구조물이 전혀 손상되지 않은 완전한 상태에서의 변형에너지를 나타낸다.

3. 위상필드 균열모델을 이용한 좌굴하중 근사화

3.1. 향상된 전단 변형 쉘이론

원통형 구조물의 변위장 $u = {u, v, w}^{T}$ 은 전단 변형 쉘이론에 따라 다음과 같이 표현된다.

(5)

u = r + z \cdot d

여기서, $r$ 과 $d$ 는 각각 $r = {\{u_{0}, v_{0}, w_{0}\}}^{T}$ 과 $d = {\{β_{x}, β_{y}, 0\}}^{T}$ 를 나타낸다. 그리고 균열에 따른 응력의 특이성(singularity)을 나타내기 위한 향상된 변위장 $\hat{u} = {\hat{u}, \hat{v}, \hat{w}}^{T}$ 은 다음과 같이 향상시킨다.

(6)

\hat{u} = r + z \cdot d + \sum_{I = 1}^{2} \sum_{k = 1}^{4} ({\hat{r}}_{I}^{k} + z {\hat{d}}_{I}^{k}) Q^{T}

여기서, $Q = {\{Q_{I}^{k} (r, θ), Q_{I}^{k} (r, θ), R_{I}^{k} (r, θ)\}}^{T}$ 를 나타내고 각각의 특이 함수는 다음과 같이 정의된다.

(7)

\{Q_{I}^{k} (r, θ)\} = {\sqrt{r} \sin \hat{θ}, \sqrt{r} \cos \hat{θ}, \sqrt{r} \sin \hat{θ} \sin (θ), \sqrt{r} \cos \hat{θ} \sin (θ)}

(8)

\{R_{I}^{k} (r, θ)\} = {\sqrt{r} \sin (\hat{θ}), 0, 0, 0}

여기서, $\hat{θ} = θ / 2$ 이다.

3.2. 록킹-프리(locking-free) 자연 요소 근사화

내부 균열이 있는 원통형 구조물의 좌굴거동을 자연요소법(NEM)으로 근사화 하기 위해서 Fig. 3에 도시한 바와 같이 중립면을 M개의 삼각형과 N개의 그리드 점으로 분할한다. 그리고, 식 (6)의 향상된 변위장 $\hat{u}$ 는 다음과 같이 N개의 라플라스(Laplace) 보간함수 $ψ_{J} (x)$ 로 근사화하게 된다.

(9)

{\hat{u}}^{h} = \sum_{J = 1}^{N} (r_{J} + z \cdot d_{J}) ψ_{J} (x, y) + \sum_{I = 1}^{2} \sum_{k = 1}^{4} (\hat{r_{I}^{k}} + z \hat{d_{I}^{k}}) Q^{T}

https://cdn.apub.kr/journalsite/sites/jcoseik/2025-038-06/N0040380601/images/Figure_jcoseik_38_06_01_F3.jpg

Fig. 3.

Geometry transformation between the curved neutral surface and a planar NEM grid

곡면의 중립면 상에서 Laplace 함수를 정의하고 행렬식을 유도하는 수치작업은 어려운 수학적 유도를 요구하기 때문에 Fig. 3에 나타낸 바와 같이 2차원 사각형 형상의 NEM 그리드로 좌표변환 $T_{C} : (x, y) \to (ζ_{1}, ζ_{2})$ 을 하게 된다. 이에 따라 라플라스 보간함수 $ψ_{J} (x)$ 는 이 그리드에서 $(ζ_{1}, ζ_{2})$ 의 함수로 ${\hat{ϕ}}_{J} (ζ)$ 정의되며, 함수의 형태는 그리드 점의 위치에 따라 Fig. 3에 도시하였다.

한편, 선형 좌굴문제는 가상일의 원리에 따라 다음과 같이 표현된다.

(10)

δ U - δ W_{e x t} = 0

여기서, 가상 변형에너지 $δ U$ 그리고 가상일 $δ W_{e x t}$ 은 다음과 같이 정의된다.

(11)

δ U = \int_{- h / 2}^{h / 2} \int_{ϖ} (1 - ϕ)^{2} δ ϵ^{T} σ d A d z

(12)

δ W_{e x t} = \int_{ϖ} δ {\{\begin{array}{l} w_{, x} \\ w_{, y} \end{array}\}}^{T} p [\begin{matrix} r_{11} & 0 \\ 0 & r_{22} \end{matrix}] δ {\{\begin{array}{l} w_{, x} \\ w_{, y} \end{array}\}}^{T} d A

여기서, 𝜙와 $p$ 는 위상필드 그리고 구조물에 원주방향으로 작용하는 압력을 나타낸다.

식 (10)에 자연요소 근사화 (9)에 도입하고 수치적 처리를 거치게 되면 다음의 선형 고유치 행렬방정식에 도달한다.

(13)

[(K_{σ} + \sum_{e = 1}^{M} K_{s}^{e}) - λ K_{G}] (\bar{b}, \bar{r}) = 0

위 식에서 $\{λ_{J}\}$ 와 $\{b_{J}\}$ 는 각각 임계 좌굴하중(critical buckling load)과 좌굴모드(buckling mode)를 나타낸다. 위 식에 포함되어 있는 강성행렬 $K_{σ}, K_{s}^{e}$ 와 기하 강성행렬 $K_{G}$ 는 다음과 같이 정의된다.

(14)

K_{σ} = \int_{- h / 2}^{h / 2} \int_{ϖ} (B + H)^{T} D (B + H) d A d z

(15)

K_{s}^{e} = \int_{- h / 2}^{h / 2} \int_{ϖ} (\hat{B} + H_{s}^{T}) {\hat{D}}_{s} (\hat{B} + H_{s}) d A d z

(16)

K_{G} = \int_{ϖ} {(B_{w} + H_{w})}^{T} p [\begin{matrix} r_{11} & 0 \\ 0 & r_{22} \end{matrix}] (B_{w} + H_{w}) d A

여기서, $B$ 와 $H$ 는 면내 변형율 그리고 특이성 면내 변형율을 정의하는 편비분을 라플라스 보간함수 그리고 특이성 함수의 곱으로 정의되는 행렬, $D$ 는 면내 변형율과 면내응력을 상관짓는 3 × 3 물성치 행렬이다. 그리고 $\hat{B} + H_{s}$ 그리고 ${\hat{D}}_{s}$ 는 면외 전단 변형률과 면외 전단응력에 상응하는 행렬들이다. 한편, $B_{w}$ 와 $H_{w}$ 는 변위성분 $w$ 의 편비분 연산자와 라플라스 보간함수 그리고 특이성 함수에 곱하여 생성되는 행렬이다. 보다 상세한 정의는 문헌(Cho, 2023a)을 참고하기 바란다.

한편, 굽힘변형 지배적인 박판 구조물의 수치해석에서 발생되는 수치해석 정확도의 저하, 즉 록킹현상(locking phenomenon)(Pitkaranta, 1992)을 방지하기 위해 Lee 등(2014)이 제안한 MIT3+ shell 개념을 적용하였다. 이 기법에 따라 자연요소법 근사화로 구한 면외 변형율을 삼각형 요소 내 특정 지점에서 재 보간(re-interpolation)하였다.

3.3. 가상회전을 이용한 가상필드 근사화

위상필드 $ϕ (x)$ 는 구조물 변위장 근사에 사용된 그리드와 라플라스 보간함수를 이용하여 근사하게 된다. 한편, 내부 균열의 균열선은 그리드(grid)의 수평, 수직 혹은 경사선 중 어느 하나와도 일치하지 않을 수도 있다. 이러한 경우에는 균열선 상에 그리드의 절점들이 놓이지 않기 때문에 위상필드로 균열을 표현할 수 없다. 이러한 문제를 해결하기 위해 Fig. 4와 같이 균열선을 가까운 그리드 선과 일치하도록 위상필드 좌표계 $(x_{p}, y_{p})$ 를 다음의 좌표변환 $T_{R}$ 을 이용하여 가상적으로 회전시키게 된다.

(17)

T_{R} : (x_{p}, y_{p}) \to ({\bar{x}}_{p}, {\bar{y}}_{p})

(18)

\{\begin{array}{l} {\bar{x}}_{p} \\ {\bar{y}}_{p} \end{array}\} = \{\begin{array}{l} x_{0} \\ y_{0} \end{array}\} + \bar{R} \cdot \{\begin{matrix} \cos \bar{θ} \\ \sin \bar{θ} \end{matrix}\}, \bar{θ} = \tan^{- 1} (\frac{y_{p}}{x_{p}}) - (α + θ_{G})

여기서, ${\bar{R}}^{2} = {(x_{p} - x_{0})}^{2} + {(y_{p} - y_{0})}^{2}$ 이고 $θ_{G}$ 는 일치시키고자 하는 그리드 선과 균열라인 사이의 각도이다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2025-038-06/N0040380601/images/Figure_jcoseik_38_06_01_F4.jpg

Fig. 4.

Virtual rotation of inclined central crack

4. 수치실험 및 해석결과

앞의 3절에서 소개한 위상필드 정식화와 자연요소법을 기반으로 한 균열된 그래핀-보강 원통형 구조물의 선형 좌굴해석 기법에 대한 알고리듬을 Fortran으로 코딩하였다. 프로그램은 Fig. 5에 나타낸 바와 같이 자연요소 그리드 생성, 라플라스 보간함수 정의, 록킹-프리 전단변형률 보간, 균열 특이함수 정의, 강성행렬 적분, 고유치 행렬식 생성 및 계산으로 구성되어 있다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2025-038-06/N0040380601/images/Figure_jcoseik_38_06_01_F5.jpg

Fig. 5.

A flowchart for the buckling analysis of cracked graphene-reinforced cylindrical structure

균열된 윈통형 구조물을 개상으로 개발한 위상필드 정식화 기반의 선형 좌굴해석 기법의 안정성을 위해 NEM 그리드의 밀도에 따른 수렴성을 분석하였다. 탄성계수 E = 200.0GPa와 프와송비 𝜈=0.3를 가지고 내부 균열을 가진 경계가 고정된(clamped) 원통형 강재 판넬(panel)을 선정하였다. 구조물 그리고 내부 균열의 기하학적 치수는 R = 0.5m, L = S = 1.0m, h = 0.025m이고 균열 상대길이와 각도는 각각 $a / S$ = 0.7, 𝛼 = 45°이다. 원주방향으로의 단축 좌굴(uni-axial buckling; $r_{11}$ = 1.0, $r_{22}$ = 0.0)로 설정하였다.

Table 1에 나타나 있는 임계 좌굴하중은 ${\hat{λ}}_{1} = λ_{1} L^{2} / D$ 로 무차원화하였으며, 여기서 $D = E h^{3} / 12 (1 - ν^{2})$ 이다. 그리고 ${\hat{λ}}_{1}^{r e l}$ 은 그리드 밀도 25 × 25로 계산한 임계 좌굴하중에 대한 상대 차이(relative difference)를 의미한다.

그리드 밀도가 증가할수록 임계 좌굴하중이 점진적으로 줄어들고 특정한 값으로 수렴함을 확인할 수 있다. 그리드 밀도 21 × 21에서 상대 차이가 2% 미만인 수렴성 결과를 바탕으로 나머지 수치실험은 이 그리드 밀도를 적용하였다.

Table 1.

Sensitivity of ${\hat{λ}}_{1}$ of cracked isotropic cylindrical panel to the grid density (L/R = 2, S/L = 1, R/h =20, a/S = 0.7, 𝛼 = 45°)

Items	Grid density
Items	9 × 9	13 × 13	17 × 17	21 × 21	25 × 25
${\hat{λ}}_{1}$	51.01	45.05	42.59	42.01	41.93
${\hat{λ}}_{1}^{r e l} (%)$	21.65	7.42	1.56	0.19	-

다음으로 본 연구의 해석기법의 신뢰성을 검증하기 위해 내부 균열을 가진 경계가 고정된 세라믹 원통 판넬의 임계 좌굴하중을 계산하였다. 참고문헌(Torabi and Ansari, 2020)에서 사용한 재료물성치와 기하학적 치수는 각각 E = 380.0GPa, 𝜈 = 0.26, R = 0.1m, L = S = 0.2m 그리고 h = 0.005m이다. 내부 균열의 상대길이는 a/S = 0.3이고, 좌굴하중은 위 수렴성 예제와 같이 무차원화하였으며, 단 탄성계수와 프와송비는 알루미늄 물성치(E = 70.0GPa, 𝜈 = 0.3)를 사용하였다. Table 2에 균열여부와 균열각도에 따른 결과를 참고문헌과 비교하여 나타내었다. 본 연구를 통해 구한 임계 좌굴하중이 참고문헌에서의 결과와 큰 차이가 없음을 확인할 수 있다.

Table 2.

Comparison of critical buckling loads ${\hat{λ}}_{1}$ of cracked ceramic cylindrical panel (L/R = 2, S/L = 1, R/h = 20, a/S = 0.3)

Method	Intact	Crack angle 𝛼(deg)
Method	Intact	0	45	90
Torabi and Ansari (2020)	896.21	624.52	826.08	857.26
Present	912.06	641.29	809.91	865.75

다음으로 경사진 내부균열을 가진 그래핀이 강화된 원통형 구조물의 좌굴거동을 분석하였다. 기저 재료(matrix material)는 수렴성 예제의 등방성 강재이며 재료 물성치와 기하학적 치수는 앞의 수렴성 예제와 동일하게 설정하였다. 반면, 보강재인 GPL의 재료 물성치는 참고문헌(Cho, 2023b)를 참고하기 바란다. 특별한 언급이 없는 한 원주방향으로의 단축 좌굴(uni-axial buckling; r₁₁ = 1.0, r₂₂ = 0.0)로서 상대 균열길이는 a/S = 0.3이고 GPL의 패턴과 함량은 각각 FG-U 그리고 $g_{G P L}^{*}$ = 0.8%로 설정하였다. Fig. 6은 GPL의 함량 이 FG-U분포로 보강된 경우에 있어 균열각도에 따른 2차 좌굴모드를 보여주고 있다. 좌굴모드가 균열각도에 따라 달라짐을 확인할 수 있다. 참고로 1차 좌굴모드는 균열각도에 따른 차이가 그다지 심하지 않아 2차 모드를 선택하였다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2025-038-06/N0040380601/images/Figure_jcoseik_38_06_01_F6.jpg

Fig. 6.

Second buckling modes:
(a) 𝛼 = 0°, (b) 𝛼 = 30°, (c) 𝛼 = 45°, (d) 𝛼 = 90°

Fig. 7은 서로 다른 균열각도에 있어 GPL의 함량 $g_{G P L}^{*}$ 이 임계 좌굴하중 ${\hat{λ}}_{1}$ 에 미치는 영향을 나타내고 있다. 그림으로부터 임계 좌굴하중이 GPL의 함량에 비례하여 증가함을 확인할 수 있다. 이것은 모재인 강재보다 탄성계수가 상대적으로 큰 GPL의 비중이 증가할수록 구조물의 강성이 증가하기 때문이다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2025-038-06/N0040380601/images/Figure_jcoseik_38_06_01_F7.jpg

Fig. 7.

Variation of ${\hat{λ}}_{1}$ to the GPL mass fraction $g_{G P L}^{*}$ for different crack angles

Fig. 8은 서로 다른 GPL의 보강패턴에 있어 균열 경사각도가 임계 좌굴하중에 미치는 영향을 나타내고 있다. 그림으로부터 임계 좌굴하중은 균열 경사각도와 더불어 증가하다가 𝛼 = 60° 이후부터는 증가세가 둔화되고 있다. GPL의 보강 패턴에 있어서는 FG-X가 가장 높은 임계 좌굴하중을 반면 FG-O가 가장 낮은 임계 좌굴하중을 각각 보여주고 있다. 그리고 FG-U와 FG-Λ는 중간 범위의 임계 좌굴하중을 나타낸다. 이러한 임계 좌굴하중의 상대적인 차이로 부터 GPL의 보강 위치가 구조물의 중립면에서 멀어질수록 구조물의 강성이 상대적으로 증가함을 확인할 수 있다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2025-038-06/N0040380601/images/Figure_jcoseik_38_06_01_F8.jpg

Fig. 8.

Variation of ${\hat{λ}}_{1}$ to the crack angle 𝛼 for different GPL distribution patterns

Fig. 9는 구조물의 서로 다른 상대두께 L/h에 있어 균열 상대길이가 임계 좌굴하중에 미치는 영향을 나타내고 있다. 임계 좌굴하중은 균열의 상대길이와 더불어 감소하다가 특정 상대길이 이후로 부터는 감소세가 둔화되고 있다. 이것은 균열의 길이가 증가할수록 구조물의 강성이 감소하기 때문이다. 한편, L/h의 값은 구조물의 두께 h가 감소할수록 증가한다. 따라서, 동일한 균열 경사각도에 있어 L/h가 작을수록 임계 좌굴하중은 증가하게 된다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2025-038-06/N0040380601/images/Figure_jcoseik_38_06_01_F9.jpg

Fig. 9.

Variation of ${\hat{λ}}_{1}$ to the relative crack length a/S for different structure thickness

다음으로 면내 압축하중이 작용하는 방향에 따른 임계 좌굴하중을 계산하고 Fig. 10에 나타내었다. 그림에서 r₁₁ = 1과 r₂₂ = 1은 각각 원주방향 그리고 축방향으로의 단축 좌굴을 나타내고 r₁₁ = r₂₂ = 1은 양축 좌굴을 의미한다. 그림으로부터 축방향으로의 단축 좌굴하중이 원주방향 단축 좌굴하중보다 현저히 높음을 알 수 있다. 이는 원통형 구조물에서 축 방향으로의 구조강도가 원주 방향으로 값보다 높기 때문이다. 한편, 양축 좌굴강도는 원주방방향 좌굴하중보다 낮게 나오는데 이는 양방향에서 압축하중이 작용하기 때문에 그만큼 좌굴에 필요한 임계하중이 줄어들기 때문이다. 한편, 그래핀 함량에 따른 임계 좌굴하중의 증가추세는 거의 유사함을 알 수 있다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2025-038-06/N0040380601/images/Figure_jcoseik_38_06_01_F10.jpg

Fig. 10.

Variation of ${\hat{λ}}_{1}$ to the buckling axis for different GPL mass fractions.

5. 결 론

본 논문에서는 위상필드 정식화를 기반으로 한 좌굴해석기법으로 내부 균열을 가진 그래핀이 보강된 원통형 구조물의 임계 좌굴하중을 분석하였다. 수치해석 기법은 2차원 자연요소법(NEM)을 기반으로 개발되었으며, 균열선과 그리드 라인과의 불일치를 해결하기 위해 위상필드 좌표를 가상적으로 회전시키는 방법을 도입하였다. 균열에 따른 거동의 특이성을 묘사하기 위해 균열선단 특이함수를 추가하여 변위장 근사를 보완하였으며, 록킹에 따른 해석정도 저하를 해결하기 위해 전단 변형율을 MITC 개념을 활용하여 재보간하였다.

개발된 좌굴해석기법의 신뢰성과 안정성은 참고문헌과의 비교검증 그리고 그리드 밀도에 따른 수렴성을 통해 입증하였다. 그래핀이 강화된 균열된 원통형 구조물의 임계 좌굴하중은 그래핀의 함량에 비례하여 증가하였으며, 그래핀의 보강 패턴에 있어서는 FG-X, FG-U, FG-Λ, FG-O 순으로 감소하였다. 한편, 임계 좌굴하중은 균열각도와 더불어 증가하였고, 균열 상대길이에 대해서는 감소하였다. 하지만, 두 경우 공히 특정 각도와 상대길이 이상에서는 증가 그리고 감소세가 줄어들었다. 그리고 단축 좌굴에 있어 원주방향 보다 축 방향 임계 좌굴하중이 현저히 높았으며, 양축 좌굴은 단축 좌굴에 비해 다소 낮은 임계 좌굴하중을 나타내었다.

Acknowledgements

이 논문은 2025학년도 홍익대학교 학술연구진흥비에 의하여 지원되었음. 이 성과는 정부(과학기술정보통신부)의 재원으로 한국연구재단의 지원을 받아 수행된 연구임(No. RS-2023-00240618, 2020R1A2C1100924).

References

Ambati, M., Gerasimov, T., De Lorenzis, L. (2015) A Review on Phase-Field Models of Brittle of Brittle Fracture and a New Fast Hybrid Formulation, Comput. Mech., 55, pp.383~405.

10.1007/s00466-014-1109-y

Cho, J.R. (2019) Near-Tip Grid Refinement for the Effective and Reliable Natural Element Crack Analysis, Struct. Eng. & Mech., 70, pp.279~287.

Cho, J.R. (2023a) Investigation of Buckling behavior of Cracked FG Cylindrical Panels Reinforced by Graphene Platelets, Symmetry, 15, 2162.

10.3390/sym15122162

Cho, J.R. (2023b) Free Vibration Analysis of Functionally Graded Porous Cylindrical Panels Reinforced with Graphene Platelets, Nanomater., 13, 1441.

10.3390/nano1309144137176986PMC10180390

Cho, J.R., Lee, H.W. (2006) A Petrov-Galerkin Natural Element Method Securing the Numerical Integration Accuracy, J. Mech. Sci. & Technol., 20, pp.94~109.

10.1007/BF02916204

Cho, J.R., Oden, J.T. (2000) Functionally Graded Material: A Parametric Study on Thermal-Stress Characteristics using the Crack-Nicolson-Galerkin Scheme, Comput. Methods Appl. Mech. & Eng., 188, pp.17~38.

10.1016/S0045-7825(99)00289-3

Halphin, J.C., Kardos, J.L. (1976) The Haplin-Tsai Equations: A Review, Polym. Eng. Sci., 16, pp.344~352.

10.1002/pen.760160512

Lee, K., Ghayesh, M.H. (2023) A Review on the Mechanics of Graphene Nanoplatelets Reinforced Structures, Int. J. Eng. Sci., 186, 103831.

10.1016/j.ijengsci.2023.103831

Lee, Y., Lee, P.S., Bathe, K.J. (2014) The MITC3+ Shell Element and Its Performance, Comput. & Struct., 281, pp.12~23.

10.1016/j.compstruc.2014.02.005

Liew, K.M., Lei, Z.X., Zhang, L.W. (2015) Mechanical Analysis of Functionally Graded Carbon Nanotube Reinforced Composites: A Review, Compos. Struct., 120, pp.90~97.

10.1016/j.compstruct.2014.09.041

Natarajan, S., Annabattula, R.K. (2019) Phase Field Modeling of Crack Propagation in Functionally Graded Materials, Compos. Part B, 169, pp.239~248.

10.1016/j.compositesb.2019.04.003

Pitkaranta, J. (1992) The Problem of Membrane Locking in Finite Element Analysis of Cylindrical Shells, Numnerische Math., 61 pp.523~542.

10.1007/BF01385524

Ramanathan, T., Abdala, A.A., Stankovich, S., Dikin, D.A., Herrera-Alonso, M., Piner, R.D., Adamson, D.H., Schniepp, H.C., Chen, X., Ruoff, R.S., Nguyen, S.T., Aksay, I.A., Prud’Homme, R.K., Brinson, L.C. (2008) Functionalized Graphene Sheets for Polymer Nanocomposites, Nature Nanotechnol., 3, pp.327~331.

10.1038/nnano.2008.96

Seifi, R., Googarchin, H.S., Farrokhi, M. (2015) Buckling of Cracked Cylindrical Panels under Axially Compressive and Tensile Loads, Thin-Walled Struct., 94, pp.457~465.

10.1016/j.tws.2015.05.011

Shen, S.H. (2009) Nonlinear Bending of Functionally Graded Carbon Nanotube-Reinforced Composite Plates in Thermal Environments, Compos. Struct., 91, pp.9~19.

10.1016/j.compstruct.2009.04.026

Shi, G., Arby, S., Gibson, C.T., Meng, Q., Zhu, S., Ma, J. (2018) Graphene Platelets and Their Polymer Composites: Fabraication, Structure, Properties, and Applications, Adv. Funct. Mater., 28, 1706705.

10.1002/adfm.201706705

Sukumar, N., Moran, B., Belytschko, T. (1998) The Natural Element Method in Solid Mechanics, Int. J. Numer. Meth. Eng., 43, pp.839~887.

10.1002/(SICI)1097-0207(19981115)43:5<839::AID-NME423>3.0.CO;2-R

Torabi, J., Ansari, R. (2020) Numerical Investigation on the Buckling and Vibration of Cracked FG Cylindrical Panels based on the Phase-Field Formulation, Eng. Fract. Mech., 228, 106895.

10.1016/j.engfracmech.2020.106895

Yang, B., Yang, J., Kitipornchai, S. (2017) Thermoelastic Analysis of Functionally Graded Graphene Reinforced Rectangular Plates based on 3D Elasticity, Meccanica, 52, pp.2275~2292.

10.1007/s11012-016-0579-8

Young, R.J., Kinloch, I.A., Gong, I., Novoselov, K.A. (2012) The Mechanics of Graphene Nano-Composites: A Review, Compos. Sci. Technol., 72, pp.1459~1476.

10.1016/j.compscitech.2012.05.005

Journal of the Computational Structural Engineering Institute of Korea ISSN:1229-3059(Print) 2287-2302(Online) 한국전산구조공학회 논문집

Preview

Numerical Buckling Analysis of Cracked Graphene-Reinforced Cylindrical Structures using a Phase-Field Crack Formulation

ABSTRACT

MAIN

Fig. 1.

A cracked cylindrical panel functionally reinforced with graphene platelets

(1)

(2)

(3)

Fig. 2.

Phase field ϕ(x) within a cracked structure

(4)

(5)

(6)

(7)

(8)

(9)

Fig. 3.

Geometry transformation between the curved neutral surface and a planar NEM grid

(10)

(11)

(12)

(13)

(14)

(15)

(16)

(17)

(18)

Fig. 4.

Virtual rotation of inclined central crack

Fig. 5.

A flowchart for the buckling analysis of cracked graphene-reinforced cylindrical structure

Table 1.

Sensitivity of λ^1 of cracked isotropic cylindrical panel to the grid density (L/R = 2, S/L = 1, R/h =20, a/S = 0.7, 𝛼 = 45°)

Table 2.

Comparison of critical buckling loads λ^1 of cracked ceramic cylindrical panel (L/R = 2, S/L = 1, R/h = 20, a/S = 0.3)

Fig. 6.

Second buckling modes: (a) 𝛼 = 0°, (b) 𝛼 = 30°, (c) 𝛼 = 45°, (d) 𝛼 = 90°

Fig. 7.

Variation of λ^1 to the GPL mass fraction gGPL* for different crack angles

Fig. 8.

Variation of λ^1 to the crack angle 𝛼 for different GPL distribution patterns

Fig. 9.

Variation of λ^1 to the relative crack length a/S for different structure thickness

Fig. 10.

Variation of λ^1 to the buckling axis for different GPL mass fractions.

Acknowledgements

References

Phase field $ϕ (x)$ within a cracked structure

Sensitivity of ${\hat{λ}}_{1}$ of cracked isotropic cylindrical panel to the grid density (L/R = 2, S/L = 1, R/h =20, a/S = 0.7, 𝛼 = 45°)

Comparison of critical buckling loads ${\hat{λ}}_{1}$ of cracked ceramic cylindrical panel (L/R = 2, S/L = 1, R/h = 20, a/S = 0.3)

Second buckling modes:
(a) 𝛼 = 0°, (b) 𝛼 = 30°, (c) 𝛼 = 45°, (d) 𝛼 = 90°

Variation of ${\hat{λ}}_{1}$ to the GPL mass fraction $g_{G P L}^{*}$ for different crack angles

Variation of ${\hat{λ}}_{1}$ to the crack angle 𝛼 for different GPL distribution patterns

Variation of ${\hat{λ}}_{1}$ to the relative crack length a/S for different structure thickness

Variation of ${\hat{λ}}_{1}$ to the buckling axis for different GPL mass fractions.