본문 바로가기

CFD

Immersed Boundary Method - 3. Equations of Motion

 

앞선 글에서 action $S$를 최소화하기 위해서 변분 $\wp S$를 유도하였고, 이 변분의 제약인 $dJ/dt=0$을 라그랑지안 좌표계에서 그대로 다루는 것이 어렵기때문에 라그랑지안 좌표 $\mathbf{X}$의 변분 $\wp \mathbf{X}$에 대응하는 pseudo-velocity $\mathbf{v}$를 도입하여서 incompressibility가 오일러리안 좌표계에서

$$\nabla \cdot \mathbf{v} = 0$$

와 같이 표현될 수 있음까지 유도하였다. 이제 라그랑지안 변수의 적분으로 표현된 $\wp \mathbf{S}$를 오일러리안 변수를 도입하여서 나타낼 준비가 된 것이다. 라그랑지안 변수를 오일러리안 변수로 표현하는데 Dirac 델타 함수를 사용하는데, 먼저 다음의 두 식은 쉽게 이해가 된다. 여기서 대문자 $\mathbf{X}$는 탄성체가 실제로 있는 위치를 나타내는 것이고 소문자 $\mathbf{x}$는 오일러리안 공간 변수를 나타내는 것이다.

$$\wp \mathbf{X}(q,r,s,t)=\int \mathbf{v}(\mathbf{x},t) \delta (\mathbf{x}-\mathbf{X}(q,r,s,t))d\mathbf{x}$$

$$\frac{\partial ^2 \mathbf{X}}{\partial t^2} (q,r,s,t) \cdot \wp \mathbf{X}(q,r,s,t) = \int \frac{D \mathbf{u}}{Dt} \cdot \mathbf{v}(\mathbf{x},t) \delta (\mathbf{x} - \mathbf{X}(q,r,s,t))d \mathbf{x}$$

그리고 action의 변분은 다음과 같이 표현되었었다.

$$0=-\wp S = \int_{0}^{T} \int \left( M \frac{\partial ^2 \mathbf{X}}{\partial t^2} - \mathbf{F} \right) \cdot \wp \mathbf{X} dq dr ds \ dt $$

이 action의 변분 식에 위의 두 식을 대입하면 다음을 얻는다.

$$0=\int_{0}^{T} \int \int \left( M(q,r,s) \frac{D \mathbf{u}}{Dt}(\mathbf{x},t) - \mathbf{F}(q,r,s,t) \right) \cdot \mathbf{v}(\mathbf{x},t) \delta (\mathbf{x} - \mathbf{X}(q,r,s,t)) d\mathbf{x} \ dq dr ds \ dt $$

여기서 $M(q,r,s)$와 $\mathbf{F}(q,r,s,t)$는 각각 mass density와 elastic force density로, 모두 라그랑지안 변수이다. 이 두 라그랑지안 변수만 없애면 $q, \ r, \ s$에 대한 적분도 없앨 수 있다. 이 두 변수도 디랙 델타 함수로 다음과 같이 오일러리안 변수로 바꿀 수 있다. 

$$\rho (\mathbf{x},t) = \int M(q,r,s) \delta (\mathbf{x} - \mathbf{X}(q,r,s,t))dq dr ds$$

$$\mathbf{f}(\mathbf{x},t) = \int \mathbf{F}(q,r,s,t) \delta(\mathbf{x} - \mathbf{X}(q,r,s,t)) dq dr ds$$

이 두 개는 딱봐도 오일러리안 질량밀도와 오일러리안 탄성력밀도임이 뻔하지만, 조금 더 명확하게 하기 위해서 양변을 임의의 영역 $\Omega$에 대해서 적분한다. 여기서 $\int_{\Omega} \delta (\mathbf{x} - \mathbf{X}(q,r,s,t))d\mathbf{x}$의 값은 $\mathbf{X}(q,r,s,t) \in \Omega$의 여부에 따라서 0 또는 1이기때문에 적분의 결과는 다음과 같다.

$$\int_{\Omega} \rho (\mathbf{x},t) d\mathbf{x}= \int_{\mathbf{X}(q,r,s,t) \in \Omega} M(q,r,s) dq dr ds $$

$$\int_{\Omega} \mathbf{f} (\mathbf{x},t) d\mathbf{x} = \int_{\mathbf{X}(q,r,s,t) \in \Omega} \mathbf{F}(q,r,s,t) dq dr ds $$

이러한 $\rho$와 $\mathbf{f}$의 정의를 사용하여서 변분이 정의된 적분식에서 라그랑지안 변수를 모두 없애주면 다음을 얻는다.

$$0=-\wp S = \int_{0}^{T} \int \left( \rho (\mathbf{x},t) \frac{D\mathbf{u}}{Dt} (\mathbf{x},t) - \mathbf{f}(\mathbf{x},t) \right) \cdot \mathbf{v}(\mathbf{x},t) d\mathbf{x} \ dt $$

유체역학을 전공하는 사람 입장에서는 뭔가 익숙한 식이 나오고 있는 것이 느껴질 것이다. 이제 여기서 Hodge decomposition을 생각할 것이다. 사실 논문에서는 호지 분해라고 나와있긴 한데, 내가 호지 분해라는 거를 처음 봐서 내용을 찾아보니까 이 내용은 오히려 헬름홀츠 분해에 해당하는 것 같긴 한데, 아무튼 논문에서는 호지 분해라고 하긴 했다. 

임의의 벡터장은 어떤 gradient와 divergence-free 벡터장의 합으로 나타내질 수 있다는 것이 그 내용이다. 그러면 바로 위의 적분 식에서 괄호 안에 있는 것을 임의이 벡터장이라고 생각하면 다음과 같이 쓸 수 있다.

$$\rho \frac{D\mathbf{u}}{Dt} - \mathbf{f} = - \nabla p +\mathbf{w} $$

그리고 이 때 다음이 성립해야 한다.

$$\nabla \cdot \mathbf{w} = 0$$

사실 이 내용 자체는 쉽게 이해가 된다. 어떤 스칼라의 gradient는 curl-free이고, 임의의 벡터장이 curl-free와 divergence-free로 분해될 수 있다는 것은 매우 직관적이다. 다만 이게 헬름홀츠 분해인지 호지 분해인지 그 이름이 다소 헷갈린다.

여기서 $\mathbf{w} = 0$이라는 것을 보일 건데, $\mathbf{v}$는 $\mathbf{v}(\mathbf{x},0)=\mathbf{v}(\mathbf{x},T)=0$의 제약이 있었기 때문에, $\phi(0)=\phi(T)=0$인 $\phi(t)$에 대해서 다음과 같은 선택이 가능하다.

$$\mathbf{v}(\mathbf{x},t)=\phi(t) \mathbf{w}(\mathbf{x},t)$$

여기서 $\phi$에 대해서도 수많은 선택이 가능한데, $t \in (0,T)$에 대해서 $\phi(t)>0$를 선택하도록 한다. 

위에서 정의한 curl-free와 divergence-free 벡터장을 변분의 적분에 대입하고 $\mathbf{w}$를 사용하면 적분은 다음과 같아진다.

$$0=\int_{0}^{T} \phi(t) \int (- \nabla p(\mathbf{x},t) + \mathbf{w}(\mathbf{x},t))\cdot \mathbf{w}(\mathbf{x},t) d\mathbf{x} \ dt$$

여기서 $\nabla \cdot \mathbf{w}=0$이므로 $\nabla p$에 대한 적분은 부분적분하면 0이 됨을 쉽게 알 수 있고 따라서 적분은 다음과 같아진다.

$$0=\int_{0}^{T} \phi (t) \int \left| \mathbf{w}(\mathbf{x},t) \right| ^2 d\mathbf{x} dt$$

그런데 우리의 선택에 따라서 $\phi$는 $(0,T)$에서 0보다 크기 때문에 $\mathbf{w}=0$이라는 결론을 내릴 수 있다. 이제 마지막으로 incompressible elastic material의 지배 방정식을 완성하는 과정에서, 점성과 관련된 항을 추가하기만 하면 된다. Principle of least action은 보존되는 계에서 적용되는 원리인데, 운동량은 점성에 의해서 소멸되기때문에 principle of least action만으로는 점성에 의한 소산을 유도할 수 없고, 따로 추가해야 하는 것이다. 이 지배 방정식을 모두 모으면 다음과 같다.

$$\rho \left( \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} \right) + \nabla p = \mu \nabla ^2 \mathbf{u} + \mathbf{f} $$

$$ \nabla \cdot \mathbf{u} = 0$$

$$\rho (\mathbf{x} ,t) = \int M(q,r,s) \delta (\mathbf{x} - \mathbf{X} (q,r,s,t) ) dr dr ds $$

$$\mathbf{f} (\mathbf{x} ,t) = \int \mathbf{F} (q,r,s,t) \delta (\mathbf{x} - \mathbf{X}(q,r,s,t) ) dq dr ds $$

$$ \frac{\partial \mathbf{X}}{\partial t } (q,r,s,t) = \mathbf{u} (\mathbf{X}(q,r,s,t),t) = \int \mathbf{u} (\mathbf{x},t) \delta (\mathbf{x} - \mathbf{X}(q,r,s,t)) d \mathbf{x}$$

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

이 방정식에서 $M(q,r,s)$는 주어져야 하는 라그랑지안 질량 밀도이고, $E[\mathbf{X}]$ 또한 주어져야 하는 재료의 탄성 포텐셜 에너지 범함수이다. 이는 탄성에너지 모델에 따라서 달라진다. 그리고 $\wp E / \wp \mathbf{X}$는 $E$의 Frechet 미분이다.