위 그림과 같이
원기둥 주변에 나타나는 유동을 시뮬레이션을 해보면서 Gerris의 사용방법을 익혀보도록 하자.
이러한 유동을 Benard–von Karman vortex street또는 von Karman vortex street이라고 부른다.
우선 Gerris 제작자인 Stephane Popinet이 만든 소스 코드(cylinder.gfs)가 있으며,
이를 다른 포스팅에 설명되어있는 방법으로 실행만 하면 시뮬레이션을 실행할 수 있다.
그러나 개인적으로 이 소스코드에는 불필요한 정보들이 좀 많다고 생각해서...
초보자가 이해하기 쉽고 꼭 필요한 기능만 담은 버젼으로 수정하였다.
아래의 코드를 확장자가 gfs인 파일로 저장하고, 이것을 중심으로 설명을 진행해보겠다.
(나는 SDol_cylinder.gfs라는 파일로 저장하였다.)
이 포스팅에서는 이 코드를 하나하나씩 이해하면서, Gerris코드를 공부해도록 하자.
우선 가장 윗 부분에 위치한 코드부터 하나씩 살펴보자.
위의 코드는 별로 중요한 것이 아니고,
내가 시뮬레이션을 할 때 조절하는 변수들을 Define이라는 함수를 사용해서
시뮬레이션을 시작하기에 앞서 지정한 것이다.
뭐, 간단히 설명하자면 CylinderRadius라는 문자에 0.0625라는 숫자를 입력시키는 것이다.
한 가지 주의할 것은 모든 과정에서 단위를 동일하게 맞춰주어야 한다는 것이다.
나는 이 시뮬레이션에서 m, s, kg의 단위 시스템을 동일하게 사용해줄 것이고,
거의 모든 경우에 이렇게 하는 것이 편하리라 생각한다.
쉽게 말해, 길이에 대한 값은 무조건 m을 단위로 동일하게 할 것이고,
특정 값에 대해서만 mm를 사용한다던지 cm를 사용하는 일이 없어야 한다는 것이다.
즉, 위에서 CylinderRadius = 0.0625는 이후에 원기둥의 반지름을 정의하는데에 사용될 것이며
이는 0.0625 m을 나타낸다.
이렇게 단위를 늘 일치시켜야 하는 이유에 대해 더 자세히 알고 싶다면,
다른 포스팅을 더 참고하도록 하자.
Define은 C에서 사용하는 기능으로...
Gerris에서 기본적으로 제공되는 기능이 아니기에
이를 사용하기 위해서는 이 후에 시뮬레이션을 돌릴 때 명령어에 -m을 추가해주어야 한다.
이에 대해서는 마지막 즈음에 다시 언급하도록 하겠다.
여기서 2과 1이라는 숫자는 시뮬레이션을 실행할 영역 지정과 관련이 있다.
2이라는 숫자는 2개의 '기본' 정사각형 영역을 이용해서 시뮬레이션할 영역을 만든다는 것이다.
(2D시뮬레이션이므로 정사각형이고, 3D였다면 정육면체를 말하는 것임)
그리고 1은 이 2개의 '기본' 정사각형 사이에 1개의 겹치는 경계가 있다는 것을 의미한다.
아래의 그림을 보면 이해가 빠를 것이다.
|
위의 그림이 지금 내가 만들고자 하는 시뮬레이션 영역이다.
위처럼 2개의 정사각형을 옆으로 연결한 후,
왼쪽에 정사각형의 중앙에 원기둥을 위치시킨 후 왼쪽에서 일정한 유량이 공급되는 시뮬레이션이다.
정사각형의 한 변의 길이를 Lref라고 부르는데, 위처럼 특별히 지정하지 않았을 때는 이 길이는 1이다.
특별히 지정하고자 할 때는 PhysicalParams라는 기능을 이용하면 되지만,
이 포스팅에서는 특별히 변경하지 않도록 하겠다.
그리고 앞서 말한 것과 같이, 이 포스팅에서는 m, s, kg의 단위 시스템을 사용할 것이기에
위의 1이라는 숫자를 1 m로 생각하도록 하자.
이 포스팅에서 하고자 하는 시뮬레이션은 아래의 그림을 보면, 이해가 빠를 것이다.
실린더가 왼쪽 정사각형의 중앙에 위치하거나 정사각형들이 저렇게 연결되어 있는 것은
이 후에 나올 코드로 정의하게 된다.
그렇다면, 정사각형 1개만을 영역으로 사용하게 되면 어떻게 정의해줘야 할까?
2 1가 아닌, 1 0이라고 하면 될 것이다.
정사각형 1개의 경우는 경계가 없을테니 말이다.
그리고 다음에 나온 명령어로 GfsSimulation는 비압축성 Navier-Stokes 방정식을 푸는 것을 의미한다.
점성을 0로 설정하였을 경우, 자동적으로 오일러 방정식을 풀게 된다.
정확히 어떤 방정식을 푸는지를 보고 싶다면, 다른 포스팅을 참고하자.
앞으로 Gfs라는 말이 자주 등장할 텐데, 이건 별건 아니고...
Gerris flow solver의 약자이다.
그리고 GfsBox GfsGEdge {} {이라는 명령어는 그냥 시뮬레이션 영역과 경계때문에 선언해준 것인데,
위의 그림에서 표현한 정사각형이 Gfsbox이고 정사각형 사이의 경계가 GfsGEdge이다.
그리고 괄호 사이에 아무것도 없이 표현한 {}은 기본옵션으로 사용한다는 것을 의미하고,
대부분의 경우에 이 안에 옵션을 특별히 지정할 일은 없을 것이다.
그리고 뒤에 따라오는 괄호인 {은 여기서부터 시뮬레이션에 대한 정보를 입력하기 시작한다는 것을 의미하고,
시뮬레이션에 대한 모든 정보를 입력하면 이후에 }를 이용해서 괄호를 닫아줄 것이다.
GfsBox GfsGEdge {} {이라는 명령어는 관습적으로 늘 저렇게 사용하는 걸로 기억해주면 좋을 것이다.
그 다음으로 나오는 명령어는
으로 시뮬레이션하는 시간을 조절하는 기능이다.
그리고 앞으로의 대부분의 Gerris명령어는 GfsTime이라고 쓸 수도 있고,
Gfs를 생략해서 Time이라고 쓸 수도 있다.
그래서 Gerris홈페이지의 설명을 보면 GfsTime이라고 되어 있지만,
대부분의 사람들이 Time이라고 명령어로 쓰지... GfsTime이라고 명령어를 쓰지는 않는다.
'end ='는 언제 시뮬레이션을 끝낼지를 정의하며, 위에서 이미 FinalTime을 7로 입력해두었었다.
그리고 시뮬레이션 사이의 시간 간격(time step)인 dt는 특별히 뭘 더 입력하지 않는 한,
CFL안정조건에 의해 자동적으로 설정된다.
시뮬레이션 사이의 시간 간격이 이 후에 설정할 시뮬레이션 정확도에 의해 자동적으로 변하기 때문에,
여기서는 특별히 이것을 설정할 필요가 없으리라 생각된다.
그리고, 앞서 말한대로 시간에 대한 값은 모두 s를 기준으로 입력하기로 했기에,
7이라는 숫자는 7 s를 의미한다.
그 다음에 오는 명령어는
으로 처음 시작할 때의 mech의 조밀함을 나타낸다.
물론, adaptive mesh refinement기능을 사용하지 않을 경우에는 시뮬레이션이 끝날 때가지
이 mesh상태를 유지하지만....
Gerris의 엄청난 장점인 adaptive mesh refinement기능을 사용하지 않을 이유가 거의 없기에...
이 후에 이 mesh구성은 컴퓨터가 알아서 변경할 것이므로...
Refine기능은 초기 mesh구성으로 생각하면 될 것이다.
Adaptive mesh refinement에 대해서는 뒤에서 더 자세히 다루도록 하겠다.
여기서 6이라는 숫자는 한 박스 안을 26개의 작은 mesh로 채우는 것이다.
아래의 그림이 26수준의 mesh이다.
<Image from very very simple tutorial>
즉, 6을 7로 바꾸면 27수준의 mesh가 되는 것으로 숫자를 하나 올릴 때마다
얼마나 엄청난 차이가 있을지 짐작할 수 있을 것이다.
그 다음에 오는 명령어는
으로 우리의 원기둥을 생성하는 명령어이다.
물론 2D시뮬레이션이기에 원기둥은 원으로 표현이 되며,
위의 식이 원을 기술하는 식임을 쉽게 알 수 있을 것이다.
그리고 원의 반지름은 이미 정의한 CylinderRadius에 저장이 되어있다.
좀 더 다양한 형상을 만드는 것에 대해서는 링크를 참고하도록 하고...
비행기 날개와 같은 좀 더 복잡한 형상에 대해서는 다른 소프트웨어로 stl파일을 제작하고
이걸 gts파일로 바꾸어 준 후 이를 Gerris에서 읽어주는 방식을 사용한다.
Gerris에서는 이 stl파일을 제작하는 데에 있어 무료 소프트웨어인 Blender를 사용할 것을 추천하며,
이를 사용해 gts파일을 제작하는 방법은 링크에 자세히 설명이 되어 있다.
그 다음에 나오는 명령어는
으로 초기조건을 나타낸다.
U는 x방향의 유체의 속도를 가르키며,
위와 같이 설정해둔 것은 x방향의 초기 유속을 1로 설정한다는 것이다.
그리고 앞서 이야기한대로 이 시뮬레이션에서
나는 단위에 대해 m와 s를 기준으로 할 것이기에
여기서 유속 1은 1 m/s를 의미한다.
참고로 이 Init기능을 아예 적어주지 않을 수도 있는데,
이런 경우에는 모든 초기값이 default로 설정된 값으로 입력이 된다.
속도에 대한 default값은 모든 방향의 속도가 0인 것이다.
예를 들어, 위에서 y방향의 속도(V)가 선언되지 않았는데,
이 경우에는 Gerris가 y방향의 초기속도는 0인 것으로 자연스럽게 인지한다.
그 다음에 나오는 명령어는 Gerris가 제공하는 아주 멋진 기능이다.
이 기능은 컴퓨터가 판단하여 mesh가 더 필요한 곳은 더 mesh를 넣어주고,
필요없는 곳은 mesh를 없애는 과정으로...
매우 효율적인 CFD계산이 가능하게 하는 것이다.
이것이 adaptive mesh refinement이다.
이런 복잡하고 중요한 기능을 위의 명령어 한 줄로 사용할 수 있는 것이다.
AdaptVorticity는 vorticity를 기준으로 하여
vorticity의 변화가 크면 알아서 mesh를 더 넣고, 반대의 경우에는 mesh를 줄이는 것이다.
예를 들어 아래와 같이 말이다.
Vorticity말고도 속도, 압력이나 온도같은 것들이
이런 adaptive mesh refinement을 적용하는 기준이 될 수 있으며...
이럴 때는 AdaptGradient기능을 사용한다.
그리고 cmax는 허용가능한 최대 cell cost를 의미하는 것으로,
한 cell에서 내가 지정한 어떤 값의 정규화된 값이...
이 값보다 크면 refinement을 실행하라는 명령어이다.
여기서, 이 값은...
AdaptVorticity에서는 이 값이 vorticity이고,
AdaptGradient에서는 내가 지정한 값의 구배값이 될 것이다.
그렇기에 이 값을 줄이면 줄일수록, 더 정밀한 계산이 이루어지게 될 것이다.
더 자세한 설명이 보고 싶다면,
링크에서 '3.4 절'부분을 자세히 읽어보자.
아무튼 위 명령어에서 istep은 얼마나 자주 adaptive mesh refinement을 수행할 것인지를 나타내며,
'istep = 1'은 매 단계마다 adaptive mesh refinement을 수행하라는 것을 나타낸다.
maxlevel은 컴퓨터가 채워넣을 수 있는 최대한의 mesh조밀도를 나타낸다.
6이라는 숫자는 최대 26수준의 조밀한 mesh까지 컴퓨터가 넣을 수 있다는 걸 나타낸다.
즉, 숫자가 증가할 수록 정밀해지지만... 계산속도는 느려질 것이다.
보다 자세한 설명은 위에 Refine에 대한 설명에 나와 있다.
다시 말하지만, 이 기능은 매우 유용하고 필수적이므로...
꼭 사용하기를 추천한다.
다음은 유체의 점도(viscosity)와 밀도(density)를 정의하는 것이다.
먼저, SourceViscosity이 정의하는 것은 dynamic viscosity이다.
그리고 PhysicalParams안의 alpha는 밀도의 역수를 가리킨다.
역시나 앞서 언급한대로, 단위는 m, s, kg을 기준으로 할 것이기에...
점성계수 0.78125은 0.78125 kg/m∙s를 의미하고,
밀도의 역수 1e-3은 103 kg/m3의 밀도를 의미한다.
여기서 e는 10의 대한 지수를 나타내며, 프로그래밍에서 이와같이 자주 쓰인다.
예를 들어, 1e-1은 10-1이다.
다음은 시뮬레이션이 어떻게 진행되어 가는 지...
시간에 대한 정보를 표시하라는 명령어이다.
이것은 시뮬레이션을 얼마나 오랫동안 하였는지...
또는 얼마나 시뮬레이션이 남았을 지 짐작하기 위해서 사용한다.
이 명령어도 그냥 복사, 붙여넣기 수준으로 해서 사용하면 좋을 것이다.
'istep = 10'은 10번 계산을 할 때마다 결과를 출력하는 것을 의미한다.
아래처럼 말이다.
일단 중요한 것만 설명을 하면...
step은 지금까지 계산 횟수를 의미하고,
t는 계산완료된 시뮬레이션의 시간...
dt는 step사이의 시간 차이...
real은 지금까지 실제로 계산하는 데 걸린 시간이다.
위에서 step을 제외하고, 모든 단위는 s이다.
즉, 위에서 30단계까지 계산해서...
시간 0.05757 s까지 시뮬레이션을 완료하였고,
여기까지 계산하는데 1 s정도가 소요되었음을 의미한다.
그 다음 기능은 결과를 동영상으로 만들어주는 기능이다.
ppm이란 사진 이미지를 의미하는 데...
ppm2mpeg라는 기능을 사용해서 이미지를 동영상으로 만들어주는 것을 의미한다.
이렇게 동영상을 만드는 기능을 사용하기 위해서는 ffmpeg라는 기능이 우분투에 설치되어 있어야 한다.
만약 설치되지 않았다면, 다음의 명령어를 이용해서 설치하도록 하자.
아무튼 명령어는 2단계마다 (istep =2) vorticity에 대한 동영상(v=Vorticity)을
최대갑이 10이고 (빨간색으로 표시) 최소값으로 -10 (파란색으로 표시)으로 하여 제작하라는 명령어이다.
위에서 v는 vorticity외의 다른 물리량이 될 수도 있다.
즉, 다른 물리량이 시간에 따라 변화하는 동영상을 제작할 수 있다는 것이다.
예를 들어, 속도의 절대값을 v = Velocity이며, 수평방향의 속도만 보고 싶으면 v = U이다.
다른 물리량도 이와 같이 대입하면 된다.
그 다음 명령어도 내가 아주 좋아하는 명령어로...
내가 지정한 solid(원기둥)에 얼만큼의 힘이 가해지는 지를 자동으로 계산해서,
정해준 파일에 저장하는 거다.
즉, 위의 명령어는 매 단계마다 계산을 하여 (istep=1)
force.dat이라는 파일에 solid가 받는 힘을 계산하는 명령어가 될 것이다.
이 파일을 엑셀 비슷한 프로그램으로 열어주면, 아래와 같이 나오는데...
여기서 1열은 위에 적힌 것처럼 시간(T)을 의미하고,
2, 3, 4열은 압력에 의한 힘을 각각 x, y, z방향으로 분해해서 표시하는 것이고,
5, 6, 7열은 점성에 의한 힘을 각각 x, y, z방향으로 분해해서 표시하는 것이고,
나머지는 비슷한 방법으로 회전력을 계산한 것이다.
나는 이후에 시뮬레이션 결과를 사용해서 항력계수(drag coefficient)를 계산할 생각이기에,
x방향의 힘인 2와 5번째 열 정보만이 필요할 것이다.
더불어 시간에 대한 정보인 1열도 필요할 것이다.
그리고 아래는 시뮬레이션을 한 모든 결과를 저장한 파일을 만들기 위한 명령어로...
이것에 대해서 설명하기전에 맨 마지막 줄에 }로
앞서 'GfsSimulation GfsBox GfsGEdge'줄에서 열어주었던 {괄호를 닫아주는 점에 주목하자.
그리고 첫번째와 두번째 명령어는 각각 100번의 계산이 끝날 때마다,
'Output-계산횟수.gfs'과 'Output-계산횟수.txt'파일에 결과를 저장하는 것이다.
gfs파일로 저장된 시뮬레이션 결과 파일은 아래의 명령어를 이용해서
gfsview프로그램으로 볼 수 있다.
예를 들어, step 3000번째의 파일인 Output-3000.gfs를 읽어서
압력분포를 그려보면 아래와 같다.
gfsviewer는 쉬우며 상당히 다양한 기능을 제공하니,
시간을 내서 한 번 살펴보도록 하자.
그리고 같은 3000번의 step에서 저장된 txt파일인 Output-3000.txt를
엑셀 프로그램으로 열어보면 아래와 같다.
적혀있는대로 1, 2, 3번째 열은 x, y, z좌표..
그리고 4,5,6,7열은 각각 좌표에 해당하는 압력, MAC projection 압력, x방향과 y방향 속도를 가리킨다.
드디어 마지막 명령어이다.
이 명령어들은 위에서 말한 시뮬레이션 사각형 또는 박스의
경계조건을(boundary condition) 설정하고 어떻게 연결이 되어있는 지를 표시하기위함이다.
설명하기 전에 특별하게 변경하지 않았을 경우,
1이라고 명명된 박스의 중앙에 좌표축의 원점이 있다는 점을 명심해두자.
그리고 두번째 박스는 첫번째 박스의 오른쪽에 있다는 것을 '1 2 right'이라는 명령어로 표현하였다.
이외에도 left, top, bottom명령어를 사용하여 두번째 박스를 왼쪽, 위, 또는 아래에 위치시킬 수 있다.
그리고 첫번째 GfsBox는 첫번째 박스의 경계조건으로 왼쪽면에 디리클레(Dirichlet) 경계조건에 의해
x방향 유속 (U)을 1으로(InletSpeed에 입력됨) 설정해주었다.
역시나 이는 단위 시스템을 m과 s로 정했기에, 1 m/s를 의미한다.
디리클레 경계조건은 단순하게 그 값 자체를 경계조건으로 사용하는 것을 의미하며,
노이만(Neumann) 경계조건은 값의 구배값이 경계조건으로 사용되는 것 의미하지만
이 예제에서는 사용되지 않았다.
이 포스팅에서는 일정한 유속을 공급해주었지만,
유속을 시간(t)에 대한 함수로 표시해주면
임의의 어떤 유속의 변화도 경계조건으로 설정이 가능하다.
그리고 두번째 줄의 GfsBox는 오른쪽 면으로 자유롭게 유체가 흘러갈 수 있게 한 것이다.
그렇다면, 아무것도 지정하지 않은 면들의 경계조건은 어떻게 된 것인가?
예를 들어, 첫번째 박스의 윗부분과 같은 곳 말이다.
Gerris에서는 경계조건을 입력하지 않을경우,
default값으로 고체 표면이 거기 있다고 생각하고 no-slip 경계조건을 할당한다.
시뮬레이션을 실행하기 위해서는 아래의 명령어로 실행한다.
기본적인 시뮬레이션의 경우는 -m명령어가 필요가 없지만,
이 예제에서는 Define명령어를 사용해주었기에 -m를 명령어에 포함시켜주어야 한다.
계산을 마치면 아래와 같은 vorticity에 대한 동영상을 얻을 수 있고,
(내 컴퓨터로 20분정도 걸림...)
동영상에서 빨간색은 vorticity가 양의 값이고, 파란색은 vorticity가 음의 값을 의미한다.
양수인 vorticity는 반시계방향으로 유체가 회전한다는 걸 의미하며,
음의 vorticity는 시계방향의 회전을 나타낸다.
팁을 주자면...
오른손을 이용하여 오른손 엄지손가락이 화면에서 나오는 방향을 가르키게하면,
나머지 손가락들이 반시계방향인데...
이걸 양의 방향을 vorticity를 기억하는 데에 사용할 수 있다.
그리고, 아래 그래프와 같이 시간이 변함에 따라 항력계수가 변화함을 볼 수 있다.
이에 따라 계산된 항력계수는 1.5으로,
상응하는 레이놀즈 수(Re=UL/ν=160)에 대한 실제 실험 결과와 유사하다는 것을 알 수 있다
(링크와 비교!).
항력계수를 계산할 때는 힘은 본 예제에서는 x방향의 힘만을 고려할 것이다.
즉, force.dat파일에 저장된 2와 5번째 열을 합친 값이다.
그리고 아래의 공식을 이용하면, 2D상황에 대한 항력계수를 계산할 수 있다.
최대한 자세히 적느라 힘들었다...
그래도 이 튜토리얼이 우리 나라 사람으로 하여금,
Gerris에 대해 보다 쉽게 접근할 수 있게하여 수행하는 일에 도움이 될 수 있기를 바란다.
'Software > CFD' 카테고리의 다른 글
| [Gerris] 시뮬레이션에 사용할 형상 제작 (gts파일 만들기) (6) | 2015.06.18 |
|---|---|
| Gerris에서 푸는 Navier-Stokes방정식 (0) | 2015.06.03 |
| [Gerris] 시뮬레이션 실행하기! (0) | 2014.09.27 |
| Gerris 설치하기 (6) | 2014.09.27 |
| 무료 유체역학 시뮬레이션 (CFD) 프로그램 Gerris! (1) | 2014.09.24 |