지난 글의 힘과 에너지의 관계에 이어서, 이번에는 이런 힘이 가해질 때 실제 물체[여기서는 비압축성 탄성체]가 어떻게 움직이는지를 소개하는 Equations of motion이다. 힘이 가해질 때 운동을 기술하는 방법이 몇 가지 있을 수 있는데, 그거는 살짝 뒤에 소개하기로 하고 먼저 비압축성 탄성체라는 말을 수학적으로 표현하면서 시작한다.
Material의 label을 나타내는 라그랑지안 좌표계 $(q,r,s)$에서 $\mathbf{X}(q,r,s,t)$가 그 label에 해당하는 실제 점의 물리적 위치일 때, 라그랑지안 좌표 공간에서 다음과 같은 구간에 의해 bound되는 아주 작은 직육면체를 생각하자.
$$\left[q, q+dq \right] \times \left[r, r+dr \right] \times \left[s, s+ds \right] $$
이 때 라그랑지안 공간에서의 이 직육면체의 부피는 아주 쉽게 다음과 같이 쓸 수 있다.
$$dqdrds$$
그런데 이 라그랑지안 공간에서의 부피는, 실제 물리적 공간[오일러리안 공간]에서 다르게 표현된다. 간단하게 생각하면 라그랑지안 공간에서의 변위 $dq$에 대응하는 오일러리안 공간에서의 변위는 다음과 같다.
$$\frac{\partial \mathbf{X}}{\partial q} dq$$
즉, $r$과 $s$ 방향까지 생각한다면 라그랑지안 공간에서의 직육면체는 오일러리안 공간에서 다음의 세 벡터에 의해서 생기는 평행육면체인 것이다.
$$ \frac{\partial \mathbf{X}}{\partial q} dq, \ \frac{\partial \mathbf{X}}{\partial r} dr, \ \frac{\partial \mathbf{X}}{\partial s} ds $$
그리고 어떤 세 벡터 $\vec{a}$, $\vec{b}$, $\vec{c}$에 의해 생기는 평행육면체의 부피는 다음과 같은 것은 쉽게 알 수 있다.
$$\left| \vec{a} \cdot \left(\vec{b} \times \vec{c} \right) \right|$$
따라서 오일러리안 공간에서의 미소부피는 위의 미분으로 표현된 세 벡터의 스칼라 삼중곱으로 다음과 같이 표현할 수 있다.
$$dV = \left| \frac{\partial \mathbf{X}}{\partial q} dq \cdot \left( \frac{\partial \mathbf{X}}{\partial r} dr \times \frac{\partial \mathbf{X}}{\partial s} ds \right) \right| = \left| \frac{\partial \mathbf{X}}{\partial q} \cdot \left( \frac{\partial \mathbf{X}}{\partial r} \times \frac{\partial \mathbf{X}}{\partial s} \right) \right| dq dr ds$$
여기서 마지막의 $dq,\ dr,\ ds$는 스칼라니까 절대값 밖으로 뺀 것이다. 그리고 절대값 기호 안에 벡터의 스칼라 삼중곱이 있는데, 벡터로 구성된 행렬의 행렬식은 스칼라 삼중곱과 같다는 것은 그냥 기본이다. 따라서 미소부피와 자코비안의 관계를 다음과 같이 쓸 수 있다.
$$J(q,r,s,t) = \det \left( \frac{\partial \mathbf{X}}{\partial q}, \frac{\partial \mathbf{X}}{\partial r}, \frac{\partial \mathbf{X}}{\partial s} \right) $$
$$dV = J \ dq dr ds$$
그리고 이 때 탄성체의 전체 부피는 다음과 같다.
$$V = \int dV = \int_{Q} J(q,r,s,t) dr dr ds \quad \mathrm{where} \quad (q,r,s) \in Q$$
그런데 우리가 다루는 탄성체는 incompressible하니까 어떤 Q의 선택에 대해서도 $dV/dt = 0$이어야 한다. 따라서 이는 다음으로 이어진다.
$$\frac{\partial J}{\partial t} = 0$$
여기까지는 그냥 단순한 incompressibility를 라그랑지안 좌표에서 수학적으로 표현한 것이고, 이제부터 실제 운동 방정식을 유도하는 과정이다.
지금까지 어떤 형상 $\mathbf{X}(q,r,s,t)$에 대해서 탄성에너지 $E[\mathbf{X}]$를 정의할 수 있었고, 에너지의 변분을 통해 힘 $\mathbf{F}$를 얻을 수 있었다. 그런데 이러한 힘을 받는 비압축성 탄성체가 어떻게 움직이는지를 기술하는 운동방정식을 아직 모른다.
힘과 운동 사이를 잇는 방법은 뉴턴 역학적으로 운동량 보존을 직접 쓰는 방법이 가장 간단하다. 사실 이 글을 어느 정도 이해할 수 있는 사람은 부연설명을 안해도 바로 무슨 소리인지 이해하겠지만, 모든 힘을 다 쓰고 가속도의 모든 성분을 쓴다는 것이다. 대표적인 예시로 나비에 스토크스 방정식도 유체에 작용하는 모든 힘과 유체의 가속도의 모든 성분을 그냥 쓴 것이다.
그런데 지금 다루는 유체와 고체가 같이 있는 경우에는 이 방법이 어렵다. 왜냐하면 유체는 오일러리안 좌표에서 정의되고 탄성체는 라그랑지안 좌표로 정의되는데, 이 두 좌표를 같이 다루기가 어려운 것이다. 일단 대표적으로 두 좌표계에서는 가속도를 정의하는 방법 자체가 다르다. 그래서 IBM을 고안한 Charels Peskin 아저씨는 에너지의 관점에서 운동을 기술하는 방법을 선택했다.
이는 라그랑주 역학적인 접근을 한다는 것인데, 라그랑주 역학의 핵심 철학인 해밀턴의 원리[또는 최소 작용의 원리]는 자연은 항상 작용이 최소화되는 방향으로 운동하도록 evolve한다는 것이다. 이는 수학적으로 다음과 같이 표현된다. 라그랑지안 $L(t)$에 대해서 시간 $(0,T)$ 동안 작용[action] $S$는 다음과 같이 정의된다.
$$S= \int_{0}^{T}{L(t)}dt$$
이 $S$가 최소화되도록 운동한다는 것이고, 이 때 라그랑지안 $L(t)$는 운동에너지에서 포텐셜에너지를 뺀 스칼라로 정의된다. 운동에너지를 $T$라고 하고 포텐셜에너지를 $U$라고 하면 간단하게 다음과 같다.
$$L(t) = T-U$$
그러면 이제 운동에너지랑 포텐셜에너지를 쓰면 되는데, 먼저 mass density $M(q,r,s)$에 대해서 운동에너지는 다음과 같은 것이 쉽게 이해된다.
$$T=\frac{1}{2} \int M(q,r,s) \left| \frac{\partial \mathbf{X}}{\partial t} (q,r,s,t) \right|^2 dq dr ds$$
그리고 여기서 다루는 포텐셜에너지는 탄성에너지로 이전 포스트에서 정의했었다. 그러면 라그랑지안은 다음과 같이 쓸 수 있다.
$$L(t) = \frac{1}{2} \int M(q,r,s) \left| \frac{\partial \mathbf{X}}{\partial t} \right| ^2 dq dr ds - E[\mathbf{X}(,,,t)]$$
여기서 $S$가 최소가 되도록 하는 경우를 찾기 위해서 $S$에 변분을 취해서 0이 되는 경우를 찾으면 된다. 그리고 변분과 시간미분은 서로 관련이 없다는 것을 생각하면 쭉 쉽게 이해할 수 있다.
$$\wp S = \int_{0}^{T} \wp {L(t)} dt $$
이제 운동에너지와 포텐셜에너지의 변분을 각각 계산하면 되는데, 먼저 운동에너지의 변분은 다음과 같다.
$$\wp T = \int M(q,r,s) \frac{\partial \mathbf{X}}{\partial t} \cdot \wp \left( \frac{\partial \mathbf{X}}{\partial t} \right) = \int M(q,r,s) \frac{\partial \mathbf{X}}{\partial t} \cdot \frac{\partial}{\partial t} (\wp \mathbf{X} )$$
그리고 탄성에너지의 변분은 탄성력의 음수로 이전 포스트에서 보였다. 따라서 action의 변분을 쓰면 다음과 같다.
$$\wp S = \int_{0}^{T} { \int M(q,r,s) \frac{\partial \mathbf{X}}{\partial t} \cdot \frac{\partial}{\partial t} (\wp \mathbf{X} )} dq dr ds \ dt +\int_{0}^{T} \int {\mathbf{F} \cdot \wp \mathbf{X} }dq dr ds \ dt$$
여기서 우변의 첫 번째 항을 시간에만 대해서 부분적분하면
$$\int_{0}^{T} M(q,r,s) \frac{\partial \mathbf{X}}{\partial t} \cdot \frac{\partial}{\partial t} (\wp \mathbf{X})dt = \left[ M(q,r,s) \frac{\partial \mathbf{X}}{\partial t} \cdot \wp \mathbf{X} \right]_0^T - \int_{0}^{T} M(q,r,s) \frac{\partial ^2 \mathbf{X}}{\partial t^2} \cdot \wp \mathbf{X} dt $$
그런데 마지막 상태와 초기 상태의 perturbation은 없으므로
$$\wp \mathbf{X}(T) = \wp \mathbf{X}(0) =0$$
이고, 따라서 다음이 성립한다.
$$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$$
이 적분식은, 모든 가능한 임의의 perturbation 변위 $\wp \mathbf{X}$에 대해서 적분이 0이 되도록 운동한다는 것을 의미한다. 여기서 perturbation 변위 $\wp \mathbf{X}$에 아무런 제약이 없다면 다음과 같은 뉴턴의 제 2법칙으로 귀결된다.
$$ M \frac{\partial ^2 \mathbf{X}}{\partial t^2} = \mathbf{F} $$
그런데 $\wp \mathbf{X}$가 아무런 제약이 없지 않고, 이 글의 처음에서 설명한 incompressibility를 만족해야 한다는 제약이 있다. 그래서 이 최종 적분식에서 $dJ/dt=0$을 고려해야 하는데, 이는 3 by 3 행렬의 시간미분을 적분식에 반영해야 해서 다소 복잡한 면이 있다. 그런데 유체역학을 조금 안다면, 유체의 incompressibility는 오일러리안 속도 $\mathbf{u}$에 대해서 $\nabla \cdot \mathbf{u}=0$으로 간단하게 정의할 수 있다는 것을 안다. 그래서 라그랑지안으로 표현된 변분의 적분에서 라그랑지안으로 표현된 incompressibility 제약을 고려하는 대신, Charles Peskin의 아이디어는 변분 자체를 오일러리안으로 표현하는 것이다. 그러면 incompressibility를 고려하기 쉬워지기 때문이다. 그래서 오일러리안 좌표에서의 변수를 도입하는데, 한 가지 새로운 것은 perturbation $\wp \mathbf{X}$에 대응하는 가짜 속도 $\mathbf{v}$를 도입하는 것이다. 다음과 같이 말이다.
$$\mathbf{u} \left( \mathbf{X}(q,r,s,t), t \right)=\frac{\partial \mathbf{X}}{\partial t} (q,r,s,t) $$
$$\mathbf{v} \left( \mathbf{X}(q,r,s,t), t \right) = \wp \mathbf{X} (q,r,s,t) $$
여기서 $\mathbf{u}$는 그저 단순히 실제 오일러리안 속도이고 $\mathbf{v}$는 material point가 시간 $t$에서 $\mathbf{X}$ 위치에 있을 때 겪는 perturbation velocity이다.
그리고 오일러리안 실제 속도에 대해서 유체역학에서 익숙한 물질미분과 가속도를 사용할 것이다. 다음과 같이 말이다.
$$\frac{D \mathbf{u}}{Dt} = \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} $$
$$ \frac{D \mathbf{u}}{Dt} \left( \mathbf{X}(q,r,s,t), t \right) = \frac{\partial ^2 \mathbf{X}}{\partial t^2} (q,r,s,t) $$
여기서부터는, 결국 $\wp J=0$이 $\nabla \cdot \mathbf{v} =0$이랑 같다는 것을 보이는 과정이다. 그런데 중간중간 행렬에 대한 identity를 사용한 유도가 나오는데, 목적만 잘 기억하고 따라가면 쉽게 이해할 수 있다.
먼저 정방행렬 $A$에 대해서 다음의 identity를 사용할 것이다.
$$ \wp \log \left( \det (A) \right)= \mathrm{tr} \left( (\wp A) A^{-1} \right)$$
이를 유도하기 위해서 먼저 여인수행렬의 전치행렬인 수반행렬을 통해 역행렬을 구하는 것부터 시작한다.
$$\left( A^{-1} \right)_{ji} = \frac{1}{\det A } \left( \mathrm{adj}\ A \right)_{ji} = \frac{1}{\det A} C_{ij}$$
여기서 $\mathrm{adj}$ 는 수반행렬, adjoint를 말하고 $C_{ij}$는 여인수, cofactor를 말한다. 여기까지는 그냥 기본인데, 조금만 생각하면 행렬식이 어떤 특정 행렬 성분에 대해서 선형이고 어떤 특정 성분을 고정하면 성분과 여인수의 곱에 그 성분과 관련 없는 나머지 값들의 합이라는 것을 쉽게 알 수 있다. 좀 말이 복잡한데, 식으로 보면 쉽다.
$$\det A = A_{ij} C_{ij} + D$$
라는 것인데, 이는 어떤 특정한 성분 $A_{ij}$를 고정시켰을 때, 이 성분과 여인수의 곱과 $D$의 합이라는 것이고 여기서 $D$는 $A_{ij}$와 곱해지지 않는 다른 값들이다. 그렇다면 다음까지도 이해할 수 있다.
$$\frac{\partial \det A}{\partial A_{ij}} = C_{ij}$$
이므로 결국 다음이 성립한다.
$$ \left( A^{-1} \right)_{ji} = \frac{1}{\det A} \frac{\partial \det A}{\partial A_{ij}}$$
이제 원하는 identity를 유도하기 시작하면, 먼저 $\log (\det A)$의 변분을 볼 것인데, 어떤 다변수함수 $f$의 미소 변화는 $x_1, \ x_2, \ \cdots , \ x_n$의 모든 성분에 대한 편미분과 그 방향의 미소 변형의 곱의 합으로 나타낼 수 있는 것은 쉽게 이해할 수 있을 것이다. 수식으로 보면 더 쉬운데,
$$df = \frac{\partial f}{\partial x_1} dx_1 + \frac{\partial f}{\partial x_2} dx_2 + \cdots \frac{\partial f}{\partial x_n}dx_n$$
이라는 것이다. 이거랑 똑같이 생각하면 다음을 바로 이해할 수 있다.
$$\wp \log (\det A) = \sum_{ij} (\wp A)_{ij} \frac{\partial \log (\det A)}{\partial A_{ij}}$$
그리고 우변의 편미분은 로그미분이니까 다음과 같이 쉽게 처리할 수 있고,
$$ \sum_{ij} (\wp A)_{ij} \frac{\partial \log (\det A)}{\partial A_{ij}} = \sum_{ij} (\wp A)_{ij} \frac{1}{\det A} \frac{\partial \det A}{\partial A_{ij}}$$
여기의 우변에다가 위에서 짧게 설명했던 여인수[행렬식의 편미분]를 통한 역행렬의 계산을 적용하면,
$$ \sum_{ij} (\wp A)_{ij} \frac{1}{\det A} \frac{\partial \det A}{\partial A_{ij}} = \sum_{ij} (\wp A)_{ij} \left( A^{-1} \right) _{ji}$$
가 되고 여기서 우변은 정확히 trace의 정의이다.
$$ \sum_{ij} (\wp A)_{ij} (A^{-1}) _{ji} = \mathrm{tr} \left( (\wp A) A^{-1} \right) $$
따라서 최종적으로 다음이 성립한다.
$$ \wp \log \left( \det (A) \right)= \mathrm{tr} \left( (\wp A) A^{-1} \right)$$
이제 원래 목적이었던 incompressiblity를 단순히 하기 위해서 다음의 벡터를 도입한다.
$$\mathbf{a} = (q,r,s)$$
그렇다면 다음과 같이 쓸 수 있다.
$$J = \det \left( \frac{\partial \mathbf{X}}{\partial \mathbf{a}} \right)$$
그리고 앞서 $\mathbf{v}$를 정의했던
$$\mathbf{v} \left( \mathbf{X}(q,r,s,t),t \right) = \wp \mathbf{X} (q,r,s,t)$$
의 양변을 $\mathbf{a}$로 미분하면, $\mathbf{v}$는 오일러리안 좌표에서 정의된 속도이므로 chain rule에 의해서 오일러리안 좌표 $\mathbf{x}$에 대해서 다음과 같이 표현된다.
$$\frac{\partial \wp \mathbf{X}}{\partial \mathbf{a}} = \frac{\partial \mathbf{v}}{\partial \mathbf{x}} \frac{\partial \mathbf{X}}{\partial \mathbf{a}} $$
여기서 소문자 $\mathbf{x}$는 오일러리안 좌표에서의 공간 변수를 의미하고 대문자 $\mathbf{X}$는 material point가 실제로 있는 물리적 위치의 좌표를 의미한다. 그런데 여기서 변분을 먼저 하든 $\mathbf{a}$에 대해서 먼저 미분을 하든 상관이 없으니까 다음과 같이 쓸 수 있다.
$$\left( \wp \frac{\partial \mathbf{X}}{\partial \mathbf{a}} \right) \left( \frac{\partial \mathbf{X}}{\partial \mathbf{a}} \right) ^{-1} = \frac{\partial \mathbf{v}}{\partial \mathbf{x}} $$
여기서 양변에 trace를 취하면 아까 유도했던 행렬식의 로그의 변분과 트레이스의 관계에 의해서 다음이 성립한다.
$$\wp \log \left( J(q,r,s) \right) = \nabla \cdot \mathbf{v} \left( \mathbf{X}(q,r,s,t),t \right)$$
따라서 라그랑지안에서 정의했던 incompressibility $(\wp J = 0)$은 다음과 같음이 보여졌다.
$$\nabla \cdot \mathbf{v} = 0$$
이제 비압축성 조건을 오일러리안으로 성공적으로 옮겼으니까 작용의 변분 $\wp S$를 오일러리안으로 옮기면 되는데, 원래 그거까지 다 해서 Equation of Motion에 대한 내용을 close하는 게 두 번째 글의 예상 scope였는데, 쓰다보니까 너무 길어져서 여기까지만 하고 세 번째 글에서 본격적으로 라그랑지안으로 표현된 작용을 오일러리안으로 옮겨서 Equation of Motion까지 마무리하겠다.
'CFD' 카테고리의 다른 글
| Immersed Boundary Method - 3. Equations of Motion (0) | 2026.02.19 |
|---|---|
| Immersed Boundary Method - 1. 간단한 힘과 에너지의 관계 (1) | 2026.01.26 |
| Windows 10에 OpenFOAM 설치 (1) | 2024.07.03 |
