1. 서 론

최근 예기치 못한 극한의 폭발하중에 노출된 주요 사회기반구조물에 대한 성능 및 안전성 평가를 위해, 실험적 연구의 한계성으로 수치해석적 연구가 활용되고 있다(^{Nam et al., 2008}; ^{Tham, 2009}). 특히 폭발해석의 경우 폭발에 의한 폭풍파의 유체(fluid) 영역과 표적재인 구조체의 고체(solid) 영역이 연성되어 있는 다물리계 특성을 고려해야 하며, 수치해석의 효율성 및 정확한 폭발하중을 확보하기 위해 유체영역에서의 적합한 메쉬 크기가 선정되어야 한다(^{Draganic and Varevac, 2018}; ^{Yun and Park, 2013}). 본 연구에서는 철근 콘크리트(Reinforcement concrete, 이하 RC) building 구조물에서 발생되는 다양한 폭발 시나리오에 따른 폭발해석 및 손상평가가 수행된다. 이를 위해 먼저 폭발하중에 대한 이론적 고찰 및 폭발량에 따른 수치 모델링이 이루어지고, 폭발 시나리오는 TNT 질량 20 kg 및 100 kg 이 각각 RC building 구조물 내 총 7 지점 위치에 적용되며, 구조물의 폭발손상 거동을 효과적으로 모사하기 위해 고체영역의 Lagrange 기법 및 유체영역의 Euler 기법이 커플링된 Euler-Lagrange 커플링 기법이 적용된다. 이 때 폭발하중에 의한 고변형률 속도가 고려된 콘크리트 및 강재 재료구성모델이 제안되며, 이는 ^{ANSYS-AUTODYN solver(2005)}에 연결되어 폭발 수치해석이 수행된다.

2. 폭발 해석

2.1 폭발 하중

폭발에 의해 발생되는 시간에 따른 압력하중의 경우, 평상시 대기압 $P_0$를 유지하고 있지만, 폭발 발생 후 일정시간이 지나면 최대압력 이후 다시 대기압과 같은 압력을 갖게 되는데, 이 구간을 정압기(positive phase duration)라고 하며, Fig. 1의 빗금 친 부분의 면적만큼 충격하중으로 작용하게 된다. 이후 대기압보다 낮은 압력을 갖게 되는 부압기(negative phase duration)를 거쳐 평상시 압력으로 돌아오게 된다. 구조물 설계 시에는 정압단계의 최대압력이 구조물의 거동에 중요한 영향을 미치고, 부압단계에 의한 충격량은 정압단계의 충격량에 비해 상대적으로 작기 때문에 고려하지 않으며, 정압단계의 압력-시간 곡선은 Eq. (1)과 같이 지수함수로 표현할 수 있다(^{Bulson, 1997}). 정압단계의 초과압력을 삼각형의 압력-시간 곡선으로 가정하면 단순화된 시간-압력 곡선식이 Eq. (2)와 같이 얻어진다. 이때 TNT 폭발일 경우 최대압력 $P_{so}^{+}$ 은 Eqs. (3) and (4)와 같이 결정할 수 있다(^{Hetherington and Smith, 1994}).

여기서, $p\left(t\right)$는 시간 $t$에서의 초과압력, $P_0$는 대기압, $P_{so}^{+}$는 최대압력, $T^{+}$는 정압단계의 지속시간, $b$는 압력하강의 적합상수를 의미하며, $T^{-}$는 부압단계의 지속시간, $i_{so}^{+}$과 $i_{so}^{-}$는 각각 양의 충격량과 음의 충격량을 나타낸다. TNT의 양으로 환산한 폭풍파의 환산거리는 $z=r/w^{1/3}$로 계산되고 $r$은 폭발지점에서 측정점까지의 거리(m), $w$는 TNT 양으로 환산한 폭발물의 장약 질량(kg)를 나타낸다.

2.2 지배방정식

폭발 해석은 평형조건에 기초한 유한요소 해석과는 달리 질량, 운동량, 그리고 에너지 보존방정식에 근거한다. 이 때 지배방정식은 관성효과가 고려되어야 하고, 보존방정식은 매 시간단계에서 계산이 일제히 수행되어야 하며, 이는 Table 1에 나타나 있다. $\rho$은 밀도, $u_i$은 속도, ${\sigma }_{ij}$은 Cauchy 응력, $f_i$은 외력, $e$는 내부에너지, $s_{ij}$는 편차응력, $p$는 정수압력, 그리고 $\dot{{ϵ}_{ij}}$는 변형률속도를 나타낸다.

2.3 Temporal Domain (Explicit Method)

Implicit 시간적분법은 Newton-Raphson 반복절차를 통해 다소 큰 시간간격을 적용할 수 있기 때문에 빠른 수렴성을 가지는 장점이 있으나, 외력과 내력의 평형상태 유지를 위하여 매 시간간격에서 시스템 연립방정식을 풀어야 하는 단점이 있다. 상대적으로, 폭발하중과 같이 순간적인 하중작용과 극도의 비선형성을 수반하는 동역학 문제를 해석하는 경우, 평형상태 유지를 위한 시간간격이 매우 짧아지는 특징이 있기 때문에 강성 및 질량 행렬에 대한 계산량이 기하급수적으로 증가되어 해석에 어려움이 있다. 이와 같은 문제를 해결하기 위한 explicit 시간적분법은 충돌과 같은 극도의 비선형성을 수반하는 문제의 경우, 외적인 작용하중과 내부 구조와의 평형을 부여하기 위하여 강성 및 질량 행렬을 매 시간간격마다 계산되는 것이 불필요하므로 적합한 수치해석기법이다. Explicit 수치해석을 위한 유한요소 공간상에 이산화된 운동방정식(discretized equation of motion)이 Eq. (5)에 나타나 있다.

여기서, $M_d$는 질량 행렬이며 explicit 시간적분의 장점을 살리고 CPU 연산처리 속도 향상을 위해 대각화 행렬이 적용되어진다. $f_{\mathrm{ext}}$는 외력벡터, $f_{\mathrm{int}}$는 내력벡터를 의미하고, 6육면체 요소의 1점 적분시 발생되는 모래시계 모드(hourglass mode)를 제거하기 위하여 안정화 벡터 $f_{\mathrm{stab}}$가 추가되어진다. 운동방정식의 시간적분은 explicit 차분방식인 중앙차분(central difference)방식이 적용된다. 즉, 시간이 $t_n$인 상태에서 각 노드에서의 시간적분은 Eq. (6), Eq. (7), Eq. (8)과 같이 수행된다.

또한 충돌문제 해석 시 물리량들의 급격한 변화가 있는 충격면에서 초래되는 수치발산(numerical oscillation) 현상을 제거하기 위하여 불연속인 충격파 전면에 수학적 연속성이 부여되어야 한다(^{von Neumann and Richtmyer, 1950}). 따라서 수치적 안정성을 도모하기 위한 정수압력 항 $q$를 추가하여 해석이 수행되며 이는 Eq. (9)에 나타나 있다.

여기서, $\rho$는 밀도, $V$는 현재 상태의 요소부피, $b_1$과 $b_2$는 사용자 정의 상수, ${\dot{ϵ}}_V$는 체적 변형률 속도, $c$는 음속을 의미한다.

2.4 Spatial Domain

순간적으로 급격히 작용하는 폭발하중에 대한 구조물의 응답을 예측하기 위해 공간이산 방식에 따라 크게 Lagrange 및 Euler 기법으로 나눌 수 있다. 이중 연속체 고체역학에서 일반적으로 사용되고 있는 메쉬에 기반한 Lagrange 기법은 다른 기법들에 비해서 계산 효율성이 뛰어나고, 재료에 고정된 좌표를 통해 시간 변화에 따른 재료거동을 기술할 수 있는 장점이 있어 재료의 변형거동을 조사하는데 적합한 반면 충돌 현상 시 발생할 수 있는 국부적 대변형으로 인한 메쉬 얽힘(tangling) 현상이 발생된다. 이러한 메쉬 왜곡문제는 요소 소진(element erosion)과 같은 수치적 알고리즘을 도입하여 일부 해결이 가능하지만 메쉬의 재정렬이나 소진과정에서 불가피하게 발생하는 메쉬 의존성(mesh-dependency) 문제, 에너지 손실문제, 그리고 접촉(contact)면 불규칙 형상 문제는 수치해의 안정성에 여전히 큰 영향을 미친다(^{Dohrmann et al., 2000}).

Euler 기법은 유한요소 메쉬를 재료의 변형과 무관한 공간에 고정시켜 시간 변화에 따른 재료거동 해석을 수행하므로 고체 영역의 문제보다는 폭발가스와 폭풍파의 전파와 같은 유체 영역의 문제에 적합하다(^{Gebbeken and Ruppert, 2000}). Euler 기법은 기본방정식 내에 이류항(convective terms)이 정식화되어 수치해의 안정성이 저해되는 문제가 있으며, 순간적인 충격 및 충돌 문제를 정교한 Euler 메쉬를 통해 계산해야 할 경우, 미소 시간간격이 적용되어 해의 정밀도를 유지하기 위해 급격한 계산량 증대가 초래되며, 모델링 시 구조화된 메쉬만을 사용해야 하는 단점이 있다.

2.5 Euler-Lagrange Coupling Method

폭풍파의 대기전파, 유체-고체 상호작용, 높은 폭압과 같은 물리적인 현상이 복합적으로 발생되는 폭발 문제에서는 하나의 기법을 통해 복잡한 물리적인 현상을 수치적으로 재현하기는 사실상 불가능하며, 각 현상을 효과적으로 모사할 수 있는 수치해석 기법을 적합하게 선택하는 것이 중요하다. 폭풍파와 같은 대기전파의 경우 Euler 기법을 적용하면 매우 효과적으로 유체흐름을 모사할 수 있고, 반면에 폭풍파의 충격압력을 받는 콘크리트 구조물의 거동은 일반적으로 재료의 거동표현에 적합한 Lagrange 기법을 적용할 경우 효과적으로 모사될 수 있다. 본 연구에서는 Euler 기법과 Lagrange 기법을 커플링하는 복합적 수치해석 기법을 적용하여 비선형 동적해석이 수행된다. Euler-Lagrange 커플링 알고리즘은 Euler 영역과 Lagrange 영역이 겹치는 부분에서의 유체-고체의 상호작용이 고려되어야 한다. 상호작용 발생에 따라 규칙적이고 범위가 제한된 물리적 표면을 Lagrange 물체에 미리 정의해야 하며, 상호작용이 발생되는 접촉면에서 Lagrange 영역은 Euler 영역에 기하학적 구속조건으로 작용되고, Euler 영역은 Lagrange 영역에 압력경계조건으로 작용하게 된다. 이러한 표면은 요소절점들의 연결을 통해서 만들어지는 상호작용 다각형을 통해서 구현되어지며, 이는 Fig. 2에 도시적으로 나타나 있다. 이 때 작은 Euler 메쉬를 모델링할 경우 경계면에서의 정확한 재료 흐름을 모사할 수 있지만, 대량의 추가계산이 요구되기 때문에, 효율적인 Euler- Lagrange 커플링을 위해 적절한 메쉬 크기가 요구된다(^{Lu and Wang, 2006}).

3. RC building 구조물 폭발 해석

3.1 Euler-FCT Method

폭발 지점과 구조물과의 거리가 일정하게 떨어져 있는 경우 효율적인 폭발하중 정의를 위해 Euler-FCT (flux corrected transport) 기법(^{Zalesak, 1979})이 적용된다. 이는 대기 조건에서 사전 기폭시켜 필요한 정보를 미리 확보하여, 구조체 모델링이 완료된 후 확보된 정보를 입력하여 해석시간을 효과적으로 단축시킬 수 있다. 구체적인 기폭방법은 TNT 질량에 해당하는 구면체의 반지름(r)과 폭발지점에서 구조물까지의 거리(R)를 1차원 상에 모델링한 후, 사전 폭발시켜 R 지점에서의 밀도, 에너지, 속도, 그리고 압력 데이터가 확보되며, 이후 2D 및 3D 도메인에 맵핑(mapping)되어 최종적으로 Fig. 3(c)와 같이 폭발하중이 입력된다.

3.2 상세 제원

본 연구에서는 TNT 질량 20 kg의 폭발력을 받는 2-span, 3-bay를 가지는 3층 높이의 RC building 구조물에 대하여, 총 7 지점에서의 폭발 해석이 수행되고, 이에 대한 구조물 폭, 길이, 높이, 그리고 단면상세가 Fig. 4에 나타나 있다. 이때, RC building 구조물의 재료적 특성은 콘크리트 압축강도 24 MPa, 철근 항복강도 335 MPa, 주철근 지름 24 mm, 띠철근 10@200 mm이며, 슬래브 두께는 150 mm이다.

3.3 폭발하중, 기하학적 모델링 및 경계조건

Euler-FCT 기법을 이용한 시간에 따른 폭발압력을 계산하기 위해, TNT 질량 20 kg 일 경우 1차원에서 반지름(r) 143 mm, 100 kg 일 경우 반지름(r) 245 mm로 산정하여 기폭시켰으며, 폭발지점에서 구조물까지 거리(R)인 300 mm 떨어진 지점에서의 시간에 따른 폭발압력이 Fig. 5에 각각 나타나 있다. RC building 구조물의 폭발해석을 위한 3D 모델은 콘크리트의 경우 Solid 164 element, 철근의 경우 Beam 161 element가 각각 적용되며, 총 절점 및 요소 개수는 각각 603,316개 및 443,268개로 구성된다. TNT 폭발 위치는 외부 1 지점 및 내부 6 지점인 총 7 지점으로, 이격거리는 각 부재에서 300 mm씩 떨어진 지점에 위치하고 있으며, 이 때 TNT 질량은 20 kg 및 100 kg 일 때 각각 해석이 수행된다. 경계조건은 외부 대기를 형성하는 Euler 도메인의 지면을 제외한 면에 대해 Flow-out으로 설정하고, RC building 구조물의 기둥 아래 면은 Fixed support가 적용되며, 철근과 콘크리트의 접촉은 완전부착, 외부 대기와 구조물은 완전 커플링된다. 또한, RC Building 구조물의 주부재인 슬래브의 시간에 따른 처짐 측정 위치는 1-18번, 기둥의 경우 25-36번, 거더의 경우 37-54번, 그리고 부부재인 벽체의 측정 위치는 19-24번 위치로 설정되었으며, 각 지점별 시간에 따른 처짐값을 구하여 구조물의 거동 및 손상을 파악하고자 한다. 해석을 위한 pre-processing, 폭발지점, 측정위치, 그리고 경계조건이 Fig. 6 및 Table 2에 상세히 나타나 있다.

3.4 재료구성모델

본 연구에서는 탄성-점소성-손상이 연계된 RHT 콘크리트 모델이 적용된다(^{Riedel et al., 1999}). RHT 모델은 일반적인 취성재료를 분석하기 위한 모델로서 압력경화, 변형률 경화, 변형률 속도 경화, 누적손상, 변형률 연화 등과 같은 여러 특징을 도입할 수 있는 장점이 있다. 또한, 이 모델에서 Eqs. (10), (11) and (13)로서 제시되는 탄성한계면( $Y_{elastic}$), 재료파괴면( $Y_{fail}$), 파손된 재료의 잔류강도면( $Y_{fric}$)과 같은 세 강도 면을 도입하여 충격과 충돌을 받는 콘크리트재료의 비선형적 손상거동을 효과적으로 묘사할 수 있다.

Eq. (10)에서 $F_{elastic}$은 파괴면 강도에 대한 탄성강도의 비, $F_{ca}\left(p\right)$은 정수압력 상태의 탄성편차응력을 제한하는 함수이다. Eq. (11)에서 $F_{rate}\left(\dot{ϵ}\right)$은 변형률속도함수, $R_3\left(\theta \right)$는 날줄비(meridian ratio) $Q_2$와 제 2, 3 불변량의 함수로서 모델의 세 번째 불변량의존성을 정의하며 $Y_t$는 Eq. (12)에 나타나있다.

여기서, $f_c$는 압축강도, $A$는 파괴면 상수, $N$은 파괴면 지수, $p$는 $f_c$로 정규화한 압력, $p_{spall}^{}=p^{}(f_t/f_c)$, $\theta$는 콘크리트입자 사이의 각도, $\dot{ϵ}$는 변형률속도이다. Eq. (13)에서 $B$는 잔류파괴면 상수, $M$은 잔류파괴면 지수이다.

폭발하중상태의 보강 철근 비선형거동을 정확히 묘사하기 위해서는 철근의 변형률경화, 변형률속도경화, 열적 연화효과가 고려되어야 한다. 본 연구에서는 Eq. (14)과 같이 재료의 항복응력( $Y$)을 속도 의존적, 탄성적인 거동으로 잘 묘사할 수 있는 Johnson-Cook 모델(^{Johnson and Cook, 1983})을 채택하여 콘크리트에 보강된 철근을 모델링 하였다.

여기서, $Y_0$는 초기항복강도, ${ϵ}_{}$와 ${ϵ}^{*}$는 각각 소성 변형률과 정규화된 유효변형률속도를 의미하고, $C$, $D$, $n$, $m$은 재료상수, $T_h$는 $T_h=(T-T_{ref})/(T_{melt}-T_{ref})$로서 구해지는 상응(homologous) 온도이며, 여기서 $T_{ref}$와 $T_{melt}$는 각각 기준 및 융점온도를 나타낸다.

4. 폭발 해석결과 및 손상평가

TNT 질량 20 kg 및 100 kg에 대한 RC building 구조물의 폭발 해석 결과, 20 kg의 경우 폭발 지점이 각각 4, 6, 그리고 7지점에서, 100 kg 일 경우 3-7지점에서 주부재인 슬래브 및 보 부재에서 손상이 발생되었고, 나머지 폭발 지점의 경우, 안전성에 영향을 미치는 않는 수준의 손상이 발생되었으며, 손상이 나타는 폭발 지점에서의 최대 처짐 및 변위비가 Figs. 7 and 8에 상세히 나타나 있다. 또한, ^ASCE(1999)에서 제시하는 RC 구조물의 부재별 변위비( $\delta /L$)에 따른 파괴기준이 Table 3에 나타나 있고, TNT 질량 20 kg 및 100 kg 일 때의 파괴기준에 따른 손상 정도가 Tables 4 and 5에 상세히 나타나 있다.

먼저 TNT 질량 20 kg 일 경우 폭발 지점 4, 6, 그리고 7 지점에서 슬래브, 벽체, 기둥, 그리고 보 부재 중 슬래브에서만 주요 손상이 발생되었다. 각 지점별 슬래브 최대 처짐은 각각 283.6 mm, 233.0 mm, 그리고 349.1 mm가 발생되었고, 이때 각각의 변위비는 4.73 %, 3.88 %, 그리고 5.82 %로 이는 중간 손상, 가벼운 손상, 그리고 중간 손상에 해당된다. TNT 질량 100 kg 일 경우는 폭발 지점 3, 4, 5, 6, 그리고 7 지점에서 슬래브 및 보 부재에서 주요 손상이 발생되었다. 이 중 최대 처짐은 3 지점 슬래브의 경우 398.5 mm, 5 및 6 지점 슬래브의 경우 각각 339.4 mm 및 473.1 mm, 보 부재의 경우 각각 333.9 mm 및 430.7 mm 로 발생되었고, 이때 각각의 변위비는 6.64 %, 5,66 %, 7.18 %, 5.57 %, 그리고 7.18 %로 전체 중간 손상에 해당된다. 또한, 폭발 지점 4 및 7 지점에서는 슬래브에서만 최대 처짐 각각 751.9 mm 및 892.7 mm가 발생되었고 이때 변위비는 각각 12.53 % 및 14.8 %로 심각한 손상에 해당된다. 최종적으로, TNT 질량 20 kg 및 100 kg에 대한 전체 폭발 지점 1-7 지점에서의 손상 분포가 Figs. 9 and 10에 나타나 있으며 damage 1.0일 경우 완전파괴를 의미한다.

5. 결 론

비예측 극한하중인 폭발하중을 받는 RC building 구조물의 폭발손상평가를 위한 수치해석이 수행되었다. 해석의 효율성 및 정확성을 높이기 위해 Euler-FCT 기법을 이용한 사전 폭발하중이 정의되고, Euler-Lagrange 커플링 기법을 이용하여 유체-구조 연성 해석이 수행되며, 이때 고변형률 속도가 고려된 콘크리트 및 강재 재료구성모델이 적용되었다. 해석 결과, TNT 질량 20 kg 일 경우 3곳의 폭발손상 지점 중, 4번, 6번, 그리고 7번 지점에서 주부재인 슬래브에서만 중간 및 가벼운 손상이 발생되었고, 100 kg 일 경우 5곳의 폭발손상 지점 중, 5번 및 6번 지점에서 슬래브 및 보 부재 동시에 중간 손상이 발생되었으며, 3번 지점의 경우 슬래브에서만 중간 손상이 발생되었다. 나머지 2곳인 4번 및 7번 지점에서는 슬래브에서 심각한 손상이 발생되었다. 이는 향후 빌딩 구조물 뿐 만 아니라, 주요 사회기반구조물에 대한 폭발손상 및 안전성능 평가에 활용가능하다.