Natural Convection CFD with Elmer

Elmer tutorial 문서의 자연 대류 현상을 모사하는 예제를 참고하여, 따라해 보았다.

모델링은 Gmsh에서 했다.  *.geo 파일 내용은 아래와 같다.

cl__1 = 1;
Point(1) = {0, 0, 0, 1};
Point(2) = {100, 0, 0, 1};
Point(3) = {300, 0, 0, 1};
Point(4) = {400, 0, 0, 1};
Point(5) = {400, 400, 0, 1};
Point(6) = {0, 400, 0, 1};
Line(1) = {2, 3};
Transfinite Line {1} = 50Using Progression 1;
Line(2) = {3, 4};
Transfinite Line {2} = 10Using Progression 1;
Line(3) = {4, 5};
Transfinite Line {3} = 30Using Progression 1;
Line(4) = {5, 6};
Transfinite Line {4} = 30Using Progression 1;
Line(5) = {6, 1};
Transfinite Line {5} = 30Using Progression 1;
Line(6) = {1, 2};
Transfinite Line {6} = 10Using Progression 1;
Line Loop(8) = {4, 5, 6, 1, 2, 3};
Plane Surface(8) = {8};
Recombine Surface {8};
Physical Line(9) = {1};
Physical Line(10) = {2, 3, 4, 5, 6};
Physical Surface(11) = {8};

간단히 말해 2차원 단면으로 했고, 0.4m 사이즈의 변으로 된 정사각형 도메인이다.
열원은 바닥 가운데 중에서 0.2m 영역으로 할 수 있도록 했다.
Gmsh로  만든 *.geo 파일을 2D매쉬를 형성해서 *.msh파일을 생성한 후,
ElmerGrid에서 이것을 변환하였고, ElmerGUI를 이용하여 빠르게 case.sif 입력파일을 생성했다.

만들어진 case.sif 파일은 아래와 같다.

Header
  CHECK KEYWORDS Warn
  Mesh DB "." "."
  Include Path ""
  Results Directory ""
End

Simulation
  Max Output Level = 5
  Coordinate System = Cartesian
  Coordinate Mapping(3) = 1 2 3
  Simulation Type = Transient
  Steady State Max Iterations = 1
  Output Intervals = 1
  Timestepping Method = BDF
  BDF Order = 1
  Timestep intervals = 500
  Timestep Sizes = 0.1
  Solver Input File = case.sif
  Post File = case.ep
Coordinate Scaling = Real 0.001
End

Constants
  Gravity(4) = 0 -1 0 9.82
  Stefan Boltzmann = 5.67e-08
  Permittivity of Vacuum = 8.8542e-12
  Boltzmann Constant = 1.3807e-23
  Unit Charge = 1.602e-19
End

Body 1
  Target Bodies(1) = 11
  Name = "Body 1"
  Equation = 1
  Material = 1
  Body Force = 1
  Initial condition = 1
End

Solver 1
  Equation = Result Output
  Single Precision = True
  Procedure = "ResultOutputSolve" "ResultOutputSolver"
  Save Geometry Ids = True
  Output File Name = case
  Output Format = Vtu
  Exec Solver = After Timestep
End

Solver 3
  Equation = Heat Equation
  Variable = Temperature
  Procedure = "HeatSolve" "HeatSolver"
  Exec Solver = Always
  Stabilize = True
  Bubbles = False
  Lumped Mass Matrix = False
  Optimize Bandwidth = True
  Steady State Convergence Tolerance = 1.0e-5
  Nonlinear System Convergence Tolerance = 1.0e-7
  Nonlinear System Max Iterations = 20
  Nonlinear System Newton After Iterations = 3
  Nonlinear System Newton After Tolerance = 1.0e-3
  Nonlinear System Relaxation Factor = 1
  Linear System Solver = Iterative
  Linear System Iterative Method = BiCGStab
  Linear System Max Iterations = 500
  Linear System Convergence Tolerance = 1.0e-10
  BiCGstabl polynomial degree = 2
  Linear System Preconditioning = Diagonal
  Linear System ILUT Tolerance = 1.0e-3
  Linear System Abort Not Converged = False
  Linear System Residual Output = 1
  Linear System Precondition Recompute = 1
End

Solver 2
  Equation = Navier-Stokes
  Procedure = "FlowSolve" "FlowSolver"
  Variable = Flow Solution[Velocity:2 Pressure:1]
  Exec Solver = Always
  Stabilize = True
  Bubbles = False
  Lumped Mass Matrix = False
  Optimize Bandwidth = True
  Steady State Convergence Tolerance = 1.0e-5
  Nonlinear System Convergence Tolerance = 1.0e-7
  Nonlinear System Max Iterations = 20
  Nonlinear System Newton After Iterations = 3
  Nonlinear System Newton After Tolerance = 1.0e-3
  Nonlinear System Relaxation Factor = 1
  Linear System Solver = Iterative
  Linear System Iterative Method = BiCGStab
  Linear System Max Iterations = 500
  Linear System Convergence Tolerance = 1.0e-7
  BiCGstabl polynomial degree = 2
  Linear System Preconditioning = Diagonal
  Linear System ILUT Tolerance = 1.0e-3
  Linear System Abort Not Converged = False
  Linear System Residual Output = 1
  Linear System Precondition Recompute = 1
End

Equation 1
  Name = "Convection"
  Convection = Computed
  Active Solvers(3) = 1 3 2
End

Material 1
  Name = "Air (room temperature)"
  Heat Conductivity = 0.0257
  Heat Capacity = 1005.0
  Density = 1.205
  Relative Permittivity = 1.00059
  Viscosity = 1.983e-5
  Sound speed = 343.0
  Heat expansion Coefficient = 3.43e-3
End

Body Force 1
  Name = "Buoyancy"
  Boussinesq = True
End

Initial Condition 1
  Name = "Initial"
  Velocity 2 = 1e-9
  Velocity 1 = 0
  Temperature = 293
End

Boundary Condition 1
  Target Boundaries(1) = 10 
  Name = "Wall"
  Noslip wall BC = True
  Temperature = 293
End

Boundary Condition 2
  Target Boundaries(1) = 9 
  Name = "Plate"
  Noslip wall BC = True
  Temperature = 353
End

처음에는 수렴에 실패했었기 때문에, 몇가지 조건들을 수정해 가면서 감을 잡아 나갔다.
예를 들어, Time Step을 어느정도로 끊어줘야 적절한지 Pysical Sense가 좀 부족한 상태라서, 처음에는 2초로 줬다가 수렴에 실패했다.  그 다음에 1초로 줄여서 수렴에 성공하였고, 십여개의 해석 프레임을 얻어내서 그것을 Paraview로 확인해서 Transient Speed를 보고 대충 감을 잡은 후 0.1초로 확정했다.

그리고 정밀도가 좀 떨어지더라도 빠르게 수렴해서 신속하게 결과를 우선 확인해 볼 수 있도록 하기 위하여, Navier-Stokes Eq Solver 설정에서 Linear System Convergence Tolerance 값을 1.0e-7 으로 크게 키워주었다.  이렇게 하면 유체해석 부분에서 수렴이 좀 덜 된 상태라도 그냥 다음 스텝으로 넘어가게 할 것이다.

바운더리 조건은, 전부 No-Slip으로 해 줬고, 발열원은 80도씨, 벽면은 20도씨 정도로 K단위로 상수로 줬다.
열전달 이론에 밝다면 좀 더 엄밀하게 Heat Flux 또는 시간이나 위치에 따라 변화하는 온도 조건 등으로 더 상세하게 줄 수 있을 것이다.

초기조건의 경우, 아주 살짝 초기 속도를 주었는데 Tutorial에 따르면 이렇게 함으로써 기하학적 대칭 형상을 살짝 깨뜨려주기 때문에 해석이 용이하게 해 주는 효과가 있다고 한다.

500프레임의 해석 시간은 약 25분 정도 걸린 것 같다.
매쉬도 적당히 성글기 때문에 꽤 빨랐다.

결과는 아래와 같다.

본 해석의 경우, 3가지 물리 현상이 고려된 것으로 볼 수 있다.

우선 Navier-Stokes Eq을 사용한 유체의 거동.
또 Heat Eq을 사용한 대류 열전달 현상.
마지막으로 Boussinesq 모델을 적용한 부력 작용이 그것이다.
부력의 경우, 온도에 따라 해당 물질의 밀도차이가 발생하므로 그것을 고려하여 Force를 발생시키는 것이다.

복사 열전달 현상은 옵션에서 뺐기 때문에 무시되었다.  이 경우에는 거의 의미가 없기 때문이다.

본 해석에서 도메인을 0.4x0.4m 사이즈로 주고, 가운데 열원을 80도씨로 0.2m 사이즈로 준 이유는 다름이 아니고, 요새 유행하는 염가형 자작 3D 프린터들이 대부분 그정도 사이즈이기 때문이다.
3D 프린터 오픈소스 사이트인 RapRap에 가서 한 번 훑어보던 중, 재미있는 발명품이 하나 보이던데, 바로 Heat Plate가 그것이다.  적층된 플라스틱 재료가 좀 더 잘 접착되면서 쌓아올라갈 수 있도록 하기 위해서는 어느정도 녹진녹진하게(?) 열을 좀 주는 것이 좋다고 한다.

그래서 상용 제품들 보면 여러가지 수단으로 바닥에 다양한 열원과 Thermal Sensor를 달아서 온도제어까지 실시하고 있다.
RapRap에서 제시된 여러 열원들 중에서, PCB 상의 회로 패턴을 아주 가늘고 길게 형성시켜 Copper의 전기저항을 최대한 끌어내어 열이 발생하도록 한 것이 있었는데, 아무래도 동박 재질이다 보니 저항값이 낮은지라 잘 될까 싶었는데, 실험한 것을 보니 대충 60도씨 정도까지 온도가 올라가더라는 것이다.  고정관념에 구애되지 않은 독창성이 돋보였다.
게다가 그 아이디어는 여러 개인 제작자들의 손을 거쳐가면서 점진적으로 개량되어가고 있었다.

아무튼 이런 식의 Heat Plate가 3D 프린터의 챔버 안에 설치되어 적절히 온도제어가 이루어진다면, 문제는 그 챔버 안에서 열이 어떻게 확산되어가는가 하는 구체적인 정보가 있다면 더 도움이 될 것이다.

본 해석 모델은 매우 간단하게 예제를 따라한 것에 불과하지만, 이제 하나씩 조건을 덧붙여가면서 좀 더 유용한 결과를 도출해 낼 수 있지 않을까 한다.
예를 들어, Heat Plate 위에 올라앉은 프린트 물체도 넣어 보고, 그것의 초기 온도와 열전도 현상까지 연성시켜 보고, 또 그 위에서 움직이는 토출기(Extruder)의 형상 및 발열도 함께 고려해 볼 수 있을 것이고, 아울러 열원이 되는 PCB Plate의 전기저항에 의한 발열까지 모사해서 한꺼번예 연성시켜볼 수 있을 것이다.  물론 3D 상태로 해서 보면 상당히 유용한 정보를 얻어낼 수 있을 것으로 예상된다.

이 결과를 통해 최적의 프린트 속도가 얼마나 되어야 할지를 반복실험 보다는 이론적인 결과에서 도출된 예상치를 확립할 수 있지 않나 한다.  또 내부 대류에 의한 열 분포 양상을 알 수 있으므로, 적층되는 층수에 따라 프린트 속도를 변화시켜 간다던가 하는 등의 소프트웨어적 발전도 기대할 수 있을 것이다.

충분히 조사해 본 것은 아니지만, 이렇게 3D 프린터를 위해서 공개적으로 CAE를 한 예가 거의 없는 것 같다.  아무리 검색해도 안 나온다.



마지막으로, 이상 검토한 조건을 기반으로, 충분히 더 잘게 쪼개어서 해상도를 높여서 돌려본 결과도 올려둔다..