Loading
2015. 6. 19. 00:54 - 성돌

[Gerris] 움직이는 원기둥 주변에 발생하는 유동 시뮬레이션 하기!




이전 튜토리얼에서는 원기둥이 정지해있을 때
원기둥을 향해 일정한 유량을 공급해주었을 때 나타나는 유체역학 현상을 시뮬레이션 하였다.


여기서는 원기둥이 움직일 때 나타나는 현상을 시뮬레이션 해보자.

이 포스팅에서는 원기둥을 일정한 속도로 움직일 것이기에,
일정한 시간이 지난 후 규칙적인 흐름 패턴에 도달했을 때는
원기둥이 가만히 있고 일정한 유속이 오는 경우와 동일한 유체역학적인 현상이 발생해야 한다.

즉, 여기서도 von Karman vortex street가 관측이 되어야 한다는 것이다.

그러나
시뮬레이션 관점에서 보았을 때,
이 경우에는 실린더가 계속 움직이기에 그에 따라 mesh가 같이 움직여주어야 하기에...
사실 CFD관점에서 아주 쉽지는 않은 시뮬레이션이다.

그러나 Gerris는 매우 쉽게 이러한 시뮬레이션을 할 수 있게 해 두었다.
원한다면, 포스팅 맨 밑에 있는 동영상부터 확인해보자.





이 포스팅에서는 이 시뮬레이션을 Gerris로 할 수 있는 명령어를 설명할 것이다.

아래의 코드를 gfs파일로 저장하고, 이에 대해서 설명하도록 하겠다.
(나는 MovingCylinder.gfs란 파일명으로 저장하였다.)


그리고 움직이는 물체에 관한 시뮬레이션에 경우,
물체의 형상을 미리 gts파일로 가지고 있어야 한다.

내가 사용한 원기둥의 형상은 아래에서 다운받을 수 있다.

MyCylinder.gts


이 gts파일은 stl파일로부터 제작이 가능하며,
stl파일은 무료 프로그램인 Blender를 이용해서 만들 수 있다.

그리고 터미널을 이용해서 쉽게 stl파일에서 gts파일로 전환이 가능하다.

이에 대해서는 다른 포스팅에 자세하고 설명을 해놨기에,
찬찬히 읽으면 쉽게 이해할 수 있으리라 생각한다.



많은 부분들이 이전 튜토리얼과 겹치기에 이전 튜토리얼에서 설명하지 않은 부분들만 설명하도록 하겠다.
설명되지 않은 부분은 이전 튜토리얼을 참고하자.

가장 먼저 눈에 띄는 부분은 아래이다.


가장 중요한 것은 이전의 GfsSimulation과는 달리 GfsSimulationMoving이 사용되었다는 점이다.

움직이는 물체가 시뮬레이션에 있을 때, 반드시 이렇게 먼저 선언을 해주어야한다.

그렇지만
다른 포스팅에서 설명된 것처럼 다루는 방정식은
Naiver-Stokes 방정식으로
GfsSimulation 동일하다.

이전 포스팅을 충분히 공부하였다면,
위의 5 4의 숫자가 5개의 기본 정사각형 영역(box)을 설정하고
이 정사각형 영역끼리 4개의 경계를 갖는 걸 나타낸다는 걸 쉽게 알 수 있을 것이다.




여기서는 Refine에 if기능을 응용해보았다.

Refine은 앞서 설명한바와 같이 초기에 mesh를 얼마나 조밀하게 할 건지를 설정는 것인데...

위에서 한 것처럼 내가 원하는 특정 부분은 더 조밀하게 하고,
어떤 부분은 덜 조밀하게 할 수 있다.

불필요한 곳에 mesh를 많이 넣어주는 것은 계산낭비이기 때문이다.

위의 조건으로 mesh를 넣어주면 아래와 같이 초기 mesh가 형성된다.


맨 왼쪽에 원은 이후에 있을 SolidMoving명령어에 의해 설정된 것이다.

즉, 원 (원기둥) 근처에만 조밀하게 mesh를 분포해놓은 것을 알 수 있다.

위에서 명령어로 if명령어가 사용되었는데 아래와 같은 구조이다.



즉, (조건)이 성립하면 A를 실행하고, 성립하지 않는 부분은 B를 실행한다는 것이다.

if안에 있는 fabs는 C명령어로 절대값을 의미하고
&&은 and 조건을 의미한다.

즉,
fabs(y) < 0.2 && x < -0.1 && x > -0.4 조건은...
-0.2 < y < 0.2 and -0.4 < x < -0.1인 영역을 나타내는 것이다.

아직 정의되지 않았지만, 시뮬레이션 영역에서 원점
(x=y=0)
1번 정사각형 영역(box)이 위치한 곳으로...
이후에 맨 왼쪽에 있는 정사각형 영역을 1번 영역으로 설정해 줄 것이다.

정리하면, 위의 조건이 성립하면 29수준의 mesh를 형성하고
(Resol에 9라는 값을 입력해 두었음을 다운 받은 코드에서 확인하자)
그렇지 않을 경우
27수준의 mesh를 형성하는 것이다.




이 명령어가 바로 움직이는 물체를 정의해주는 명령어들이다.

첫번째 줄은 물체와 물체의 초기위치 및 크기를 설정하는명령어이며,
두번째 줄은 물체의 움직임을 설정하는 명령어이다.

우선 첫번째 줄에서 istep은 몇번째 step마다 solid의 움직임을 계산하여 줄 것인지에 관련이 있고,
MyCylinder.gts파일을 불러 읽는 것을 의미한다.

tx와 scale은 앞서 포스팅에서 설명한 바와 같이 solid의 위치와 크기를 지정해주는 것으로
각각 절대값을 의미하는 것이 아니라 박스 한 칸의 길이(Lref)에 대한 비율을 나타낸다.

즉 tx=-0.25은 원점으로부터
Lref × tx 만큼 solid의 위치를 평행이동시키라는 것이고,
이 포스팅에서는
Lref을 변경시키지 않았고 m, kg, s시스템을 사용하기에
default값인
Lref = 1 m이기에...
1 m × -0.25 = -0.25 m만큼 x방향으로 평행이동을 시키라는 말이 된다.

같은 방식으로 원기둥의 지름은 scale명령어에 의해서 0.125 m가 된다.
앞서 Define에 의해서 CylinderDiameter를 0.125로 설정해주었음에 주목하자.

그리고  level은 solid의 움직임의 정밀도를 설정해주는 것이다.


Dirichelet조건은 U의 값 자체를 설정하는 것으로 
U=SphereV로 만들어주는 명령이 된다.
앞서 Define명령어에 의해 SphereV는 1로 설정이 되어있고,
m, s, kg시스템을 사용하기에 이는 1 m/s를 의미한다.

추가로 더 설명을 하자면,
SurfaceBc의 default값은 Neumann=0인 조건으로...
위에서는 y방향 속도인 V에 대한 조건을 입력하지 않았기에...
이는 default값에 의해 자동적으로
Neumann=0으로 설정된다.


Neumann 경계 조건은 표면에서 바깥쪽방향으로 미분하여준 값(normal derivative)을
경계조건으로 설정하는 것을 의미한다.

즉, 이 조건에서는 slip이 발생할 수 있고 굳이 slip을 없애고 싶다면
아래의 조건을 사용해서 V속도를 0으로 만들자.
SurfaceBc V Dirichlet 0



나머지 코드들은 앞서 포스팅에서 모두 다 설명했던 것들이기에 생략하도록 하겠다.

시뮬레이션 실행을 터미널에서 아래의 명령어를 통해서 할 수 있다.


그런데 이번에는 계산량이 많기에 cpu갯수가 충분하다면,
병렬 연산(parallel computing)을 하여 계산 속도를 높여 보도록 하자.
(자세한 설명은 다른 포스팅을 참고하자)

예를 들어, 내 컴퓨터는 cpu가 4개이기에...
4개를 모두 사용하는 짓은 하지 말고... (컴퓨터에 부담이 가니까)
그의 절반인 2개만 사용해보도록 하겠다.

아래의 명령어들을 순차적으로 실행함으로 cpu 2개로 계산을 더 빨리 수행할 수 있다.


우분투의 경우 한 가지 팁을 주자면,
위의 명령어를 복사해서 터미널에서 마우스 휠버튼을 클릭하면 붙여넣기가 된다.

그리고 위의 병렬 연산에 대해 더 자세히 알고싶다면, 링크를 참고하자.




여기까지 과정을 마치고 시뮬레이션을 하면,

아래에 보여주었던 것과 같은 동영상을 얻을 수가 있다.


그리고 앞서 포스팅에서 설명한 바와 같이 항력계수를 계산하면 아래 그래프와 같이 되고...
정상상태(steady-state)에서 Re가 1250일때 Cd=1.58인것으로 계산이 된다.

이는 실제 실험결과와 비교해보았을 때, 타당한 값임을 알 수 있다.