Immersed Boundary Method - 1. 간단한 힘과 에너지의 관계

 

 

나는 유체-구조 상호작용[Fluid-Structure Interaction, FSI] 시뮬레이션을 주로 하고 있는데, 이는 간단하게 말하면 유체와 고체를 동시에 수치해석하는 것이다. 이를 위해서는 당연히, 유체 격자와 고체 격자가 각각 필요할 것이고, 유체 격자와 고체 격자 사이에서 유속, 변위, 힘 등의 변수를 주고받게 할 수 있는 수치적 전략이 필요하다. 

 

이러한 전략 중에 가장 쉽게 떠올릴 수 있는 것은 흐르는 유체와 변형하는 고체[여기서는 특히 탄성체] 사이의 경계가 정확히 보장되도록 고체의 움직임에 맞게 유체 격자도 움직여주는 것이다. 이 방법을 Arbitrary Lagrangian Eulerian [ALE]라고 하고, 아래의 그림은 ALE로 시뮬레이션한 결과의 일부이다. 빨간색이 고체 도메인이고 이를 둘러싼 회색이 유체 도메인인데, 고체의 휨에 맞춰서 유체 격자도 움직이는 것을 볼 수 있다. 

 

 

이 방법의 가장 큰 장점은 일단 시뮬레이션을 위해 필요한 수치해석적 지식이 일반적인 유체만 수치해석하는 CFD에 비해서 그렇게 많이 필요하지 않고, 사용하기 쉽다는 점이 있다. 다만 단점은, 유체 격자의 품질이 크게 저하될 수 있다는 점인데, 격자 품질 저하는 오차를 누적시키거나 수치적 안정성을 떨어뜨린다. 따라서 고체가 크게 변형하는 경우나, 많은 고체가 변형하는 경우에서는 아무리 수치적으로 안정성을 높이는 스킴을 사용해도 시뮬레이션이 터지는 상황이 자주 발생한다. 유체 격자를 더 fine하게 설정해서 한 cell에서의 오차 누적을 줄이는 방향으로 해결할 수도 있지만, 유체 격자를 fine하게 했을 때 수치적 안정성을 보장하기 위해서는 timestep을 더 낮춰야 해서 원하는 양의 자료를 얻기까지 매우 많은 시간이 걸린다. 예를 들어서, 격자 크기와 timestep과 유속과 관련되어서 수치적 안정성에 대한 기준이 될 수 있는 Courant-Friedrichs-Lewy [CFL] number라는 것이 있다. 만약 격자를 2배 더 fine하게 하려면 3차원 시뮬레이션의 경우에서는 x, y, z 방향 모두 2배가 되므로 격자의 수가 8배가 되고 이 때 같은 수준의 CFL number를 유지하기 위해서는 timestep이 2배 이상 더 작아져야 한다. 따라서 2배 fine하게 했을 때 시뮬레이션 시간은 최소 16배 이상 더 걸리는 것이다. 물론 이는 컴퓨터가 좋으면 해결되는 문제이긴 하다만, 만약 현재 fine한 정도로도 내가 원하는 유동 현상을 볼 수 있는데 오직 ALE의 문제를 해결하기 위해서 격자를 fine하게 한다면, 이는 사용 가능한 자원을 다소 낭비하는 셈이 될 수 있다. 

 

이러한 단점을 극복할 수 있는 것 중 한 가지 전략이 이 글의 제목이기도 한 Immersed Boundary Method [IBM]이다. IBM은 1972년에 Charles Peskin이라는 아저씨가 처음에 혈류역학 시뮬레이션을 위해서 고안한 기법이고, 그 이후에 많은 발전이 있어왔다. Immersed boundary라는 이름에서 알 수 있듯이, 이 방법의 핵심 아이디어는 고체의 경계가 유체의 cell face와 정확히 일치하지 않고, 어떤 arbitrary한 경계가 유체 도메인 내부에 잠겨있을 때, 유체 격자를 움직이지 않고 그 주변 유체 cell의 정보를 사용해서 고체의 경계가 거기에 있다고 생각하는 것이다. IBM을 조금이라도 아는 사람이면 이렇게만 말해도 대충 이해할 수는 있겠지만, 전혀 모르면 이게 뭔 소리인지 전혀 감이 오지 않을 것이다. 나도 처음에 그러했으니까. 먼저 IBM에 대한 이론을 처음부터 살펴보겠다. 

 

먼저 IBM의 equation of motion을 유도하는 과정에서 기본이 되는 힘과 에너지의 관계에 대해서 살펴본다. 

 

3차원 공간에서, 유체는 움직이지 않는 Eulerian grid로 설명할 수 있고 고체는 움직이는 Lagrangian grid로 설명할 수 있다. 여기서 고체를 설명하는 Lagrangian curvilinear coordinate를 $(q,r,s)$라고 하자. 이 $(q,r,s)$는 실제 물리적 위치를 나타내는 좌표라기보다도 그 고체의 material point에 해당하는 하나의 label이라고 생각하면 된다. 그리고 $\mathbf{X} (q,r,s,t)$를 시간 $t$에서 그 label에 해당하는 고체 point의 Eulerian grid에서의 물리적 위치이다. 이러한 coordiante system에서 먼저 고체의 탄성 에너지 $E[\mathbf{X}]=E[\mathbf{X}(q,r,s,t)]$라고 생각할 수 있다. 

 

이 때 고체가 아주 조금 움직인 상황을 생각하자. 이 경우는 수학적으로 위치에 아주 작은 perturbation이 있다고 표현할 수 있는데, 이러한 small perturbation은 주로 $\delta$ 심볼로 표현하지만 IBM에서는 Dirac delta가 매우 중요한 역할을 하기 때문에 $\delta$는 이를 위해서 사용하고 small perturbation은 Weierstrass p인 $\wp$를 사용해서 나타내겠다. 이때 위치에 small perturbation $\wp \mathbf{X}(,,,t)$ 에 의해서 생기는 Energy의 small perturbation $\wp E[\mathbf{X}(,,,t)]$을 생각할 건데, 1차까지는 small perturbation $\wp E$가 $\wp X$에 대해서 linear하다고 표현될 수 있다. 다음과 같이 말이다.

$$E[X+\epsilon \wp \mathbf{X}] = E[\mathbf{X}] + \epsilon \wp E[X] + O(\epsilon ^2 \wp ^2 )+ \cdots $$

즉, $\wp E[\mathbf{X}(,,,t)]$가 linear function이고, linear function은 어떤 두 벡터의 내적으로 쓸 수 있다. 이거는 또 무슨 말이냐면 다음과 같은 선형 함수 $L(v)$를 생각하자.

$$L(v) = a_1 v_1 + a_2 v_2 + \cdots + a_n v_n $$

그런데 이건 정확히

$$L(v) = a \cdot v $$

라고 쓰는 것과 동일하다는 소리이다. 

따라서 energy의 small perturbation $\wp E$는 어떤 $-\mathbf{F}$라는 벡터에 대해서 다음과 같이 쓸 수 있다. 

$$\wp E[\mathbf{X}(,,,t)] = \int (-\mathbf{F}(q,r,s,t)) \cdot (\wp \mathbf{X} (q,r,s,t)) dq dr ds $$

원래 내적의 정의에서는 $\Sigma$가 여기서는 연속 공간에서의 변수를 다루므로 적분으로 바뀐 것이다.

그리고 여기서 $\mathbf{F}$는 internal force density로 새각하면 되는데, 다음을 보면 더 clear해진다.

$$\mathbf{F} = -\frac{\wp E}{\wp \mathbf{X}}$$

에너지의 small perturbation에 대한 변분이니까 force density라는 것은 clear하고, 부호가 마이너스인 이유는, 탄성체에서 internal force는 에너지를 감소시키는 방향으로 작동하기 때문이다. 고무줄과 같이 변위가 커졌을 때 포텐셜 에너지가 커지면 여기서 작용하는 탄성력은 원래대로 변위를 줄이는 방향으로 작용한다는 것이라고 이해하면 쉽다. 

 

논문에서는 여기까지 설명하고, 잠시 다른 얘기를 한다. 탄성 섬유의 에너지와 이 에너지에 의해서 만들어지는 탄성력의 관계를 보여주는 예시를 설명하는데, 사실 크게 중요하지는 않지만 수학적으로 꽤나 간단해서 나도 잠깐 이 내용을 적는다. 

 

여러 탄성 섬유가 있는 시스템을 생각하고, 이 때 탄성 섬유의 에너지를 $E$라고 하고 탄성력을 $\mathbf{F}$라고 하고, 위치에 따라서 부드럽게 변하는 섬유의 방향을 $\mathbf{\tau}$라고 하자. 그리고 $q, r, s$를 각각 섬유에 대응하는 라그랑지안 curvilinear 좌표라고 하고, 이 때 $q, r$은 각 섬유 줄기 하나에 대해서 constant하자고 하자. 즉, 한 섬유 내에서는 $s$만 변하는 것이다. 그렇다면 섬유가 갖는 에너지는 섬유의 strain의 크기에 대한 함수이고, 이 arbitrary한 에너지 밀도 함수를 $\varepsilon$이라고 했을 때 다음과 같이 쓸 수 있다. 

$$E=\int \varepsilon \left( \left| \frac{\partial \mathbf{X}}{\partial s} \right| \right) dq dr ds $$

이제 여기서 힘 $\mathbf{F}$를 유도하기 위해 에너지의 변분 $\wp E$를 생각할 것이다.

먼저 변분은

$$\mathbf{X} \to \mathbf{X}+\alpha \wp \mathbf{X}$$

에 대해서 

$$ \wp E = \left. \frac{d}{d\alpha} E[\mathbf{X} + \alpha \wp \mathbf{X}] \right|_{\alpha=0}$$

이다. 따라서 위에서 적분으로 정의된 에너지의 식 안에 미분을 넣으면,

$$\wp E = \int \left. \frac{d}{d\alpha} \varepsilon \left( \left| \frac{\partial}{\partial s} \left( \mathbf{X} + \alpha \wp \mathbf{X} \right) \right| \right) \right|_{\alpha =0} dqdrds $$

이 되는데, 여기서 chain rule에 의해

$$\frac{d}{d\alpha} \varepsilon \left( \lambda (\alpha) \right) = \varepsilon '(\lambda) \frac{d \lambda}{d \alpha} \quad \mathrm{where} \quad \lambda(\alpha) = \left| \frac{\partial}{\partial s} \left( \mathbf{X} + \alpha \wp \mathbf{X} \right) \right| $$

이 된다. 그리고 여기서 벡터의 크기에 대한 미분이 등장하는데, 어떤 벡터 $\mathbf{v} (\alpha)$의 norm을 $\alpha$로 미분하면

$$\frac{d}{d\alpha} \left| \mathbf{v} \right| = \frac{\mathbf{v}}{\left| \mathbf{v} \right|} \cdot \frac{d\mathbf{v}}{d\alpha}$$

이다. $\left| \mathbf{v}(\alpha) \right| = \sqrt{\mathbf{v}(\alpha) \cdot \mathbf{v}(\alpha)} $를 생각하면 쉽다.

따라서 위의 $\lambda$의 정의에서 절대값 기호 안에 있는 항을 $\mathbf{v}$라고 하면, 

$$\mathbf{v}(\alpha) = \frac{\partial}{\partial s} \left( \mathbf{X} + \alpha \wp \mathbf{X} \right)$$

$$\frac{d \mathbf{v}}{d \alpha} = \frac{\partial}{\partial s} \left( \wp \mathbf{X} \right)$$

$$ \mathbf{v}(0) = \frac{\partial \mathbf{X}}{\partial s} \quad \mathrm{at} \quad \alpha = 0$$

이므로 

$$\left. \frac{d\lambda}{d\alpha} \right|_{\alpha=0} = \frac{d}{d\alpha} \left. \left| \frac{\partial \mathbf{X}}{\partial s} + \alpha \frac{\partial (\wp \mathbf{X})} {\partial s} \right| \right|_{\alpha=0} = \frac{\partial \mathbf{X} / \partial s}{\left| \partial \mathbf{X} / \partial s \right| } \cdot \frac{\partial (\wp \mathbf{X} )}{\partial s} $$

따라서 

$$\wp E = \int \varepsilon ' \left( \left| \frac{\partial \mathbf{X}}{\partial s} \right| \right) \frac{\partial \mathbf{X} / \partial s}{ \left| \partial \mathbf{X} / \partial s \right| } \cdot \frac{\partial (\wp \mathbf{X} )}{\partial s} dq dr ds$$

이고, 이를 내적의 앞과 뒤로 나눠서 $s$에 대해서 부분적분하면,

$$\wp E=\left[ \varepsilon ' \left( \left| \frac{\partial \mathbf{X}}{\partial s} \right| \right) \frac{\partial \mathbf{X} / \partial s}{\left| \partial \mathbf{X} / \partial s \right| } \cdot \wp \mathbf{X} \right]_{s-boundary} - \int \frac{\partial}{\partial s} \left( \varepsilon ' \left( \left| \frac{\partial X}{\partial s} \right| \right) \frac{\partial \mathbf{X} / \partial s}{ \left| \partial \mathbf{X} / \partial s \right| } \right) \cdot \wp \mathbf{X} dq dr ds $$

그런데 여기서 대괄호 안의 항은 boundary condition인데, 특히 고정단에서는 $\wp \mathbf{X} =0$이고 free end에서는 이계도함수가 0이라서 결국 0이다. 따라서

$$\wp E = - \int \frac{\partial}{\partial s} \left( \varepsilon ' \left( \left| \frac{\partial \mathbf{X}}{\partial s} \right| \right) \frac{\partial \mathbf{X} / \partial s}{\left| \partial \mathbf{X} / \partial s \right| } \right) \cdot \wp \mathbf{X} dq dr ds$$

이다. 그리고 위에 있는 에너지의 변분으로 정의되는 힘에 의해서,

$$\mathbf{F} = \frac{\partial}{\partial s} \left( \varepsilon ' \left( \left| \frac{\partial \mathbf{X}}{\partial s} \right| \right) \frac{\partial \mathbf{X} / \partial s}{\left| \partial \mathbf{X} / \partial s \right| } \right)$$

이다. 여기서,

$$T=\varepsilon ' \left( \left| \frac{\partial \mathbf{X}}{\partial s} \right| \right) = \mathrm{fibre\ tension}$$

이라고 한다면

$$\mathbf{\tau} = \frac{\partial \mathbf{X} / \partial s}{\left| \partial \mathbf{X} / \partial s \right| } = \mathrm{unit\ tanget\ to\ the\ fibres.}$$

이고, 따라서

$$\mathbf{F} = \frac{\partial}{\partial s} \left(T \mathbf{\tau} \right)$$

이 된다. 이 결과는 탄성 섬유에 의해서 생성되는 탄성력은 섬유에 평행한 방향의 힘만 가지고, $\mathbf{\tau}$와 $\partial \mathbf{\tau} / \partial s$에 의해서 생기는 평면과 수직인 방향의 힘은 없다는 것을 의미한다. 

 

여기까지 아주 간단한 힘과 에너지의 관계에 대해서 정리했고, 다음의 글은 '이러한 힘이 작용할 때 실제로 탄성체가 어떻게 움직이는가?'에 대해서 정리할 것이다. '그냥 실제로 운동량의 balance대로 움직이는 것이 아닌가?'라고 생각할 수도 있는데, IBM을 유도하는 과정에서는 라그랑주 역학적인 접근이 필요하고, 라그랑주 역학의 핵심 아이디어는 운동량 보존이 아니라 다른 곳에 있다. 다음 글이 그에 대한 내용이 될 예정이다.