1. 서 론

외부나 내부의 응력변화에 따른 흙의 거동을 예측하기 위하여 사용되는 수치해석 방법은 흙을 연속체로 간주하는 유한요소법(Finite element method, FEM)과 유한차분법(Finite difference method, FDM), 그리고 흙을 이산체로 간주하는 개별요소법(Discrete element method, DEM)이 대표적이다. 유한요소법과 유한차분법은 지반을 요소(element)나 격자(grid)로 나눈 뒤 탄소성모델을 비롯한 흙의 구성식을 이용하여 다양한 연구 분야뿐만 아니라 실무에서도 널리 사용되고 있다. 하지만 사면 붕괴, 산사태, 연약지반 융기, 액상화로 인한 측방유동 등과 같은 지반의 대변형과 파괴 후 흙의 흐름과 같은 해석에서 격자를 사용할 경우 격자가 꼬이거나 뒤틀리기 때문에 계산되는 변위에 한계가 있을 뿐만 아니라 변형으로 인한 계산 정확도가 현저히 떨어지게 된다. 한편 격자를 사용하지 않는 개별요소법은 지반구조물의 대변형 해석이 가능하지만, 유한요소법이나 유한차분법과 같이 흙의 구성관계를 적용하기가 어렵다. 특히 일부 입력변수는 흙의 물성치와는 관계없기 때문에 변수 결정이 애매모호하다는 단점이 있다. 따라서 유한요소법이나 유한차분법 그리고 개별요소법과 같은 전통적인 수치해석 방법을 이용하여 지반공학에서 일어나는 대변형 문제를 해석하는 데는 한계가 있다.

입자법(Particle method)은 연속체역학의 수치 시뮬레이션법으로 연속체의 거동을 유한개의 입자 운동으로 계산하는 방법이다. 대표적인 입자법에는 Lucy(1977)가 개발한 SPH(Smoothed Particle Hydrodynamics)법과 Koshizuka et al.(1995)이 개발한 MPS(Moving Particle Semi-implicit)법이 있다. 입자로 해석 영역을 구성하므로 공간이 새롭게 형성되는 대변형 문제 해석에 적합하여 유체의 분열 또는 합체와 같은 비선형 자유표면을 시뮬레이션하는 유체역학에서 오래전부터 연구되었다(Daly et al., 1965). 또한 입자법은 파도나 지진해일의 거동 해석, 유체와 선박과 같은 구조물과의 상호작용 등 유체의 흐름을 해석하는 유체역학과 고체의 파괴나 균열과 같은 연속체의 대변형 등 고체역학 분야에 널리 응용되고 있다. 하지만 이와 같은 장점을 가진 입자법을 유체와 고체의 성질을 동시에 가지고 있는 흙의 대변형 문제에 적용하기 위한 연구는 최근 일본의 몇몇 연구자에 의해 시작되었다(Bui et al., 2008; Maeda and Sakai, 2004; Naili et al., 2005). 이들 대부분의 연구자들은 Liu and Liu(2003)가 유체 및 고체역학에 적용하기 위하여 개발한 SPH법을 바탕으로 지반공학 분야에 응용하기 위하여 흙의 구성식을 추가하는 연구를 진행하고 있다. Maeda et al.(2006)은 흙과 물의 상호작용을 고려한 모래지반의 침투파괴 해석 연구를 실시하였다. 최근 Bui et al.(2008)은 유체역학을 위한 SPH법을 기초로 탄완전소성모델인 Drucker-Prager모델을 비롯한 흙의 대표적인 탄소성모델인 수정 Cam-Clay모델까지 적용하였으며, 흙과 물의 상호작용 관계까지 고려하였다. 이와 같이 Bui et al.(2008)이 개발한 지반공학용 SPH법은 사면붕괴, 토석류 해석의 가능성을 보여 주었다. 한편 MPS법은 선박공학, 토목공학, 원자력공학 등의 공학문제에서 비압축성 유체를 해석하기 위하여 SPH법을 변형한 방법으로 입자간의 상호작용을 모델링하는 방식에 차이가 있다.

본 연구에서는 토사의 대변형 해석을 위해 MPS법에 기초하여 Mohr-Coulomb 파괴기준을 포함한 3차원 유동해석법을 개발하였으며, 모래기둥 붕괴실험과 점토의 일축압축시험을 실시하여 본 해석법의 적용성을 검토하였다. 이와 같은 토사 대변형 해석법은 기존의 흙의 파괴로 인한 토사의 변형 정도나 유동 범위를 수치해석 방법으로 정확하게 예측하여 지반구조물 설계에 반영하거나 기존 구조물을 보강하여 지반의 대변형으로 인한 피해를 방지할 수 있다.

2. 입자법(Particle method): MPS법

입자법의 기본적인 과제는 격자를 사용하지 않고 어떻게 미분방정식을 이산화할 것인가이다. 본 연구에서는 비압축성 유체의 흐름 예측을 위하여 유체공학 분야에서 개발된 수치해석법인 MPS법을 사용하였다. MPS법은 구배(Gradient), 발산(Divergence), 라플라시안(Laplacian)과 같은 미분연산자에 대응하는 입자간 상호작용모델을 이용해서 연속체의 지배방정식을 이산화하였다. 각 입자가 고유의 물리량을 가지고 임계거리 내에서 있는 다른 입자와 물리량을 주고받는다. 이 때 주고받는 물리량의 크기는 Fig. 1에서 나타나듯이 입자간 거리에 따라 가중치를 갖는다. 입자법에 사용되는 지배방정식, 주요 연산자 및 응력 계산방법은 다음과 같다.

2.1 지배방정식

비압축성 점성유동에 대한 지배방정식은 연속방정식(식 (1))과 Navier-Stokes 방정식(식 (2))이다. Navier-Stokes 방정식은 속도변화를 구하기 위한 식으로 압력항, 점성항, 그리고 외력항의 합으로 구성된다.

 (1)

 (2)

여기서, ρ는 밀도, t는 시간, 는 속도벡터, 는 구배연산자, P는 압력, 는 동점성계수, 는 외력을 나타낸다.

2.2 입자간 상호작용 모델

2.2.1 Kernel function

Kernel function은 식 (3)과 같이 입자간의 거리 에 따라 입자간 상호작용에 가중치 을 두기 위한 함수이다. 즉, 입자가 가까울수록 가중치가 커지고 반대로 두 입자의 거리가 충분히 멀어 임계거리 이상이 되면 가중치를 무시할 수 있으므로 0으로 본다.

 (3)

식 (3)은 Koshizuka and Oka(1996)가 제안한 Kernel function으로 두 입자 사이의 거리가 매우 가까울 때 상호작용 가중치 값이 무한대로 증가하는 특징이 있다. 본 연구에서는 유동이 크게 격렬하지 않은 시뮬레이션에서 보다 안정적인 압력을 얻기 위해 식 (3)을 대신하여 Jeong(2008)이 제안한 식 (4)와 같은 Kernel function을 사용하였다. 이 식을 이용하면 입자간 거리가 가까워 질 때 가중치가 특정 값으로 수렴한다.

 (4)

여기서, 임계거리 의 값은 입자가 평균거리에 임의의 값 2.1, 3.1 등의 값을 곱해 계산한다.

2.2.2 구배모델(Gradient model)

입자 와 가 각각 물리량 와 를 가지고 있다고 할 때 Fig. 2와 같이 입자 와 주변입자 사이의 구배벡터를 로 계산한다. 여기에 Kernel function으로 가중치를 곱해주고 입자가 여럿일 경우 전체의 합을 초기 배치에서 입자수 밀도인 로 나누어 식 (5)와 같이 구배를 계산한다.

 (5)

여기서, d는 계산 공간의 차원이다. 의 입자수 밀도를 나타내는 는 다음과 같이 정의한다.

 (6)

비압축성 조건에서는 밀도가 일정하므로 입자수 밀도도 일정해야 한다. 가 그 기준이 되는 입자수 밀도로 입자의 초기 배치상태에서 결정되며 계산 과정에서 변하지 않는 고정 값이다. 는 입자 의 임계거리 안에서 입자의 배치상태를 나타내는 값으로 임계거리 내에 입자가 많을수록, 그리고 임계거리 내에서도 중심입자 가까이에 입자가 많이 분포해 있을수록 값이 커진다. 값은 전체 입자 배치의 대푯값으로 이와 비교하여 값이 크다면 입자가 빽빽한 것으로 작다면 느슨한 상태인 것으로 판단할 수 있다. 따라서 는 전체 입자배치를 잘 대표할 수 있는 값으로 정해야 하며, 본 연구에서는 시료의 중심에 있는 입자를 의 기준입자로 사용하였다.

2.2.3 라플라시안 모델(Laplacian model)

라플라시안()은 식 (7)과 같이 물리적으로 확산을 의미하므로 Fig. 3에서 나타나듯이 입자법에서 라플라시안 모델은 입자 i가 갖는 물리량의 일부를 주위 입자 j로 분배해주는 역할을 한다.

 (7)

여기서, 는 식 (8)과 같이 해석해와 통계적인 확산의 정도를 일치하게 하기 위한 계수이다.

 (8)

는 시간과 입자의 위치에 따라 값의 변화가 매우 작으므로 고정된 값을 사용할 수 있다.

2.2.4 비압축성 모델

MPS법의 비압축성 유체의 움직임은 두 단계로 나눈

다. 1단계에서는 점성항과 외력항(중력 등)의 계산을 통

해 입자의 중간속도()와 중간 위치()를 식 (9)와 같이 계산한다.

 (9)

2단계에서는 먼저 압력의 포아송 방정식을 식 (10)과 같이 계산한다.

 (10)

는 1단계에서 계산한 중간위치를 이용하여 구한 입자수 밀도이다. 식 (10)의 우변은 입자수 밀도의 차이로 압력의 확산()을 나타내고 있다. 좌변은 라플라시안 모델인 식 (7)로 나타낼 수 있어 연립 1차방정식으로 압력을 구할 수 있다. 구한 압력을 식 (10)에 대입하여 속도 수정치 를 다음과 같이 계산한다.

 (11)

는 식 (5)를 사용하여 구할 수 있다.

최종적으로 시각 n+1에서의 입자의 속도와 위치는 다음과 같이 구해진다.

 (12)

 (13)

입자법의 전체적인 응력 계산 흐름도는 Fig. 4와 같다.

2.3 응력 계산

2.3.1 수직 및 수평 응력

입자 i+1의 수직응력 을 계산하기 위해 입자 i+1은 식 (14a)와 같이 자신의 바로 위 입자 i의 응력 에만 영향을 받는다고 가정한다. 모든 입자는 자신의 무게만큼 바로 아래 입자에 응력을 전달한다. 따라서 표면에 있는 입자의 응력은 0이다. 그리고 모든 입자는 자신에게 영향을 주는 입자의 응력을 전달받는다.

 (14a)

 (14b)

여기서, 는 응력, 는 입자의 무게, A는 단면적을 나타낸다. 입자의 무게는 입자의 모양을 입자간 평균거리(dis)를 지름 혹은 한 변으로 가지는 구형 또는 정사각형으로 가정하여 구할 수 있다. 식 (14b)는 흙의 자중 이외에 외력이 작용할 때 응력을 구하는 식으로 F는 변형속도에 따른 힘을 나타내며, A는 힘을 받는 면적을 나타낸다. 입자의 수평응력은 입자의 수직응력에 횡토압계수 k를 곱하여 구할 수 있다.

2.3.2 Mohr-Coulomb 파괴기준

흙 입자의 파괴 여부는 Fig. 5와 같은 Mohr’s Circle을 기초로 하여 계산하였다. 마찰각 와 점착력 가 주어졌을 때 응력 에서 파괴면은 다음과 같다.

 (15)

와  이 주어졌을 때 최대 전단응력 는 다음과 같이 구할 수 있고

 (16)

이 때 파괴면에서 전단응력 는 다음과 같이 구할 수 있다.

 (17)

과 에 대하여 전단파괴가 일어나는 파괴면에서의 응력 는

 (18)

이고, 이 응력 에서 전단파괴가 일어날 때 전단응력 는 다음과 같다.

 (19)

즉, 파괴면에서의 전단응력 와 파괴가 일어날 때의 전단응력 가 일치할 때 전단파괴가 일어나는 것으로 본다.

2.3.3 항복(Yielding) 및 경화법칙(Hardening rule)

일반적으로 흙은 Fig. 6과 같이 전단응력이 Mohr- Coulomb과 같은 파괴기준에 도달하기 전에도 변형과 동시에 항복으로 인한 소성변형이 발생한다(Park et al., 2005). 한편 본 연구에서는 기존의 연속체역학에서 사용되는 다소 복잡한 항복 및 경화법칙(Park et al., 2005)을 입자법에 적용하기 보다는 흙의 거동을 점성 유체로 가정한 다음 이와 같은 항복으로 인한 경화현상을 시뮬레이션하고자 하였다. 본 연구에서 흙의 파괴 이전에 발생하는 소성변형을 시뮬레이션하기 위해 사용한 동점성계수는 0.13와 같으며, 흙이 파괴기준을 만족할 경우에 흙의 점성이 소실된다.

3. 모래기둥 붕괴실험 및 입자법 시뮬레이션 비교

건조한 모래에 구속압이 작용하지 않을 경우 자중에 의해 흘러 내리는 대변형 사례를 재현하였으며, 이를 3차원 입자법을 이용하여 시뮬레이션하였다. 3차원으로 흐르는 흙의 유동을 시간에 따라 비교하기 위해 Mohr-Coulomb 파괴기준을 결정하는 강도정수인 점착력과 마찰각은 모두 0으로 설정하여 모래기둥이 강도가 없는 유체처럼 유동하는 과정을 3차원적으로 시뮬레이션하였다.

3.1 모래기둥 붕괴실험

모래기둥 붕괴실험에서는 직경 28cm, 높이 50cm인 아크릴 원기둥 내에 낙동강 모래를 넣고 Fig. 7과 같이 빠른 속도로 올렸을 때 원기둥 내 모래가 흘러내리도록 실험하였다. Fig. 7은 0.0초, 0.6초, 1.2초에 상부(top view)와 정면(side view)에서 본 모래기둥의 흐르는 모습을 촬영한 사진이다. 점성이 없으면서 구속압이 작용하지 않을 경우 어떠한 강도도 가지지 않는 모래기둥이 빠르게 흘러내리면서 1.2초에 폭 120cm, 높이(중앙 부분) 15cm로 변하였다.

3.2 모래기둥 붕괴실험 시뮬레이션

시간에 따른 모래기둥 붕괴현상을 입자법으로 시뮬레이션하였다. Fig. 8은 모래기둥이 붕괴되는 과정을 위에서 내려다 볼 경우(top view)와 앞에서 볼 경우(side view)에 모래기둥의 3차원 거동을 입자법이 시뮬레이션하고 있다. 시간이 지남에 따라 모래가 점점 흘러내리면서 주변으로 확대되어 가는 것을 볼 수 있다. 기존 유체해석에 사용되는 입자법 프로그램은 대부분 2차원 해석만 가능하지만, 본 연구에서 개발 중인 프로그램은 토사의 대변형이나 유동을 3차원적으로 시뮬레이션하므로 좀 더 정확한 경계조건을 수치해석에 반영할 수 있다. Fig. 8은 Fig. 7의 실험 결과와 비교적 잘 일치하는 것으로 나타났다. 본 해석에서는 모래기둥의 강도를 0으로 보았기 때문에 응력계산 보다는 3차원 변형 형태만을 비교하고자 하였다.

4. 점토기둥 붕괴실험 및 입자법 시뮬레이션 비교

점착력으로 자립이 가능한 점토를 이용하여 실내 붕괴실험(일축압축시험)을 실시하였다. 원기둥 모양의 점토를 성형한 다음 하부에서 하중을 가하여 붕괴시키는 실험을 실시하였으며, 이를 본 연구에서 개발한 입자법을 이용하여 파괴 이전의 소성변형을 시뮬레이션하고자 하였다. 또한 점토 내부에 작용하는 응력상태가 Mohr-Coulomb 파괴기준을 만족할 경우에는 비점성 유체와 같은 대변형이 발생하게 된다. 수치해석에 필요한 입력변수는 점토 비중 2.65, 횡방향토압계수 0.5, 압축 속도 그리고 앞서 언급한 Mohr-Coulomb 파괴기준의 강도정수이다. 본 연구에 사용한 점토를 stiff clay로 판단하여 비배수 전단강도는 Das(2010)를 참고하여 점착력 100kPa과 마찰각 0도를 사용하였다. 그 밖에도 계산 과정에 필요한 입력변수로는 의 값 혹은 값 계산에 사용할 입자의 번호가 필요한데 값이 전반적인 수치해석 결과에 큰 영향을 미치므로 주의해서 선택해야 하며, 본 연구에서는 값으로 4.879를 사용하였다.

4.1 점토기둥 붕괴실험

본 연구에 사용한 점토는 진해 인근 해역에서 채취한 해성점토로 직경 5cm, 높이 10cm의 원기둥 모양의 점토를 성형한 다음 하부에서 느린 속도(10mm/min)로 하중을 가하여 붕괴실험을 실시하였다. Fig. 9는 실험 시작 후 30, 90, 180, 360초 일 때 점토기둥의 변형 모습을 보여주고 있으며, 360초 경과하였을 때는 높이가 10cm에서 4cm까지 압축되었다.

4.2 점토기둥 붕괴실험 시뮬레이션

Fig. 10은 시간에 따른 점토기둥의 붕괴과정을 시뮬레이션하고 있으며, Fig. 9의 실험결과와 유사하게 변형되는 것을 알 수 있다. 수직응력은 어떤 지점에 작용하는 흙의 무게와 축변형에 따라 입자에 작용하는 힘(F=ma)을 합하여 식 14(b)와 같이 계산하였으며, 전단응력은 식 (17)에 따라 파괴면에서의 전단응력으로 계산하였다. Fig. 11은 축변형에 따른 점토기둥 표면에 작용하는 전단응력 변화를 나타내고 있다.

Fig. 12는 점토기둥이 변형됨에 따라 중간(Middle) 지점의 내부(center)와 외부(surface)를 비롯하여 상부(Top)와 하부(Bottom)의 내부에 작용하는 수직응력 및 전단응력의 변화를 비교하고 있다. 여기서 내부는 단면의 중앙지점에 있는 입자의 응력을, 외부는 가장 바깥쪽에 있는 입자의 응력을 나타낸다. Fig. 12(a)는 중간 지점의 내부와 외부 수직응력(Axial stress)과 전단응력(Shear stress)을 비교하고 있으며, 내부에 작용하는 응력은 60mm 변형까지 지속적으로 증가하는 경향을 보이지만, 외부에 작용하는 응력은 50mm 변형에서부터 서서히 감소하는 경향을 보였다. 이는 흙 입자가 밖으로 밀려나면서 외력에 의한 영향이 줄어들었기 때문으로 보인다. 그리고 하부에 작용하는 응력이 상부에 작용하는 응력보다 다소 큰 것을 알 수 있다. 점토기둥의 아래 부분의 응력이 점차 증가하면서 Mohr-Coulomb 파괴기준을 만족하지는 않지만 항복으로 인한 지속적인 소성변형이 발생하였다.

본 연구에서 개발한 입자법에 Mohr-Coulomb 파괴기준을 적용하였으며, 또한 흙을 점성 유체로 가정함으로써 점토기둥 내 응력 증가에 따른 항복이나 경화현상에 따른 대변형을 시뮬레이션할 수 있었다. 한편 PFC와 같은 개별요소법도 이와 유사한 대변형 해석이 가능하지만 Mohr-Coulomb 파괴기준이나 항복 및 경화법칙과 같은 소성이론을 적용하기에는 어려운 점이 있다.

5. 결 론

본 연구에서는 유체의 다양한 경계치 문제에 적용되는 입자법을 함수비 증가에 따라 고체에서 유체로 거동하는 토사의 대변형 해석에 적용하고자 하였다. 입자법은 대변형 해석이 가능하면서도 흙의 구성식 도입이 가능하다는 장점이 있다. 따라서 본 연구에서는 Mohr-Coulomb과 같은 파괴기준을 입자법에 수식화하였으며, 이를 이용하여 모래와 점토의 붕괴실험 결과를 시뮬레이션하여 개발된 입자법의 대변형 해석에 대한 적용성을 검토하였다. 그 결과는 다음과 같다.

(1)흙에 작용하는 수직응력은 입자간 속도 변화와 자중을 이용하여 계산하였으며, 수평방향으로 작용하는 응력은 횡방향 토압계수를 이용하여 계산하였다.

(2)Mohr-Coulomb의 파괴기준을 입자법에 도입한 다음 계산된 흙의 응력상태가 Mohr-Coulomb의 파괴기준을 만족할 경우 파괴로 인한 대변형이 발생하도록 하였다. 한편 파괴기준을 만족하기 이전에는 입자법의 점성항을 이용하여 항복으로 인한 흙의 소성거동을 시뮬레이션하였다.

(3)개발된 입자법은 구속압이 작용하지 않고 강도가 0인 모래기둥의 3차원 흐름이나 점착력을 가진 점토기둥의 일축압축으로 인한 3차원 대변형과 응력 변화를 계산할 수 있었다. 입자법 해석 결과는 실험 결과와 비교적 일치하였다.

(4)본 연구에서 개발한 입자법은 하중 증가에 따른 응력 계산이 가능할 뿐 아니라 이로 인한 대변형 시뮬레이션이 가능하였다. 따라서 개발된 입자법은 사면붕괴와 같은 다양한 지반공학 대변형 예측에도 활용 가능할 것으로 판단된다.