경계면 손상을 고려한 적층복합재료에 대한 멀티스케일 피로 손상 모델

Multi-scale Progressive Fatigue Damage Model for Unidirectional Laminates with the Effect of Interfacial Debonding

Abstract

본 논문에서는 복합재료의 섬유와 기지사이의 경계면 손상을 고려한 멀티스케일 점진적 피로 손상 모델을 제안한다. 먼저 점진적인 경계면 손상을 고려하기 위해 서로 다른 4개의 경계면 상태를 정의한 미소구조 모델을 도입하였다. 각각의 상태에 대한 부피분율은 피로 하중의 사이클 수가 증가함에 따라 온전한 상태의 계면에서 완전 박리 상태의 계면으로의 전환이 일어난다. 손상된 경계면의 에쉘비 텐서(Eshelby's tensor)를 계산하기 위해 선형 스프링 모델이 사용되었으며 균질화 방법을 통해 복합재료의 유효 물성을 얻었다. 또한 복합재료의 피로거동을 묘사하기 위해 교번 응력에 대한 섬유, 기지, 그리고 섬유-기지 간의 계면 각각에 대한 손상 변수들이 정의되었고 이를 chaotic firefly 알고리즘을 통해 손상 변수를 특성화 하였다. 제안된 모델은 유한요소해석프로그램 ABAQUS의 UMAT subroutine으로 구현되어 AS4/3501-6 복합재료의 단일방향 라미네이트(unidirectional laminate) 시편들([0]₈, [90]₈,[30]₁₆)을 통해 성공적으로 검증되었다.

This paper presents a multi-scale progressive fatigue damage model incorporating the model for interfacial debonding between fibers and matrix. The micromechanics model for the progressive interface debonding was adopted, which defined the four different interface phases: (1) perfectly bonded fibers; (2) mild imperfect interface; (3) severe imperfect interface; and (4) completely debonded fibers. As the number of cycles increases, the progressive transition from the perfectly bonded state to the completely debonded fiber state occurs. Eshelby's tensor for each imperfect state is calculated by the linear spring model for a damaged interface, and effective elastic properties are obtained using the multi-phase homogenization method. The fatigue damage evolution formulas for fiber, matrix and interface were proposed to demonstrate the fatigue behavior of CFRP laminates under cyclic loading. The material parameters for the fiber/matrix fatigue damage were characterized using the chaotic firefly algorithm. The model was implemented into the UMAT subroutine of ABAQUS, and successfully validated with flat-bar UD laminate specimens ([0]₈,[90]₈, [30]₁₆) of AS4/3501-6 graphite/epoxy composite.

1. 서론

복합재료는 항공 및 방산 산업 분야를 중심으로 다양한 산업에서의 활용도가 증가되고 있다. 재료의 높은 무게 대비 강성을 가진다는 특징으로 인해 기존의 사용되었던 금속재료를 대체해가고 있다. 하지만 구조물에 복합재료의 적용을 위해서는 재료의 특성을 파악하는 과정이 수반되어야 한다. 특히, 장시간의 반복되는 하중에 노출되는 구조물의 경우 피로특성에 대한 연구가 필수적이다. 오랜 기간 동안 복합재료의 피로특성은 대부분 실험을 통해 연구가 진행되어 왔다. 하지만 복합재료는 등방성을 가지는 금속재료 등과 달리 피로 하중의 방향 및 적층 구조에 따라 서로 다른 피로 파손 거동을 보이며 섬유 파손, 기지 파손, 섬유-기지 간의 박리현상 등의 다양한 손상 메커니즘이 동시에 발생한다. 이를 모두 고려하여 실험을 수행하기에는 많은 시간 및 비용이 요구되기에 이를 대체하기 위한 피로수명 예측 모델에 대한 연구도 같이 진행되어 왔다.

일반적인 피로수명 예측 모델은 크게 ‘fatigue life model’, ‘Phenomenological model’, ‘Progressive damage model’로 구분된다. 먼저, fatigue life models은 S-N 커브 및 Goodman diagram을 이용하여 피로수명을 예측하는 방법으로 Hashin and Rotem[1]에 의하여 처음으로 제안되었다. 그 이후 Reifsnider and Gao[2]가 Mori-tanaka[3] 균질화 방법을 사용하여 미소구조에서 피로 파손기준을 제안하였다. Jen and Lee[4]에서는 Tsai-Hill 파손 기준을 피로 파손기준으로 확장하여 다축 피로하중이 작용하는 평면 응력 상태의 환경에 적용하였다. 하지만 이 방법은 실제 반복하중에서 재료의 손상이 어떻게 진행되는지 설명하지 못한다는 단점이 존재하며 방대한 수의 실험을 필요로 한다.

위 방법의 한계를 극복하기 위하여 피로하중이 진행됨에 따라 재료의 강성 및 강도가 점차 감소하는 것을 피로하중 사이클 수에 대한 함수로 모델링하는 ‘phenomenological model’에 대한 연구가 진행되었다. 이 모델은 크게 강성 저하(stiffness degradation) 모델과 강도 저하(strength degradation) 모델로 구분되는데 강성 저하 모델의 경우 흔히 손상 변수라 불리는 스칼라 변수 D를 도입하여 피로 사이클에 대한 손상 변수의 기울기 dD/dN을 측정한 가능한 물성 값들을 사용하여 정의하여 강성 저하를 설명한다. 하지만 이 모델에서 정의된 손상 변수는 실제 재료의 손상 메커니즘을 반영한 것은 아니라는 점에서 서로 다른 재료에의 적용이 어렵다. Whitworth[5]는 탄소섬유 복합재료에 대하여 강도 저하 모델을 제시하였고 재료의 피로 파손이 피로 하중의 최대 응력이 잔류 강도와 같은 경우 발생한다고 가정하여 피로수명을 예측하였다. Yao and Himmel[6]은 실제 재료가 피로 하중의 초반과 종반 단계에서 급격한 강도 저하, 중반 단계에서 완만한 강도 저하를 보이는 현상을 반영한 모델을 고안하였다. Khan et al.[7]은 직조 탄소섬유 복합재료에 대하여 손상 변수를 사용한 강성저하 모델을 적용하였다. 이 모델에서는 피로 수명은 피로 하중의 최대 변형률이 정적 극한 변형률(static ultimate strain)와 같은 시점이 피로 사이클로 정의하였다.

위의 모델은 피로 실험 결과를 통해 얻은 경험식 기반으로 제안되었기 때문에 실제 재료의 손상 메커니즘을 고려한 피로 수명 예측 모델의 필요성이 대두되었다. 이를 반영하기 위해 유한요소 기반의 점진적 손상 모델이 등장하였다. Shokrieh and Lessard[8,9]가 CFRP 복합재료의 하중에 따라 잔류 강성을 측정하여 이를 하중 사이클 수와 작용하는 응력에 대한 함수로 모델링하였고 이를 유한요소해석에 적용하여 적층형 복합재료에 대하여 검증을 수행하였다. 점진적 손상 모델은 크게 세가지의 단계로 구성되어 있다: (1) 응력 해석, (2) 파손기준 적용, (3) 재료 물성 저하. 먼저 유한요소해석을 통해 FE 모델에 작용하는 응력분포를 추출하고 특정 영역에서 피로하중에 의한 파손이 발생하는지 파손기준을 적용한다. 그 이후 파손여부에 따라 서로 다른 재료 물성저하 모델을 적용하여 피로하중에 의한 재료의 손상이 증가하는 것을 구현하였다. 위의 방법은 AS4/3501-6의 평면 및 핀 고정 시편에 대하여 검증되었다. 이후 Papanikos[10]는 강성저하는 선형, 강도저하는 2차식으로 모델링 하였고 Ye delamination criterion[11] 및 Hashin 3D criterion[12]을 이용하여 시편의 피로수명을 성공적으로 예측하였다. Kenndey et al.[13]은 수정된 Puck 파손 기준을 적용하여 파손 모드에 따른 서로 다른 강도 및 강성 저하 식을 적용하였다. 이를 GFRP 단위셀(unitcell) 모델에 적용하여 피로해석을 수행하였다.

최근에는 피로하중에 의한 손상을 복합재료의 각 구성 재료 단계에서 진전되는 것을 묘사하고자 멀티스케일 접근방식을 이용한 방법들이 제시되었다. Zhang et al.[14]에서는 미소구조 모델 중 하나인 3D General Method of Cells(GMC) 모델을 woven SiC_f/SiC 재료의 피로거동을 예측하기 위해 적용하였다. 이 때 서로 다른 피로 손상 모델 및 파손 기준이 섬유와 기지 각각에 대해 적용되었다. 그리고 Sayyidmousavi et al.[15]은 다양한 적층 구조를 가진 AS4/PEEK의 피로 수명 해석을 위해 Simplified Unit Cell Micromechanics (SUCM) 모델을 사용하였다. SUCM 모델에서는 직사각형 단면을 가지는 섬유와 그 주변을 둘러싸고 있는 기지의 단위셀을 가정한다. 대표체적(RVE)의 응력과 변형율은 각 단위셀에서의 응력과 변형율로부터 계산된다. 피로 파손기준은 섬유와 기지 단위셀의 평균 응력에 대하여 적용되었다. Tamboura et al.[16]은 Mori-Tanaka 방법을 SMC 복합재료의 강성 저하를 예측하는데 이용하였다. 먼저 SMC 복합재료의 섬유 방향 분포를 SEM을 활용하여 얻을 수 있었고 이를 바탕으로 Mori-Tanaka 방법을 사용하여 각 섬유 방향에 따라 섬유의 응력을 계산하였다. 이 때 SMC 복합재료의 경우 섬유-기지 간의 박리가 주요한 파손 메커니즘이므로 박리에 의한 파손 확률이 섬유 응력의 함수로 정의하였다.

이러한 복합재료의 피로 수명 예측을 위한 멀티스케일 접근방법의 적용에도 불구하고 섬유-기지 간의 박리를 고려한 연구는 많이 수행되지 않았다. 하지만 섬유-기지 간의 박리는 적층 복합재료의 주요한 파손 메커니즘 중 하나이므로 필수적으로 고려해야 한다. 따라서 본 연구 연구에서는 섬유-기지 간의 박리를 고려한 인장-인장 피로하중에 대한 멀티스케일 피로 손상 모델을 제시한다. 먼저 서로 다른 계면의 결합 상태를 가지는 네 가지의 섬유 구성성분을 가정하고, 이를 통해 유효 강성을 위해 다상(multi-phase) 손상 모델을 적용하였다. 또한 각 구성재료인 섬유와 기지 그리고 계면에 대한 손상 변수를 도입하여 구성재료 각각의 강성 저하 및 네 가지 섬유 구성성분들의 부피분율의 변화를 구현하였다. 본 연구는 다음과 같이 구성되어 있다. 먼저 2장에서는 피로수명 해석을 위한 멀티스케일 피로 손상모델의 방법론에 대한 정리를 하였다. 다상 손상 모델의 유효 강성을 계산하기 위한 균질화(homogenization) 방법에 대하여 설명하고 피로하중 하에서 각 구성재료에 대한 손상 변수 및 손상변수 기울기에 대한 정의를 소개하였다. 그 이후 3장에서는 2장에서 언급된 모델을 ABAQUS의 user subroutine인 UMAT으로 구현한 뒤 FE simulation을 수행하는 과정, 그리고 단방향 라미네이트 시편([0]₈, [90]₈,[30]₁₆)들에 대해 검증을 진행하는 과정에 대하여 서술하였다. 마지막으로 4장에서는 결론에 대하여 언급하였다.

2. 멀티스케일 피로손상 모델

2.1 미소구조 모델

본 연구에서는 인장-인장 피로하중 사이클이 진행됨에 따라 기지와 섬유의 계면에서 박리 현상이 발생하는 것을 반영하기 위해 Lee and Pyo[17]에서 제안된 다상 손상 모델을 사용하였다. 이 모델에서는 총 다섯가지의 구성성분(phase)이 정의되었다. 초기 상태에서는 기지(phase-0)와 함께 완전히 결합된 상태의 섬유(perfectly bonded fiber, phase-1)가 존재한다. 완전 결합된 상태의 섬유는 연속된 원형 섬유로 가정되며 하중이 가해짐에 따라 섬유와 기지 사이의 계면이 약화되며 완전히 결합된 섬유의 일부는 약한 불완전한 계면(mild imperfect interface, phase-2)를 가지는 섬유로 변화한다. 하중이 더 오래 지속되게 되면 계면의 박리현상이 심화되어 강한 불완전한 계면(severe imperfect interface, phase-3)을 가진 섬유로 변화하게 되고 최종 단계에서는 원통형의 공극(void)으로 간주되는 완전히 박리된(completely debonded, phase-4) 섬유로 변하게 된다. 초기 상태의 복합재료는 공극이 없는 상태로 가정되며 인정한 섬유 간의 상호작용은 무시되었다. 위에서 설명된 다상 손상 모델을 그림으로 나타내면 Fig. 1과 같다.

Fig. 1. Micromechanics model considering the progressive interface debonding

이 미소구조 모델의 유효강성(effective stiffness)을 계산하기 위하여 Ju and Chen[18]의 균질화 방법을 사용하였다. 이균질화 방법에서 섬유들 간의 상호작용을 나타내는 항들은 무시되었고 이는 Mori-Tanaka 균질화 방법과 동일한 결과를 보인다. 유효강성은 식 (1)과 같이 계산된다.

\(\begin{aligned}C^{*}=C_{0} \cdot\left[I+\sum_{r=1}^{4}\left\{\phi_{r}\left(A_{r}+S_{r}\right)^{-1} \cdot\left[I-\phi_{r} S_{r} \cdot\left(A_{r}+S_{r}\right)^{-1}\right]^{-1}\right\}\right]\end{aligned}\)       (1)

A_r ≡ (C_r - C₀)^-1C₀       (2)

식 (1)의 C^*은 유효강성텐서이며 C⁰는 기지의 강성텐서, I는 4차 단위 텐서(identity tensor), ϕ_r은 각 구성성분의 부피분율, 그리고 C_r은 r번째 구성성분의 강성텐서이다. 섬유와 기지는 각각 선형 탄성 직교 이방성 및 선형 탄성 등방성 재료로 가정하였으며 Eshelby’s principle[19]에 의하여 섬유의 내부의 응력과 변형률은 동일하게 분포한다. S_r은 r번째 구성성분의 에쉘비 텐서를 의미하며 완전히 결합된 섬유(S₁)와 완전히 박리된 섬유(S₄)의 에쉘비 텐서식은 Mura[20]에서 얻을 수 있다. S₁ 및 S₄의 식은 (3), (4)와 같다.

\(\begin{aligned}\begin{array}{c}S_{1111}=S_{2222}=\frac{5-4 v_{0}}{8\left(1-v_{0}\right)}, \\ S_{1122}=S_{2211}=\frac{4 v_{0}-1}{8\left(1-v_{0}\right)}\end{array}\end{aligned}\)       (3) 

\(\begin{aligned}\begin{array}{c}S_{1133}=S_{2233}=\frac{v_{0}}{2\left(1-v_{0}\right)}, \\ S_{2323}=S_{1313}=0.25, \\ S_{1212}=\frac{3-4 v_{0}}{8\left(1-v_{0}\right)}\end{array}\end{aligned}\)       (4)

식 (3), (4)에서 ν₀는 기지의 푸아송 비이다.

약한 불완전한 계면과 강한 불완전한 계면을 가지는 섬유의 에쉘비 텐서(S₂, S₃)는 Qu[21]에 의해 제안된 계면에서의 손상이 있는 원통형의 첨가제(inclusion)에 대한 선형 스프링 모델을 적용하여 얻을 수 있다. 계면에 손상이 있는 첨가제의 에쉘비 텐서식은 아래와 같이 얻어진다.

\(\begin{aligned}\begin{array}{c}S_{i j k l}^{M}=\frac{1}{\Omega} \int_{\Omega} S_{i j k l}^{O}(\mathbf{x}) d V(\mathbf{x}) \\ =S_{i j k l}^{O}+\left(I_{i j p q}-S_{i j p q}^{O}\right) H_{p q r s} L_{r s m n}\left(I_{m n k l}-S_{m n k l}^{O}\right)\end{array}\end{aligned}\)       (5)

식 (5)에서 Ω는 균질 선형 탄성 재료의 탄성 영역, L은 4차 탄성 텐서, S⁰는 에쉘비 텐서를 의미한다. 또한 H는 식 (6)으로 주어진다.

\(\begin{aligned}H_{i j k l}=\frac{1}{4 \Omega} \int_{S}\left(\eta_{i k} n_{j} n_{l}+\eta_{j k} n_{i} n_{l}+\eta_{i l} n_{j} n_{k}+\eta_{j l} n_{i} n_{k}\right) d S\end{aligned}\)       (6)

이 때 H는 계면의 물성 및 첨가제의 기하학적 특성에 의해 결정되며 원통형의 첨가제에 대하여 H는 식 (7)과 같이 간단하게 얻어진다.

H_ijkl = αP_ijkl + (β - α)Q_ijkl       (7)

\(\begin{aligned}P_{1111}=P_{2222}=4 P_{2323}=4 P_{1313}=2 P_{1212}=\frac{3 \pi}{8 a}\end{aligned}\)       (8)

\(\begin{aligned}\begin{array}{c}Q_{1111}=Q_{2222}=3 Q_{1122}=3 Q_{2211}=\frac{9 \pi}{32 a}, \text { others }=0\end{array}\end{aligned}\)       (9)

식 (7)~(9)에서 a는 섬유 단면의 반지름을 의미하며, α와 β는 선형 탄성 모델에서 계면의 접선 및 수직 방향의 compliance이다. 선형 탄성 모델에서 strain field를 균일하게 하기 위해서는 α=β이어야 하지만[22,23] 본 연구에서는 실험결과와 잘 일치하는 α, β 값을 사용하였다.

Lee and Pyo[17]의 계산결과를 활용하면 경계면의 손상이 있는 첨가제의 에쉘비 텐서는 식 (10)와 같이 표현할 수 있다.

\(\begin{aligned}\left(S_{n+1}\right)_{i j k l}= & \frac{1}{256 a\left(1-v_{0}\right)^{2}}\left\{S_{I K}^{M(2 n-1)} \delta_{i j} \delta_{k l}\right. \\ & \left.+S_{I J}^{M(2 n)}\left(\delta_{i k} \delta_{j l}+\delta_{i l} \delta_{j k}\right)\right\} \\ & n=1,2\end{aligned}\)       (10)

식 (10)에서 S₂, S₃는 각각 약한 경계면 손상과 강한 경계면 손상을 가지는 첨가제의 Eshelby’s tensor이다. 2차 텐서 S^M(2n-1)_IK, S^M(2n)_IJ 은 식 (11)~(15)와 같이 주어진다.

S^M(2n-1)₁₁ = S^M(2n-1)₂₂ = S^M(2n-1)₁₂ = S^M(2n-1)₂₁

= 32a(-1 + 5v₀ - 4v²₀) - 3πμ₀α_n + 24π(1 - 2v₀)²λ₀β_n + 3π(7 - 32v₀ + 32v²₀)μ₀β_n        (11)

S^M(2n-1)₁₃ = S^M(2n-1)₂₃ = 128αv₀(1 - v₀) + 48π(1 - 2v₀)β_n{(1 - 2v₀)λ₀ - v₀μ₀} n = 1,2       (12)

S^M(2n-1)₃₁ = S^M(2n-1)₃₂ = S^M(2n-1)₃₃ = 0       (13)

S^M(2n)₁₁ = S^M(2n)₂₂ = S^M(2n)₁₂ = S^M(2n)₂₁

= 32α(3 - 7v₀ + 4v²₀) + 3πμ₀(α_n + β_n)       (14)

S^M(2n)₁₃ = S^M(2n)₂₃ = S^M(2n)₃₁ = S^M(2n)₃₂

= 4(1 - v₀)²(16a + 3πμ₀α_n), S^M(2n)₃₃ = 0       (15)

2.2 구성 재료의 피로 손상 진전 모델

이 모델에서는 적층형 복합재료의 피로하중 하에서 손상 메커니즘을 실제 구성재료 단계에서 더 정밀하게 반영할 수 있도록 미소구조단계에서 피로 손상 변수를 정의하였다. 본 연구에서는 복합재료의 delamination 등에 의한 interlaminar 피로손상은 고려하지 아니하였고 intralaminar 손상만 고려하였다. Intralaminar 손상 중 적층형 복합재료가 일반적으로 가지는 세가지의 피로 파손 메커니즘, 섬유파손, 기지 균열, 섬유-기지 간의 박리를 고려하였으며 각각의 파손 메커니즘에 대하여 피로 손상 변수를 정의하였다. 먼저 섬유의 종방향에 대하여 섬유 파손에 대한 손상 변수(D_11f), 기지균열에 대한 기지의 손상 변수(D_m), 마지막으로 섬유-기지 간의 박리에 대한 손상 변수(D_int)를 정의하였다. 각각의 손상 변수를 통해 피로하중 하에서 구성재료의 강성이 저하되는 것을 묘사하기 위해 손상 텐서를 도입하였다. 손상 텐서가 고려된 구성재료의 강성텐서는 아래와 같이 계산할 수 있다.

C_f = (M^-1_f : C_f0 : M^-T_f)       (16)

C_m = (M^-1_m : C_m0 : M^-T_m)       (17)

식 (16), (17)에서 C_f, C_m은 각 구성재료에 대한 손상 변수가 고려된 강성 텐서이고, C_f0, C_m0는 초기상태의 강성 텐서이다. 그리고 M_f, M_m은 손상 텐서이며 식 (18)과 같이 정의된다.

\(\begin{aligned}M_{f}=\left[\begin{array}{cccccc}\frac{1}{\sqrt{1-D_{f 11}}} & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1\end{array}\right]\end{aligned}\)       (18)

\(\begin{aligned}M_{m}=\frac{1}{\sqrt{1-D_{m}}}\left[\begin{array}{llllll}1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1\end{array}\right]\end{aligned}\)       (19)

피로하중이 반복됨에 따라 각각의 손상변수의 진전을 표현하기 위해 일반적으로 많이 사용되는 손상 변수의 기울기를 정의하는 dD/dN approach를 사용하였고 Mohammadi et al.[24]의 power-law 식을 적용하였다. D_11f와 D_m에 대한 손상 변수 기울기 식은 아래와 같다.

\(\begin{aligned} \frac{d D_{k}}{d N}=\frac{A_{k} \Delta Y_{k}^{B_{k}}}{\left(1-D_{k}\right)^{C_{k}}} & =\frac{A_{k}}{\left(2 E_{k}\right)^{B_{k}}} \frac{\left(\Delta \sigma_{k}^{2}\right)^{B_{k}}}{\left(1-D_{k}\right)^{2 B_{k}+C_{k}}} \\ k & =11 f, m\end{aligned}\)       (20)

식 (20)에서 dD_k/dN은 섬유와 기지의 손상 변수 기울기, E_k은 구성재료의 영률, Δσ²_k는 반복하중에 대하여 최대 응력과 최소 응력의 제곱 값의 차이, A_k, B_k, C_k는 각 손상 변수의 파라미터이다.

섬유-기지 간의 박리를 묘사하기 위해 서로 다른 4개의 구성성분에 대한 부피분율을 계산해야 한다. 이들의 부피분율의 분포를 계산하기 위해 섬유와 기지의 손상변수와 유사하게 섬유-기지 간의 박리에 대한 손상 변수 기울기 식을 식 (21)과 같이 정의하였다.

\(\begin{aligned}\frac{d D_{i n t}}{d N}=\frac{A_{\text {int }}}{\left(2 G_{12}\right)^{B_{12}}} \frac{\left(\Delta \sigma_{12}^{2}\right)^{B_{i n t}}}{\left(1-D_{\text {int }}\right)^{2 B_{i n t}+C_{i n t}}}\end{aligned}\)       (21)

이 때, dD_int/dN은 섬유-기지 간의 박리에 대한 손상변수 기울기, G₁₂는 재료의 전단 탄성계수, Δσ²₁₂ 는 반복하중에 대하여 최대 전단응력과 최소 전단응력의 제곱 값의 차이, 그리고 A_int, B_int, C_int는 손상 변수 파라미터이다.

D_int가 증가하면서 강하게 결합된 섬유에서 완전히 박리된 섬유로의 부피분율의 변화가 발생하게 된다. 각 부피분율의 계산을 위해 섬유-기지 간의 박리 발생 확률을 정의하였고 이 확률은 Weibull process를 따른다고 가정하였다[17]. 섬유-기지 간의 박리 발생 확률은 식 (22)와 같고 이를 통해 각각의 구성성분의 부피분율은 식 (23)~(29)으로 계산된다.

Pr_int = 1 - exp(-D_int)       (22)

\(\begin{aligned}\bar{\phi}_{2}=\phi \cdot P r_{i n t}\end{aligned}\)       (23)

\(\begin{aligned}\bar{\phi}_{3}=\bar{\phi}_{2} \cdot P r_{i n t}\end{aligned}\)       (24)

\(\begin{aligned}\phi_{4}=\bar{\phi}_{3} \cdot P r_{i n t}\end{aligned}\)       (25)

\(\begin{aligned}\phi_{1}=\phi-\bar{\phi}_{2}\end{aligned}\)       (26)

\(\begin{aligned}\phi_{2}=\bar{\phi}_{2}-\bar{\phi}_{3}\end{aligned}\)       (27)

\(\begin{aligned}\phi_{3}=\bar{\phi}_{3}-\phi_{4}\end{aligned}\)       (28)

\(\begin{aligned}\bar{\phi}_{2}=\phi \cdot P r_{i n t}\end{aligned}\)       (29)

이 때, 식 (22)의 Pr_int는 섬유-기지 간의 박리 발생 확률, 식(25)~(28)의 ϕ_r은 r번째 구성성분의 부피분율을 의미한다.

2.3 재료 물성 및 손상 변수 파라미터

2.3.1 재료 물성

이 연구에서는 모델 검증을 위한 재료로 AS4/3501-6 적층형 복합재료를 선정하였고 해석을 위한 재료 물성은 기존의 선행 연구로부터 획득하였다. 먼저 lamina 단계에서 강성 및 강도는 Shokrieh et al.[9]의 실험 결과로부터, 기지의 강성 및 강도, 섬유의 강도는 Kaddour et al.[25]에서 실험을 통해 측정한 물성을 활용하였다. 섬유에 대한 강성은 복합재료와 기지의 강성 정보를 이용하여 균질화 기법을 통해 역해석하여 얻었다. 복합재료와 각각의 구성재료 AS4, 3501-6의 물성정보는 Table 1과 Table 2에 정리하였다.

Table 1. Mechanical properties of the composite (AS4/3501-6)

Table 2. Mechanical properties of constituents (AS4. 3501-6)

2.3.2 손상 변수 파라미터 특성화

앞에서 정의된 피로 손상 변수의 파라미터들을 획득하기 위해 유전 알고리즘 기반의 최적화 방법인 chaotic firefly 알고리즘을 사용하였다[26]. 최적화를 위한 목적함수를 구성하기 위해서는 시뮬레이션 데이터와 실험데이터가 필요한데 실험데이터는 Shokrieh et al.[9]에서 측정한 피로 시험중에서의 특정 사이클 수에 대한 잔류 강성을 사용하였다. 잔류 강성 데이터는 종·횡 방향 및 전단 방향 각각에서 대하여 존재하였고 각 방향 별 2개의 서로 다른 하중(종방향: 80,60%, 횡방향: 60,40%, 전단방향: 59%, 40%) 하에서 측정되었다. 잔류 강성 데이터를 통해 식 (30)에서 정의된 특정 사이클 수에서 라미나 단계에서의 손상변수를 계산할 수 있다.

\(\begin{aligned}D_{k}(N)=\frac{E_{k}^{0}-E_{k}(N)}{E_{k}^{0}}\end{aligned}\)       (30)

k = 11(longitudinal), 22(transverse), 12(shear)

식 (30)에서 D_k(N)은 특정 사이클에서의 손상 변수, E⁰_k은 초기상태에서의 강성, 그리고 E_k(N)은 특정 사이클에서의 잔류 강성이다.

시뮬레이션 데이터의 경우, 단일방향 라미네이트의 경우 응력들의 값이 피로 하중 하에서 변화가 작기 때문에 일정하다고 가정하면 손상 변수 기울기 식 (20), (21)들을 사이클 수 N으로 적분하는 과정을 통해 특정 사이클 수 N에 대한 구성재료의 피손 손상 변수를 식 (31)-(33)와 같이 계산할 수 있다.

\(\begin{aligned} D_{11 f}(N)= & 1-\left(1-A_{11 f} *\left(2 B_{11 f}+C_{11 f}+1\right)\right. \\ & \left.* \Delta\left(\sigma_{11}^{f}\right)^{{ }^{B_{11 f}}} * N\right)^{\frac{1}{2 B_{11 f}+C_{11 f}+1}}\end{aligned}\)       (31)

\(\begin{aligned}D_{m}(N)=1-\left(1-A_{m} *\left(2 B_{m}+C_{m}+1\right) * \Delta\left(\sigma_{22}^{m}\right)^{B_{m}} * N\right)^{\frac{1}{2 B_{m}+C_{m}+1}}\end{aligned}\)       (32)

\(\begin{aligned}\begin{array}{l}D_{\text {int }}(N)=1-\left(1-A_{\text {int }} *\left(2 B_{\text {int }}+C_{\text {int }}+1\right)\right. \\ \left.* \Delta\left(\sigma_{12}^{2}\right)^{B_{\text {int }}} * N\right)^{\frac{1}{2 B_{\text {int }}+C_{\text {int } t}+1}} \\\end{array}\end{aligned}\)       (33)

이렇게 계산된 구성재료별 손상 변수를 이용하여 균질화 방법을 통해 라미나 단계에서의 손상 변수를 계산할 수 있다. 시뮬레이션과 실험을 통해 계산된 손상 변수들의 값들을 손상 변수 특성화를 위한 목적 함수 계산에 사용하였다. 이 때 목적함수로는 평균 제곱근 편차(root mean square error, RMSE) 함수를 사용하였다. 위의 과정을 Fig. 2에 나타냈다.

Fig. 2. The algorithm for the fatigue damage parameter characterization

Chaotic firefly algorithm을 통해 얻어진 손상 변수 파라미터들을 Table 3에 정리하였다.

Table 3. Characterizd fatigue damage parameters

3. 멀티스케일 피로 해석 시뮬레이션

3.1 피로 해석 시뮬레이션

멀티스케일 피로 손상 모델을 구조해석에 적용하기 위하여 상용 유한요소해석 프로그램인 ABAQUS의 user subroutine UMAT을 통해 유한요소해석을 진행하였다. UMAT에서는 먼저 솔버에서 획득한 응력과 변형률로부터 각 구성재료별 응력과 변형률을 계산한다. 이를 통해 각 구성재료에서의 손상 텐서를 업데이트한 뒤 강성 텐서를 계산하고 각각의 구성성분들의 부피분율을 적용하여 복합재료의 유효 응력 텐서를 얻을 수 있다. 유효응력 텐서를 통해 새로운 응력을 계산하고 이 결과를 이용하여 손상 변수의 기울기를 새롭게 계산하게 된다. 특정 사이클에서 요소의 손상 변수가 특정 임계 값 이상이 되면 파손이 난 것으로 판단하고 매우 큰 값(0.99)으로 증가시킨다. UMAT의 전체적인 계산과정은 Fig. 3에 간략하게 정리하였다.

Fig. 3. The algorithm for the multiscale fatigue damage model

피로해석의 경우 매 사이클을 모두 고려한 유한요소해석을 수행하기에는 많은 계산시간이 소요된다. 따라서 해석시간을 줄이기 위한 많은 연구가 진행되어 왔는데 이 연구에서는 일반적으로 많이 사용되는 cycle jumping scheme을 사용하였다. Cycle jumping scheme에서는 특정사이클에서 다음 특정사이클로 해석을 건너 뛰는 과정이 존재하게 된다. Cycle jump는 N_jump,n=10^n/m-10^n-1/m으로 주어지며 m은 cycle jump parameter, n은 step number이고 N이 증가할수록 cycle jump가 증가하게 된다. 서로 다른 cycle jump parameter m에 대한 해석시간의 비교는 Fig. 4으로 나타내었다. 해석은 모두 정적하중의 90% 피로하중을 받는 0도 라미네이트 시편에 대하여 수행되었고 500 및 1000 사이클에 도달할 때까지의 시간을 측정하였다. M이 감소함에 따라 해석시간이 빠르게 감소하는 것을 확인할 수 있다.

Fig. 4. Analysis time for different cycle jump parameters

본 모델에서는 해석과정을 총 3개의 단계로 구분할 수 있다. 먼저, 피로하중에 앞서 최대 하중까지 정적 하중이 작용하게 된다(단계1). 그 이후, 피로 하중이 작용하여 최소하중까지 감소했다가 다시 최대 하중에 이르게 된다(단계 2). 이 단계에서는 반복 하중이 작용하면서 응력의 재분배가 이루어지고 요소마다 파손에 대한 해석이 이루어지게 된다. 마지막으로 특정 사이클에서 특정한 값의 사이클만큼 해석을 건너뛰는 cycle jumping 단계(단계3)가 존재한다. 이 과정에서는 손상 변수의 기울기가 일정하다는 가정하에 손상 변수의 값만 새롭게 계산하게 된다. 위의 일련의 과정을 Fig. 4에 도식화했다.

3.2 단일방향 라미네이트 시편에 대한 검증

UMAT subroutine으로 구현된 멀티스케일 피로 손상 모델을 검증하기 위해 단일방향 라미네이트에 대한 피로해석 검증을 수행하였다. 먼저 검증을 위한 실험 데이터로는 Shokrieh et al.[9]의 자료를 사용하였다. 시뮬레이션의 경우 응력비가 0.1인 인장-인장(tension-tension) 피로하중에 대해 진행하였고 각각의 라미네이트에 대한 최대 응력은 Table 4와 같다.

Table 4. The static strengths of undirectional laminates

단일방향 라미네이트 시편의 유한요소 모델은 실험의 조건과 동일하게 길이는 140 mm, 너비는 25.4 mm, 라미나 두께는 0.127 mm로 하였다. 요소로는 C3D8R 요소로 선정하였고 두께방향으로 1개의 요소가 포함되게 하여 총 28000개의 요소로 모델링하였다. 시편 모델의 한쪽 끝에는 고정 경계조건을 부여하였고, 반대쪽에는 반복하중을 가해주었다. 유한요소 모델의 치수와 경계조건은 Fig. 5와 같다.

Fig. 5. The flowchart for FE simulation

총 세 종류의 단일방향 라미네이트에 대한 해석결과 와 실험결과를 Fig. 6과 같이 S-N 선도를 통해 비교하였다. 모든 시편에 대해서 피로 수명의 시뮬레이션 결과가 실험 결과와 유사하게 얻어진 것을 확인할 수 있다. 0도 라미네이트의 경우 응력이 섬유 종방향에 대해서 매우 크게 작용하기 때문에 섬유의 종방향 손상 변수, D_11f에 의하여 파손이 발생하였고, 90도 라미네이트의 경우 기지의 손상변수, D_m에 의하여 파손이 발생하는 것을 확인하였다. 이는 실제 실험의 관찰된 파손 메커니즘과 잘 일치하는 것을 알 수 있다. 30도 라미네이트 또한 90도 라미네이트와 동일한 파손 메커니즘을 가진다. 하지만 이 시편의 경우 전단 응력도 큰 값을 가지기 때문에 섬유-기지 간의 박리에 대한 손상 변수의 기여에 의하여 90도 라미네이트보다 더 작은 피로 수명을 가진다.

Fig. 6. Dimensions and boundary conditions of the FE model​​​​​​​

Fig. 7. S-N curve comparison for flat-bar specimens ([0]₈, [90]₈, [30]₁₆)​​​​​​​

4. 결론

본 연구에서는 적층형 복합재료의 피로수명을 예측하기 위해 멀티스케일 피로 손상 모델을 제시하였다. 이 모델은 섬유-기지 간의 박리 현상에 대한 묘사를 위해 다상 손상모델 및 서로 다른 계면의 결합 상태를 가지는 네 가지의 섬유 구성성분을 가정하고 부피분율들의 변화를 통해 온전한 상태의 계면에서 완전 박리 상태의 계면으로의 전환을 묘사하였다. 각각의 구성재료 및 계면의 손상을 나타내기 위해 섬유의 종방향, 기지, 섬유-기지 간의 박리에 대해 세 가지의 손상변수를 정의하였고 이들의 기울기 모델을 도입하였다. 이 모델은 ABAQUS의 UMAT subroutine으로 구현되었으며 cycle jumping scheme과 함께 사용되어 해석의 정확성을 잃지 않으면서 해석에 소요되는 시간을 줄이고자 하였다. 이 모델은 AS4/3501-6 복합재료의 단일방향 라미네이트 시편들([0]₈, [90]₈, [30]₁₆)을 통해 해석이 수행되었으며 선행연구의 피로 수명과 피로손상을 잘 반영하는 것을 확인하였다. 추후 피로하중 하에서의 박리현상에 대한 실험적 정보가 추가적으로 얻어진다면 더 정밀한 계면의 손상 모델링을 통해 피로 수명예측의 정확성을 향상시킬 수 있으며 다양한 형태의 시편에 적용될 수 있을 것이라 기대한다.

후기

이 연구는 국방과학연구소(20210211EAD-00)와 서울대학교 공학연구원의 지원을 받아 수행된 연구 결과입니다.