Ex01

윈도우 PC에 해석 소프트웨어들 설치

Gist로 작성된 문서를 포함시켜 보았다.

기계공학용 가벼운 리눅스 배포판

CrunchBang 리눅스가 원래 Ubuntu + Openbox 조합으로 매우 가볍고 빠르면서 편리성도 갖추는데 성공한 이후, 이에 영향을 받은 여러 배포판들이 나온 것 같다.
2009년을 마지막으로 CrunchBang 프로젝트가 중단되었기 때문에 이와 같은 컨셉의 좋은 배포판에 대한 필요성이 더욱 증가하였다.

Debian + OpenBox 조합형 리눅스 배포판

• BunsenLabs : https://www.bunsenlabs.org/

BunsenLabs라는 프로젝트는 CrunchBang의 직계 후계 프로젝트로 계속 이어지고 있다.
간단히 말해 그냥 CrunchBang의 업데이트가 계속된다고 봐도 될 듯 하다.
다만 CrunchBang의 마지막 버전은 Ubuntu 대신 Debian 기반으로 전환하였기 때문에, 이것 역시 Debian 기반이다. 이점은 장점이자 약점이기도 하다. 사용목적에 따라 장단점은 달라질 것이다.
일단 설치 자체는 매우 안정적이고 신뢰도가 높다. 설치후 스크립트에 의한 환경 구성 과정이 너무 세심하게 잘 되어 있어서, APM 웹서버 구축하는 것조차도 선택에 따라 설치시에 자동으로 다 이루어진다.
개인용 웹서버와 간단한 데스크탑 정도 목적이라면 Debian 기반 시스템이 좋은 선택일 것이다.
그러나 워크스테이션으로 사용하기에는 Debian 시스템의 보수성 때문에 불편한 점이 많다.
많은 공학용 어플리케이션들이 대부분 Ubuntu 배포판은 기본으로 지원하지만, Debian 배포판은 따로 제공하지 않는 경우가 많기 때문이다. 때문에 Debian에서 억지로 사용하려면 직접 어플리케이션을 빌드해서 사용해야 할 것인데, 이는 시스템 관리에 많은 시간과 노력이 낭비된다는 것을 의미한다.
또한, 하드웨어를 자동으로 잡아주는 호환성 부분이나, 기타 여러가지 편리성 측면에서도 Debian이 아무래도 Ubuntu보다 좀 떨어질 수 밖에 없다. 설치 직후 인식되지 않는 하드웨어 설정 등에 시간이 투입될 확률이 좀 있다.

Ubuntu + OpenBox 조합형 리눅스 배포판

• SalentOS : http://www.salentos.it/
• Madbox (32bit 버전만 제공) : http://madbox.tuxfamily.org/
• GoBang : http://gobanglinux.org
• Chromixium : http://chromixium.org

따라서, 편리성과 가벼움(=리소스 절약)을 동시에 달성하기 위해서는 Ubuntu + OpenBox 조합이 최선의 선택이라고 생각된다. 설치하자마자 별다른 설정 없이 그대로 사용을 시작해도 충분한 수준이라고 생각된다.
위에 조사해 본 4가지 배포판 중에서, SalentOS 및 MadBox는 사용을 안해봤다.
MadBox는 32비트 바이너리만 제공하므로 공학해석을 위해 필요한 64비트 버전이 없기 때문에 제외해야 할 것이다.
GoBang은 시험삼아 설치해 봤는데, 설치 요령은 BunsenLabs와 크게 다르지 않으며 매우 가볍기는 하나, 안정성이 좀 떨어진다는 느낌이 든다. VirtualBox 가상머신에서 게스트 확장 설치에 실패했다. 다운되는건 아닌데 내부적으로 충돌이 좀 있어 보인다.
Chromixium은 OpenBox 기반이기는 하나, Tint2 대신 Plank Dock을 기본적으로 사용하기 때문에 훨씬 더 예쁘다. 계속 활발하게 업데이트가 이루어지고 있어서 이제는 매우 안정성이 높아졌다.
편리성을 포기하지 않으면서 동시에 메모리 절약 등 가벼움을 원한다면 현재로서는 가장 좋은 선택이라고 생각된다.
물론 가볍다고는 하지만 GoBang 수준으로 터무니없이 가볍지는 않다. GoBang이나 BunsenLabs의 경우에는 환경설정 대략 해 주고 부팅하면 메모리 사용량이 300~400MB 수준이다. 마우스 클릭시 반응속도는 뭐 말할 필요도 없다. Chromixium은 필요한 환경설정을 대충 해 주고 부팅하면 메모리 사용량이 600~700MB 수준이다. 구글 크롬북을 흉내낸 데스크탑 디자인 때문에 나름대로 산뜻한 느낌을 준다.
내가 이런 가벼운 배포판를 찾는 이유는, 가용한 PC의 메모리가 4~8GB에 불과한 경우가 많기 때문이다. 16GB~32GB 정도의 메모리를 가지고 있다면 고민할 필요가 크게 없다고 생각된다. 그때는 그냥 일반적인 Ubuntu Unity, Mint 같은 것을 사용하는 것이 더 편리하기 때문에 시간낭비를 크게 막아줄 것이다.
유한요소 해석을 할때, 대략 20만개 수준의 엘리먼트를 해석하면 4GB 메모리로도 충분히 커버 가능하다. 그런데 100만개 수준의 엘리먼트를 해석할 경우에는 거의 100% 실패한다. 메모리가 기하급수적으로 소모되기 때문에 메모리 용량은 나름대로 중요한 문제가 된다.
OS에서 기본적으로 점유하는 메모리를 아낄 수 있다면 좀 더 나을 것이고, 또 OS의 환경구성 복잡도가 낮아 단순하다면 더욱 해석에 유리할 것이라고 나름 생각한다.
또한, 리눅스 구성 자체에 재미를 느끼는 오타쿠(?)가 아닌 평범한 기계공학 엔지니어의 경우에는... 리눅스 사용 자체에 거부감을 느끼는 경우가 많다. 컴퓨팅 자체 보다는 설계 해석 행위 자체가 주관심사이기 때문이다. 컴퓨터는 그냥 도구라는 생각이 강하기 때문에, OS 구성하고 셋팅하는데 시간낭비가 없는 것이 좋겠다.
기계공학 엔지니어가 사용할 어플리케이션은 기껏해야 다음 정도다.

1. 유한요소해석 솔버
2. C,C++,Python 등의 간단한 개발환경
3. 3D CAD 관련 어플리케이션들
4. 간단한 웹서버 운용
5. SAMBA, SSH, FTP 정도의 유틸리티들
6. 웹브라우저

무거운 배포판은 필요가 없다.

Mechanical Contact Analysis with Elmer

Non-linear Contact 해석은 원래 Elmer에서 수행하기에 그다지 적합하지 않았다.

이유는 Elmer에서 지원되는 관련 기능은 가상의 접촉선을 함수나 직선으로 표현한 후,
그 이상 변형되지 않도록 Limit을 걸어서 변형량을 제한하는 기능 밖에 없었기 때문이다.

MultiBody 상황에서, 상호 접촉되는 영역을 찾아낸 후, 간섭 깊이를 노드별로 계산해서 그에 따라 Load를 걸어주는 전형적인 Hertz 접촉 해석 알고리즘이 없었던 것이다.

때문에 변형량을 인위적으로 제한하는 변칙적인 방법을 쓸 수 밖에 없었는데...

그런데 최근에 Elmer 개발팀이 왠일로 접촉 해석에 대해서 관심이 생겼는지, 그와 관련된 기능을 추가한 것 같다.
광고도 없고 공지도 없고 설명서 따위도 없이 그냥 Github에 슬그머니 업데이트 해 둔 것을 보았다.

https://github.com/ElmerCSC/elmerfem/tree/contact/fem/tests/ContactPatch3D

내역을 보니 테스트 완료 업데이트 된 것이 딱 일주일 되었다.
따끈따끈하다.

아무튼 닥치고 내려받아서 테스트 해 봤는데 잘된다.

case.sif 파일의 관련 문구들의 의미는 설명서가 전혀 없기 때문에 전부 파악하기는 아직 무리다.

다만 원래 있었던 "Mortar BC" 및 새로 추가된 것으로 보이는 "Contact BCs" 경계조건을 혼합해서 사용하고 있는 듯 하다.

여기서 "Mortar BC" 라는 것은, 자기장 해석을 위해 Elmer팀이 만들어 넣었던 기능인데 2개의 MultiBody가 접촉되어 있을 경우, 한쪽에서 발생하는 물리량을 접촉된 다른 쪽으로 전달할 수 있게 해 주도록 하는 것으로 파악된다.  (나도 아직은 충분히 이론 학습과 연습을 한 상태는 아니라서 정확한 설명인지는 모르겠다.)

아무튼 그래서 Mortar BC를 이용해서 두 개의 바디가 서로 상대운동하면서 위치가 변하는 상황에서 전자기 플럭스라던가, 온도 같은 것들을 각 스텝별로 전달하는 것이 가능한 것 같다.
이를 응용해서 이번에는 Load를 전달하는 것으로 생각된다.

Contact BC의 경우는, 구체적으로 설명된 자료가 없어서 추측만 해 볼 수 밖에 없는데....
변위에 의해 발생하는 침투량을 계산해서 발생하는 Load를 계산해 내는 전형적인 접촉 해석용 기능 아닐까 싶다.  물론 접촉되는 영역을 디텍트 해 내는 기능도 있을 것으로 추측되는데 아직 확인이 안된다.

또 경계조건의 특성은 "Initialize Dirichlet Conditions = False" 으로 선언된 것으로 보아, 기본 셋팅된 경계조건과는 좀 다른 것 같다.
( 참고 : https://en.wikipedia.org/wiki/Dirichlet_boundary_condition )


아무튼 정확한 Notation이나 구조 파악은 현재로서는 어려우나, 일단 제시된 예제를 사용해서 해 보았고, 그것을 조금 변형해서 Scanning도 해 보았다.

나중에 시간이 나면 일반적으로 사용할 수 있는 수준인지 검증해 볼 수 있을 듯 하다.
임의로 변형한 솔버 인풋 파일은 아래와 같다.
모델 데이타 자체는 원래의 것을 그대로 사용하였다.
(아쉬운 점은, 모델링된 데이타가 1/4만 그린 후 대칭성을 이용해서 조합하는 방식으로 구성되어 있는지라, 접촉 해석 관련 부분만 따로 떼어서 보기가 조금 헷갈린다.)

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

$fileid="c"

Simulation
  Max Output Level = 20
  Coordinate System = Cartesian 3D
  Coordinate Mapping(3) = 1 2 3
  !Simulation Type = Steady State
  Simulation Type = Scanning
  Timestep intervals = 200
  Timestep sizes = 1

  Steady State Min Iterations = 1
  Steady State Max Iterations = 1

  Post File = case_$fileid$.vtu
  Save Geometry Ids = Logical True

  Mesh Levels = 2

! The ElasticSolver does not really like the Dirichlet conditions at the start 
! of the nonlinear iteration. 
  Initialize Dirichlet Conditions = False
End

!Constants
!  Gravity(4) = 0 -1 0 9.82
!End

Equation 1
  Name = "Deformation"
  Active Solvers(1) = 1
  Plane Strain = Logical True
End

Body 1
  Name = "Lower block"
  Target Bodies(1) = 55
  Equation = 1
  Material = 1
End

Body 2
  Name = "Upper block"
  Target Bodies(1) = 56
  Equation = 1
  Material = 1
End

Material 1
  Name = "Ideal"
  Youngs modulus = 90.0
  Density = 10.0
  Poisson ratio = 0.25
End

Solver 1
  Equation = "Nonlinear elasticity"
  Procedure = "ElasticSolve" "ElasticSolver"
  Variable = -dofs 3 Displacement

  Nonlinear System Convergence Tolerance = 1.0e-5
  Nonlinear System Max Iterations = 10
  Nonlinear System Relaxation Factor = 1.0

  Linear System Solver = "Iterative"
  Linear System Iterative Method = "BiCGStab"
  Linear System Abort Not Converged = True
  Linear System Preconditioning = "ILU2"
  Linear System Residual Output = 100
  Linear System Max Iterations = 5000
  BiCGStabl Polynomial Degree = 4
  
  Linear System Convergence Tolerance = 1.0e-10

  Apply Contact BCs = Logical True
! Save Contact = Logical True

! Restore the linear solution
! Elasticity Solver Linear = Logical True

  Calculate Stresses = Logical True
! Optimize Bandwidth = False

  Displace Mesh = Logical True
! Calculate Boundary Weights = Logical True

! Do not include constraints when analyzing the convergence and norm of a solution
  Nonlinear System Convergence Without Constraints = Logical True
End

Solver 2
  Exec Solver = never
  Equation = "SaveLine"
  Procedure = "SaveData" "SaveLine"
  Filename = f_$fileid$.dat
End

Boundary Condition 1
  Name = "Support"
  Target Boundaries(1) = 57
  Displacement 3 = Real 0.0
End

Boundary Condition 2
  Name = "Lower surface of upper block"
  Target Boundaries(1) = 58

  Mortar BC = Integer 3
  Mortar BC Nonlinear = Logical True
  Contact Depth Offset Initial = Real 1.0e-3
  !Contact Active Set Minimum = Integer 1
  !Contact No-Slip = Logical True

! Create a strong projector for the line setting y-coordinate to zero
  Flat Projector = Logical True

! a) Use weak projector or not
  Galerkin Projector = Logical False

! b) Use more tailored projector able to do accurate integration
  Level Projector = Logical True
  Level Projector Generic = True
End

Boundary Condition 3
  Name = "Upper surface of lower block"
  Target Boundaries(1) = 59
End

Boundary Condition 4
  Name = "Pressure load the upper surface of upper block"
  Target Boundaries(1) = 60
  !Normal Surface Traction = -1.0
  Normal Surface Traction = Variable Time
    Real MATC "-1*sin(2*pi*tx/200)-1"
End

Boundary Condition 5
  Name = "Symmetry y-z"
  Target Boundaries(1) = 61
  Displacement 1 = 0.0
End

Boundary Condition 6
  Name = "Symmetry x-z"
  Target Boundaries(1) = 62
  Displacement 2 = 0.0
End


Solver 1 :: Reference Norm = 0.66961642E-02
Solver 1 :: Reference Norm Tolerance = 1.0e-6

아무튼 이 기능이 잘 발전되고, 나중에 GUI에까지 반영이 된다면 Elmer가 Mechanical CAE 하는데 훨씬 더 활용성이 높아질 것 같다.
OpenSource CAE 툴 중에서 가장 접근하기 쉬운게 바로 Elmer 이기 때문에, 매우 손쉽게 고수준 해석에 도전할 수 있는 환경이 되는 것이다.

Transient Linear Elastic with Elmer

Elmer에서 구조해석을 위한 기본적인 해석기는 실질적으로 Linear Elastic Eq. 밖에 없다. 

얼핏 보면 사용할 수 있는 범위가 굉장히 제한적으로 오해받을 수 있으나, 자꾸 살펴보다 보니 상당히 응용범위가 넓은, 잘 만들어진 해석기라는 생각이 든다.

기본적으로 이 해석기는 Implicit Time Integrator라고 생각되는데, 그렇다고 해서 Explicit Time Integration을 통해 얻을 수 있는 결과를 얻을 수 없는 것은 아니다.

 Initial Condition 또는 Boundary Condition을 Time Varient 하게 잘 주면 올바른 응답을 보여주는 것 같다. (물론 실제 물리현상과 매칭되도록 조건을 잡아주는 것은 계속 반복 검증해 가면서 맞춰가야 할 것이다. 본 글에서는 일단 기능이 된다는 점만 확인한다.)

물론 좋은 상업용 Explicit 해석기 또는, 오픈소스에서 Explicit 전용으로 개발된 "Impact"같은 툴을 잘 쓰는 것 보다는 결과를 얻어내는데 들어가는 시간이 더 많이 걸리는 것은 당연하다.

하지만 익숙한 툴을 이용하여 원하는 결과를 도출하는데 충분할 수도 있다.

(사실 "Impact" 같은 프로그램은 Java로 만들어져 있어서 대규모의 Fine Mesh 해석 같은 것은 아마 꿈도 못 꾸는 것 같다.  예제를 보면 전부 다 굉장히 매쉬 앨리먼트 개수가 적다.)

아무튼 일단 기능 검증을 위해 간단한 컨틸레버를 그렸다.




길이 100mm, 폭 10mm, 두께 5mm이고, Transfinite 옵션을 줘서 균질하게 매슁이 되도록 했다.  물론 2nd Order로 쪼갰으며, 해석시간 단축을 위해 거칠게 쪼갰다.

그리고 이것을 ElmerGUI로 아래와 같이 읽어들인다.




GUI 상에서 해석 조건을 준다.
해석 입력 파일은 아래와 같다.

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 = 3000
  Timestep Sizes = 0.0001
  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) = 1
  Name = "Body 1"
  Equation = 1
  Material = 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 2
  Equation = Linear elasticity
  Procedure = "StressSolve" "StressSolver"
  Variable = -dofs 3 Displacement
  Calculate Stresses = True
  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

Equation 1
  Name = "Equation 1"
  Calculate Stresses = True
  Active Solvers(2) = 1 2
End

Material 1
  Name = "Polycarbonate (generic)"
  Heat Conductivity = 0.205
  Youngs modulus = 2.2e9
  Mesh Poisson ratio = 0.37
  Heat Capacity = 1250.0
  Density = 1220.0
  Poisson ratio = 0.37
  Heat expansion Coefficient = 67.0e-6
End

Initial Condition 1
  Name = "InitialCondition 1"
  Displacement 3 = -0.005
End

Boundary Condition 1
  Target Boundaries(1) = 5 
  Name = "Fix"
  Displacement 3 = 0
  Displacement 2 = 0
  Displacement 1 = 0
End

이제 해석을 돌린다.  1000개의 샷을 아주 짧은 시간 간격으로 찍어내기 때문에 시간이 꽤 걸린다.

총 시간은 0.3초의 응답을 보는 것이다.

그래도 수렴이 잘 되는 것 같다.

일단 이 해석의 결과를 보면 아래 동영상으로 나온다.

처음에 가진될 때 보이는 2차 이상의 고주파 영역은 금방 감쇄되는 것이 보인다.

참고로 사용된 물성치에는 감쇄계수가 전혀 없다.  무감쇄 조건이다.

그럼에도 불구하고 1차 모드 진동만 남고 다른 진동은 다 사라진다.  

1차 모드 진동은 무한대로 감쇄하지 않고 진동이 계속된다.

이제 물성치에 감쇄조건을 추가해 본다.

Damping Coefficient 를 어느정도로 주는 것이 적절한지 전혀 감이 잡히지 않는데, 일단 아무 숫자나 넣어보고 해 봤는데 간에 기별도 안 가길래 숫자를 아주 크게 확 키워 보았다.

Material 카테고리에 아래의 문구가 들어가도록 해 준다.  (GUI에서도 지원됨)

  Damping = 1e6

그리고 더 키워 본다.

  Damping = 1e7

이상의 결과로, Damping 과도 응답(Transient)을 보는 것이 가능하며 이것은 Explicit Time Integration 하는 것과 같은 좋은 기능이라는 사실을 알았다.

대신 해석시간은 매번 전부 다 계산하는 Implicit 방식이기 때문에 꽤 소요된다는 단점이 있다.

결론적으로 Static, Eigenmodal Dynamic, Harmonical Dynamic 뿐만 아니라, Transient Dynamic 결과를 좀 효율이 떨어지더라도 얻어낼 수 있다는 것.

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를 한 예가 거의 없는 것 같다.  아무리 검색해도 안 나온다.



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

Karman Vortex CFD with Elmer

관로 내부 유동의 전형적인 예제 중에서 가장 기본적인 Karman Vortex(칼만 와류)의 형성 조건을 파악해 보기 위해 Elmer를 이용한 CFD를 실시해 보았다.

1. 우선 기본적인 와류가 발생하지 않는 층류 조건을 먼저 테스트.

모델링은 GMSH로 했다.
전체 관로 길이는 1m, 폭은 0.2m, 중간 공의 지름은 0.02m.
매쉬는 일반적인 삼각형으로 형성했는데, 조금 성글다는 생각이 든다.

정상상태 도달시의 속도 분포이다.  유체는 공기(Air)로 했는데, 이렇게 하니까 입구속도를 0.7m/s까지 높여도 칼만 와류가 형성되지 않는다.  관로가 너무 좁아서 그렇다기 보다는, 공에 관련한 레이놀즈계수가 적절히 맞춰지지 않았기 때문이라고 생각된다.

속도를 그보다 더 높여 보니까 해석시 수렴에 실패한다.  해석을 억지로 시켜내려면 시간 분할을 현재의 25e-4s보다 훨씬 더 쪼개서 해 보던가 할 필요가 있지 않나 싶다.
그리고, 매쉬가 너무 성글기 때문에 그림이 거칠게 나와서 예쁘지가 않다.

2. 위 결과를 참고하여 몇가지 조건을 수정하여 Karman Vortex 형성에 성공.

우선 모델을 조금 수정했다.
공의 크기는 D0.02m로 그대로 유지했지만, 관로 길이를 0.5m로 줄였고, 관로의 폭은 0.4m로 늘렸다.  벽면의 겅계층의 영향이 공 주변으로 미치지 않는 것이 좋겠다 싶어서다.
그리고, 그림을 예쁘게 나오게 하기 위해 매쉬 밀도를 크게 높였다.  해석 시간이 굉장히 오래 걸릴 것은 각오한다.

재료를 공기에서 물로 바꾸니까 칼만 와류 형성에 성공한다.  당연히 공기보다 물의 점성이 높아져서 레이놀즈 계수가 달라져서 그럴 것으로 생각된다.

200프레임을 해석하는데 62117.50s 즉 약 17.25h 정도 소요되었다.

(물론 오래된 Core2Duo CPU를 사용하는 노트북에서 Single Core로 해석해서 그렇다.)

2D 해석이 이정도인데 3D 해석은 엄두가 안 난다.  

아무튼 가까스로 와류 형성은 일단 해 보았다.
아래에 본 해석에 적용한 해석 입력 파일을 기록해 둔다.

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 = 200
  Timestep Sizes = 25e-4
  Solver Input File = trans4.sif
  Post File = trans4.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) = 1
  Name = "Body 1"
  Equation = 1
  Material = 1
End

Solver 1
  Equation = Result Output
  Single Precision = True
  Procedure = "ResultOutputSolve" "ResultOutputSolver"
  Save Geometry Ids = True
  Output File Name = trans4
  Output Format = Vtu
  Exec Solver = After Timestep
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 = False
  Steady State Convergence Tolerance = 1.0e-5
  Nonlinear System Convergence Tolerance = 1.0e-7
  Nonlinear System Max Iterations = 20
  Nonlinear System Newton After Iterations = 15
  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

Equation 1
  Name = "Equation 1"
  Active Solvers(2) = 1 2
End

Material 1
  Name = "Water (room temperature)"
  Heat Conductivity = 0.58
  Heat Capacity = 4183.0
  Density = 998.3
  Relative Permittivity = 80.1
  Viscosity = 1.002e-3
  Sound speed = 1497.0
  Heat expansion Coefficient = 0.207e-3
End

Boundary Condition 1
  Target Boundaries(1) = 1 
  Name = "Inlet"
  Velocity 3 = 0
  Velocity 1 = 0.5
  Velocity 2 = 0
End

Boundary Condition 2
  Target Boundaries(1) = 2 
  Name = "Outlet"
  External Pressure = 0
End

Boundary Condition 3
  Target Boundaries(4) = 3 4 5 6 
  Name = "Wall"
  Noslip wall BC = True
End

여기서, 몇가지 더 손을 봐서 좀 더 긴 영상을 만들었다.
변경된 부분은 Ball의 직경은 0.4m로 2배 더 키우고, 매쉬는 더욱 성글게 해서 해상도를 희생하는 대신 해석시간을 단축시켰고, 해석 프레임을 1000개로 늘렸다.

칼만 와류를 공학적으로 응용하는 예를 들어 보자면, 일단 유량계 중에서 칼만 와류 원리를 응용한 제품들이 있다.  유량계 내부 관로 안에서 와류를 형성시킨 후, 와류에 의해 발생하는 초음파를 압전소자로 계측해서 그 주파수를 가지고 유량을 환산하는 방식이다.
초음파로 신호가 나오는 이유는 당연히 유량계가 쪼그맣기 때문이다.

이 외에, 최근의 사례를 보면 인디고고 크라우드 펀딩 사이트에 보니깐 칼만 와류를 이용하여 풍력발전기를 발명한 사례가 있다.  MIT 친구들 같은데, 5만불 목표금액을 채우는데 성공했나 보다.

https://www.indiegogo.com/projects/vortex-bladeless-a-wind-generator-without-blades--3#/story

이름하여 Vortex Bladeless.
별 건 아니고, 원형 단면 기둥을 하나 세워서 와류가 형성되면 기둥이 흔들거릴 테니깐, 그걸로 유도기전력을 발생시켜 전력을 생산한다는 거다.
당연히 효율은 좀 떨어지겠지만, 기존의 블레이드 방식의 경우 고가의 고정밀 동력전달기구 및 회전형 발전기가 들어가기 때문에 가격이 엄청 높아진다는 점을 공략했다.
고속 회전기구가 안 들어가므로 매우 싸게 만들수 있기 때문이다.
재료도 가볍고 싼 수지 계열로 했다.
또 회전 날개가 없으니깐 새들이 죽어나가는 일도 없다는 '명분'도 만들었다.
칼만 와류 모르는 기계공학도는 없겠지만, 이런 식으로 명분 만들어서 사업화 하는건 코페르니쿠스적 발상이라고 생각된다.  또 MIT라는 브랜드를 등에 업었으니 더 쉽게 되겠지.

이 제품의 디자인을 보면, 위로 갈수록 단면 지름이 커지도록 해 놨는데, 이는 지표면에서 위로 올라갈수록 점성에 의해 바람의 속도가 증가하게 되므로, 이를 고려하여 구배값을 결정한 것 같다.

또한, 세워진 원통 막대를 지지하는 카본 목 부위의 탄성 등을 감안하여, 막대의 고유진동수를 칼만 와류에 의해 가진되는 진동수와 잘 맞춰주도록 설계함으로써, 진동에 의해 발생하는 운동에너지를 극대화하려고 노력했을 것이다.

그리고, 유량계의 경우 사이즈가 작아서 진동이 초음파대로 나오도록 설계되어 있지만, 이 녀석의 경우에는 충분히 직경이 크기 때문에 인디고고에서 이 사람들이 밝힌 바에 따르면 음파는 20Hz 미만의 낮은 주파수대로 나온다고 한다.  따라서 새나 기타 동물, 인간에게 초음파에 의한 악영향이 없을 거라고 밝히고 있다.

FGPG_V11 ::: Introduce of Fine Involute Gear Tooth Profile Generator

FGPG_V11


새로운 인볼류트 기어 설계용 프로그램을 작성해 보았다.
기존의 GPG_V08을 대체한다는 의미로 FGPG 라고 이름을 붙였다.
Fine Gear Profile Generator 의 줄임말이다.
아직 미완성본이지만, 일단 필요한 기능들은 대부분 어설프게나마 다 구현되었으므로 기록 차원에서 올려본다.

특징은 대략 다음과 같다.

1. Julia Language로 작성되었다.
2. 아직 독립 실행파일(Standalone Executable)은 빌드하지 않았다.
3. Julia를 이용하여 Python MatPlotLib의 플랏팅 능력을 활용한다.
4. Gmsh 및 Elmer를 이용하여 자동적으로 FEM 해석을 실시해 주는 Batch 파일을 만들어준다.  윈도우용 뿐만 아니라, 리눅스용 Shell Script도 만들어지므로 운영체제에 상관없이 사용 가능하다.
5. 입력하는 파라미터들은 별도의 파일로 작성해 두면, 그것을 읽어서 계산해 낸다.
6. GUI는 없다.
7. AutoCAD 또는 DraftSight에서 기어 형상을 자동으로 그려줄 수 있도록 Script 파일이 자동으로 생성된다.  여기서 dxf 파일로 저장한 후, 3D CAD 에서 모델링 할 수 있다.


우선 프로그램의 소스코드는 다음과 같다.
파일 이름은 FGPG_V11.jl 이다.  다른 이름을 붙여도 상관없다.

#########################################################################
# FGPG (Fine_Gear_Profile_Generator) V11
# 20150602
# by Dongho Kim  from Korea
# dymaxion.kim@gmail.com
# http://dymaxionkim.blogspot.kr/
#########################################################################
## Ref.
# http://tro.kr/29
# http://en.wikipedia.org/wiki/Coordinate_rotations_and_reflections
# http://dip28p.web.fc2.com/
#########################################################################
## History
# 20141114 V0.44 : SageMath Version
# 20150519 V0.1 : Trying to Conversion for Julia, Trying Interact Macro
# 20150519 V0.2 : Studying about Lack
# 20150527 V01 : Plotting Whole Gear at first
# 20150527 V02 : Trying Interact Macro
# 20150527 V03 : Debug V01
# 20150527 V04 : Adding Annotations, Saving Profile Data File, ...
# 20150527 V05 : Adding More Error Checks, Making geo File
# 20150528 V06 : Debug Making geo File
# 20150528 V07 : Align vertically Gear, Add Shaft
# 20150528 V08 : Modifying Sequences, Debug Making geo File, Make new Figure, Make svg Pictures, Save scr File, Modify Annotations
# 20150530 V09 : Modifying Shaft, Internal-External Condition, geo File
# 20150531 V10 : Adding version Variable
# 20150602 V11 : Making case.sif, case.bat, case.sh files
#########################################################################
## Todo
# Write Theory Doc
# Write User Manual
# Write Some Actual Examples
# Tip Relief to Modifying Tooth Profile
# Making Spec Table
# Making case.sif
# Profile Repairing in Undercut Condition
# Technical Error Check & Anouncing
# Change Variable Names for easy code reading
# Re-code in Functional Criteria (?)
# Wndows Batch File
# Bash Script
# Change Result File Name based in Time
# Build to Standalone Executable

# On JUNO envirenment
cd(dirname(@__FILE__))

using PyPlot
PyPlot.svg(true)

############################
# Parameters
version = "v11"
m = 1  # Module
z = 20  # Teeth Number
alpha_0_deg = 20  # Pressure Angle [Deg]
x = 0.2  # Offset Factor
b = 0  # Backlash Factor
a = 1.0  # Addendum of Tooth
d = 1.2  # Dedendum of Tooth
c = 0.3  # Radius Factor of Edge Round of Hob
e = 0.2  # Radius Factor of Edge Round of Tooth
# Center of Gear
x_0 = 0
y_0 = 0
# Segmentation Points Numbers for each curves
seg_circle = 360
seg_involute = 50
seg_edge_r = 4
seg_root_r = 20
seg_center = 2
seg_outter = 4
seg_root = 4


############################
# Input Parameters from input4fgpg.csv file
input4fgpg = readcsv("input4fgpg.csv")
input4fgpg[:,3] = float64(input4fgpg[:,3])
m = input4fgpg[1,3]
z = int64(input4fgpg[2,3])
alpha_0_deg = input4fgpg[3,3]
alpha_0 = alpha_0_deg*2*pi/360
x = input4fgpg[4,3]
b = input4fgpg[5,3]
a = input4fgpg[6,3]
d = input4fgpg[7,3]
c = input4fgpg[8,3]
e = input4fgpg[9,3]
x_0 = input4fgpg[10,3]
y_0 = input4fgpg[11,3]
seg_circle = int64(input4fgpg[12,3])
seg_involute = int64(input4fgpg[13,3])
seg_edge_r = int64(input4fgpg[14,3])
seg_root_r = int64(input4fgpg[15,3])
seg_center = int64(input4fgpg[16,3])
seg_outter = int64(input4fgpg[17,3])
seg_root = int64(input4fgpg[18,3])


############################
# Involute Curve
# alpha_m = Center Line's Slope [Rad]
alpha_m = pi/z
# alpha_is = Start Angle for Involute Curve
alpha_is = alpha_0 + pi/(2*z) + b/(z*cos(alpha_0)) - (1+2*x/z)*sin(alpha_0)/cos(alpha_0)
# theta_is = Minimum Range of Parameter to Draw Involute Curve
theta_is = sin(alpha_0)/cos(alpha_0) + 2*(c*(1-sin(alpha_0))+x-d)/(z*cos(alpha_0)*sin(alpha_0))
# theta_ie = Maximum Range of Parameter to Draw Involute Curve
theta_ie = 2*e/(z*cos(alpha_0)) + sqrt( ((z+2*(x+a-e))/(z*cos(alpha_0)))^2 - 1 )
# Condition of theta_ie
if alpha_malpha_m && alpha_m>alpha_is+theta_ie-atan(theta_ie)
    e = (z/2)*cos(alpha_0)*( theta_ie - sqrt( (1/cos(alpha_is+theta_ie-alpha_m))^2-1 ) )
end
# x_e, y_e = Location of Tooth's End Point
x_e = x_0 + m*((z/2)+x+a)*cos(alpha_e)
y_e = y_0 + m*((z/2)+x+a)*sin(alpha_e)
# x_e0, y_e0 = Location of Edge Round Center
x_e0 = m*(z/2+x+a-e)*cos(alpha_e) + x_0
y_e0 = m*(z/2+x+a-e)*sin(alpha_e) + y_0
# Parameter Range of Edge Round
theta3_min = atan((Y1[length(Y1)]-y_e0)/(X1[length(X1)]-x_e0))
theta3_max = atan((y_e-y_e0)/(x_e-x_e0))
THETA3 = linspace(theta3_min,theta3_max,seg_edge_r)
X3 = m*e*cos(THETA3)+x_e0
Y3 = m*e*sin(THETA3)+y_e0
#plot(X3,Y3,color="green",linestyle="-")


############################
# Root Round Curve of Tooth
# Condition Check
# alpha_ts = Start Angle of Root Round Curve
# THETA_s = Substitution Variable to plot Root Round Curve
alpha_ts = (2*(c*(1-sin(alpha_0))-d)*sin(alpha_0)+b)/(z*cos(alpha_0)) - 2*c*cos(alpha_0)/z + pi/(2*z)
theta_te = 2*c*cos(alpha_0)/z - 2*(d-x-c*(1-sin(alpha_0)))*cos(alpha_0)/(z*sin(alpha_0))
THETA_t = linspace(0,theta_te,seg_root_r)
if c!=0 && (d-x-c)==0
    # mc를 반지름으로 하는 원호를 그려서 대체하게 됨
    THETA_s = (pi/2)*ones(length(THETA_t))
    elseif c==0 && (d-x-c)==0
    # 루트커브 작도 생략하고, (x_is,y_is) 점을 대칭으로 이뿌리호를 그려서 대체하게 됨
    elseif (d-x-c)!=0
    THETA_s = atan((m*z*THETA_t/2)/(m*d-m*x-m*c))
end
X_t = x_0 + m*( (z/2+x-d+c)*cos(THETA_t+alpha_ts) + (z/2)*THETA_t.*sin(THETA_t+alpha_ts) - c*cos(THETA_s+THETA_t+alpha_ts) )
Y_t = y_0 + m*( (z/2+x-d+c)*sin(THETA_t+alpha_ts) - (z/2)*THETA_t.*cos(THETA_t+alpha_ts) - c*sin(THETA_s+THETA_t+alpha_ts) )
#plot(X_t,Y_t,color="green",linestyle="-")


############################
# Outter Arc
#alpha_em = 2*alpha_m-alpha_e
#x_em = x_0 + m*(z/2+x+a)*cos(alpha_em)
#y_em = y_0 + m*(z/2+x+a)*cos(alpha_em)
THETA6 = linspace(alpha_e,alpha_m,seg_outter)
X6 = m*(z/2+a+x)*cos(THETA6)+x_0
Y6 = m*(z/2+a+x)*sin(THETA6)+y_0
#plot(X6,Y6,color="red",linestyle="-")


############################
# Root Arc
#x_rm = x_0 + m*(z/2+x-d)*cos(alpha_ts)
#y_rm = y_0 - m*(z/2+x-d)*sin(alpha_ts)
THETA7 = linspace(0,alpha_ts,seg_root)
X7 = m*(z/2-d+x)*cos(THETA7)+x_0
Y7 = m*(z/2-d+x)*sin(THETA7)+y_0
#plot(X7,Y7,color="red",linestyle="-")


############################
# Combine Curves
# Not using Root Round Curve (X_t,Y_t)
if c==0 && (d-x-c)==0
    Xc = X7[1:length(X7)-1]
    Yc = Y7[1:length(Y7)-1]
else
    Xc = [ X7[1:length(X7)-1]; X_t[1:length(X_t)-1] ]
    Yc = [ Y7[1:length(Y7)-1]; Y_t[1:length(Y_t)-1] ]
end
# Not using Edge Round Curve (X3,Y3)
if e==0 || alpha_m==alpha_is+theta_ie-atan(theta_ie)
    Xc = [ Xc; X1[1:length(X1)-1] ]
    Yc = [ Yc; Y1[1:length(Y1)-1] ]
else
    Xc = [ Xc; X1[1:length(X1)-1]; X3[1:length(X3)-1] ]
    Yc = [ Yc; Y1[1:length(Y1)-1]; Y3[1:length(Y3)-1] ]
end
# Not using Outter Arc (X6,Y6)
if alpha_e == alpha_m || alpha_m==alpha_is+theta_ie-atan(theta_ie)
    Xc = Xc
    Yc = Yc
else
    Xc = [ Xc; X6[1:length(X6)-1] ]
    Yc = [ Yc; Y6[1:length(Y6)-1] ]
end
#Xc = [ X7[1:length(X7)-1]; X_t[1:length(X_t)-1]; X1[1:length(X1)-1]; X3[1:length(X3)-1]; X6[1:length(X6)-1]]
#Yc = [ Y7[1:length(Y7)-1]; Y_t[1:length(Y_t)-1]; Y1[1:length(Y1)-1]; Y3[1:length(Y3)-1]; Y6[1:length(Y6)-1]]
#plot(Xc,Yc,color="black",linestyle="-",linewidth=2)



############################
# Make Whole One Tooth
Xc2 = Xc[2:length(Xc)]
Yc2 = Yc[2:length(Yc)]
# Reflect Transform
Xcc = cos(2*alpha_m)*(Xc2-x_0) + sin(2*alpha_m)*(Yc2-y_0)
Ycc = sin(2*alpha_m)*(Xc2-x_0) - cos(2*alpha_m)*(Yc2-y_0)
# Location Transform to (x_0,y_0)
Xcc = Xcc + x_0
Ycc = Ycc + y_0
# Invert
Xcc = Xcc[length(Xcc):-1:1]
Ycc = Ycc[length(Ycc):-1:1]
# Combine
Xcc = [Xc; Xcc]
Ycc = [Yc; Ycc]
#plot(Xcc,Ycc,color="black",linestyle="-",linewidth=2)


############################
# Align to Top
align_angle = pi/2-pi/z
X_align = Xcc
Y_align = Ycc
# Rotate
X_align = cos(align_angle)*(Xcc-x_0) - sin(align_angle)*(Ycc-y_0)
Y_align = sin(align_angle)*(Xcc-x_0) + cos(align_angle)*(Ycc-y_0)
# Location Transform to (x_0,y_0)
X_align = X_align + x_0
Y_align = Y_align + y_0
#plot(X_align,Y_align,color="black",linestyle="-",linewidth=2)


############################
# Make Whole Gear
p_angle = 2*pi/z
Xccc = X_align
Yccc = Y_align
Xtemp = X_align[2:length(X_align)]
Ytemp = Y_align[2:length(Y_align)]
for i = 1:z-1
    # Rotate
    Xtemp = cos(p_angle*i)*(X_align[2:length(X_align)]-x_0) - sin(p_angle*i)*(Y_align[2:length(Y_align)]-y_0)
    Ytemp = sin(p_angle*i)*(X_align[2:length(X_align)]-x_0) + cos(p_angle*i)*(Y_align[2:length(Y_align)]-y_0)
    # Location Transform to (x_0,y_0)
    Xtemp = Xtemp + x_0
    Ytemp = Ytemp + y_0
    Xccc = [Xccc; Xtemp]
    Yccc = [Yccc; Ytemp]
end



############################
# Plot Figure
f=figure(figsize=(10,10))
grid("on")
title("FGPG (Fine_Gear_Profile_Generator) $version \n by Dymaxion.Kim")
axis("equal")


############################
# Base, Pitch, Offset, Outter, Root Circle
THETA0 = linspace(0.0,2*pi,seg_circle)
# Base Circle
base_dia = m*z*cos(alpha_0)
X_base = base_dia/2*sin(THETA0) + x_0
Y_base = base_dia/2*cos(THETA0) + y_0
plot(X_base, Y_base, color="cyan", linestyle="--")
# Pitch Circle
pitch_dia = m*z
X_pitch = pitch_dia/2*sin(THETA0) + x_0
Y_pitch = pitch_dia/2*cos(THETA0) + y_0
plot(X_pitch, Y_pitch, color="magenta", linestyle="--")
# Offset Circle
offset_dia = 2*m*(z/2+x)
X_offset = (offset_dia/2)*sin(THETA0) + x_0
Y_offset = (offset_dia/2)*cos(THETA0) + y_0
plot(X_offset, Y_offset, color="red", linestyle="--")
# Outter Circle
outter_dia = 2*m*(z/2+x+a)
X_out = (outter_dia/2)*sin(THETA0) + x_0
Y_out = (outter_dia/2)*cos(THETA0) + y_0
plot(X_out, Y_out, color="brown", linestyle="--")
# Root Circle
root_dia = 2*m*(z/2+x-d)
X_root = (root_dia/2)*sin(THETA0) + x_0
Y_root = (root_dia/2)*cos(THETA0) + y_0
plot(X_root, Y_root, color="grey", linestyle="--")


############################
# Center Line of Tooth
X_m = linspace(x_0, (m*(z/2+a+x))*cos(alpha_m)+x_0, seg_center)
Y_m = (X_m-x_0)*tan(alpha_m) + y_0
plot(X_m,Y_m,color="orange",linestyle="--")


############################
# Shaft
THETA8 = linspace(0,2*pi,2*z)
if ad  # Internal Gear
    d_shaft = m*(z+x+a)+m*(4*(a+d))
    elseif a==d  # No correct
    #d_shaft = m*(z+x-d)-m*(2*(a+d))
    d_shaft = 0
end
X_shaft = (d_shaft/2)*sin(THETA8) + x_0
Y_shaft = (d_shaft/2)*cos(THETA8) + y_0
plot(X_shaft, Y_shaft, color="black", linestyle="-", linewidth=2)


############################
# Plot Whole Gear
plot(Xccc,Yccc,color="black",linestyle="-",linewidth=2)
char_scale = 0.03*m*z
annotate(["Module = $m mm"],xy=(x_0,y_0+5*char_scale),ha="center",va="center",color="black",fontsize="x-small")
annotate(["Teeth Number = $z ea"],xy=(x_0,y_0+4*char_scale),ha="center",va="center",color="black",fontsize="x-small")
annotate(["Pressure Angle = $alpha_0_deg DEG"],xy=(x_0,y_0+3*char_scale),ha="center",va="center",color="black",fontsize="x-small")
annotate(["Offset Factor = $x"],xy=(x_0,y_0+2*char_scale),ha="center",va="center",color="black",fontsize="x-small")
annotate(["Backlash Factor = $b"],xy=(x_0,y_0+1*char_scale),ha="center",va="center",color="black",fontsize="x-small")
annotate(["Addendum Factor = $a"],xy=(x_0,y_0+0*char_scale),ha="center",va="center",color="black",fontsize="x-small")
annotate(["Dedendum Factor = $d"],xy=(x_0,y_0-1*char_scale),ha="center",va="center",color="black",fontsize="x-small")
annotate(["Radius Factor of Edge Round of Hob = $c"],xy=(x_0,y_0-2*char_scale),ha="center",va="center",color="black",fontsize="x-small")
annotate(["Radius Factor of Edge Round of Tooth = $e"],xy=(x_0,y_0-3*char_scale),ha="center",va="center",color="black",fontsize="x-small")
annotate(["Center Position [mm] = ($x_0,$y_0)"],xy=(x_0,y_0-4*char_scale),ha="center",va="center",color="black",fontsize="x-small")
legend(["Base Circle Dia = $base_dia", "Pitch Circle Dia = $pitch_dia", "Offset Circle Dia = $offset_dia",
    "Outter Circle Dia = $outter_dia", "Root Circle Dia = $root_dia", "Center Line"], loc="lower left", fontsize="x-small")
savefig("case01.svg")





############################
# Figuer 2
f2=figure(figsize=(10,10))
grid("on")
title("FGPG (Fine_Gear_Profile_Generator) $version \n by Dymaxion.Kim")
axis("equal")


############################
# Base, Pitch, Offset, Outter, Root Circle 2
THETA9 = linspace(atan(X_align[length(X_align)]/Y_align[length(Y_align)]),
                  atan(X_align[1]/Y_align[1]), int(seg_circle/10))
# Base Circle
X_base9 = base_dia/2*sin(THETA9) + x_0
Y_base9 = base_dia/2*cos(THETA9) + y_0
plot(X_base9, Y_base9, color="cyan", linestyle="--")
# Pitch Circle
X_pitch9 = pitch_dia/2*sin(THETA9) + x_0
Y_pitch9 = pitch_dia/2*cos(THETA9) + y_0
plot(X_pitch9, Y_pitch9, color="magenta", linestyle="--")
# Offset Circle
X_offset9 = (offset_dia/2)*sin(THETA9) + x_0
Y_offset9 = (offset_dia/2)*cos(THETA9) + y_0
plot(X_offset9, Y_offset9, color="red", linestyle="--")
# Outter Circle
X_out9 = (outter_dia/2)*sin(THETA9) + x_0
Y_out9 = (outter_dia/2)*cos(THETA9) + y_0
plot(X_out9, Y_out9, color="brown", linestyle="--")
# Root Circle
X_root9 = (root_dia/2)*sin(THETA9) + x_0
Y_root9 = (root_dia/2)*cos(THETA9) + y_0
plot(X_root9, Y_root9, color="grey", linestyle="--")


############################
# Plot only One Tooth
plot(X_align,Y_align,color="black",linestyle="-",linewidth=1)
plot(X_align,Y_align,"b.")
base_radius = base_dia/2
pitch_radius = pitch_dia/2
offset_radius = offset_dia/2
outter_radius = outter_dia/2
root_radius = root_dia/2
xmm = x*m
emm = e*m
bmm = b*m
cmm = c*m
annotate(["Base Radius = $base_radius mm"],xy=(0,base_dia/2),ha="center",va="bottom",color="cyan",fontsize="x-small")
annotate(["Pitch Radius = $pitch_radius mm"],xy=(0,pitch_dia/2),ha="left",va="bottom",color="magenta",fontsize="x-small")
annotate(["Offset Radius = $offset_radius mm"],xy=(0,offset_dia/2),ha="right",va="bottom",color="red",fontsize="x-small")
annotate(["Outter Radius = $outter_radius mm"],xy=(0,outter_dia/2),ha="right",va="bottom",color="brown",fontsize="x-small")
annotate(["Root Radius = $root_radius mm"],xy=(0,root_dia/2),ha="center",va="bottom",color="grey",fontsize="x-small")
annotate(["Offset = $xmm mm"],xy=(0,offset_dia/2),ha="right",va="top",color="red",fontsize="x-small")
annotate(["Radius of Hob Edge = $cmm mm"],
        xy=(X_align[seg_root+int(seg_root_r/2)],Y_align[seg_root+int(seg_root_r/2)]),
        ha="right",va="top",color="blue",fontsize="x-small")
annotate(["Radius of Edge = $emm mm"],
        xy=(X_align[seg_root+seg_root_r+seg_involute],Y_align[seg_root+seg_root_r+seg_involute]),
        ha="left",va="bottom",color="blue",fontsize="x-small")
annotate(["Backlash = $bmm mm"],
        xy=(X_align[seg_root+seg_root_r+int(seg_involute/2)],Y_align[seg_root+seg_root_r+int(seg_involute/2)]),
        ha="left",va="bottom",color="blue",fontsize="x-small")
amm = a*m
dmm = d*m
hmm = amm+dmm
tmm = 2*X_align[seg_root+seg_root_r]
#temm = 2*X_align[seg_root+seg_root_r+seg_involute+int(seg_edge_r/2)]
temm = 2*X_align[seg_root+seg_root_r+seg_involute]
annotate(["Addendum = $amm mm"],xy=(X_out9[1],Y_out9[1]),ha="left",va="top",color="black",fontsize="x-small")
annotate(["Dedendum = $dmm mm"],xy=(X_root9[1],Y_root9[1]),ha="left",va="top",color="black",fontsize="x-small")
annotate(["Tooth Height = $hmm mm"],xy=(X_offset9[1],Y_offset9[1]),ha="left",va="top",color="black",fontsize="x-small")
annotate(["Approx Tooth Thickness = $tmm mm"],xy=(0,base_dia/2),ha="center",va="top",color="black",fontsize="x-small")
annotate(["Approx Tooth End Thickness = $temm mm"],xy=(0,outter_dia/2),ha="center",va="top",color="black",fontsize="x-small")

savefig("case02.svg")





############################
# Save case.csv
writecsv("case.csv",[Xccc Yccc])


############################
# Save case.geo
fileout = open("case.geo", "w")
println(fileout, "cl__1 = 1;")
num_Point_gear = length(Xccc)
num_Line_gear = num_Point_gear-1
num_Point_shaft = num_Point_gear + length(X_shaft) -1
num_Line_shaft = num_Point_shaft-1
for i=1:num_Point_gear
    println(fileout, "Point(",i,") = {",Xccc[i],", ",Yccc[i],", ",0,", ","cl__1};")
end
for i=num_Point_gear+1:num_Point_shaft
    println(fileout, "Point(",i,") = {",X_shaft[i-num_Point_gear],", ",Y_shaft[i-num_Point_gear],", ",0,", ","cl__1};")
end
for j=1:num_Line_gear
    println(fileout, "Line(",j,") = {",j,", ",j+1,"};")
end
println(fileout, "Line(",num_Line_gear+1,") = {",num_Line_gear+1,", ",1,"};")
for j=num_Line_gear+2:num_Line_shaft
    println(fileout, "Line(",j,") = {",j,", ",j+1,"};")
end
println(fileout, "Line(",num_Line_shaft+1,") = {",num_Line_shaft+1,", ",num_Line_gear+2,"};")
# Combine
array_LineLoop = [1:num_Line_shaft+1]
array_LineLoop[num_Line_gear+2:num_Line_shaft+1] = -array_LineLoop[num_Line_gear+2:num_Line_shaft+1]
array_LineLoop = string(array_LineLoop)
array_LineLoop = array_LineLoop[2:length(array_LineLoop)-1]
println(fileout, "Line Loop(",num_Line_shaft+3,") = {",array_LineLoop,"};")
println(fileout, "Plane Surface(",num_Line_shaft+3,") = {",num_Line_shaft+3,"};")
#Physical Line(13) = {8, 9};
#Physical Surface(14) = {11};
array_PhysicalLine = [num_Line_gear+2:num_Line_shaft+1]
array_PhysicalLine = string(array_PhysicalLine)
array_PhysicalLine = array_PhysicalLine[2:length(array_PhysicalLine)-1]
#println(fileout, "Physical Line(",num_Line_shaft+4,") = {",array_PhysicalLine,"};")
# Position for Normal Force
#position_NormalForce = (length(Xcc)-1) - (seg_root+seg_root_r-2) - (seg_involute-3)
#println(fileout, "Physical Line(",num_Line_shaft+5,") = {",position_NormalForce,"};")
# Fixed
#println(fileout, "Physical Surface(",num_Line_shaft+6,") = {",num_Line_shaft+3,"};")
println(fileout, "Mesh.ElementOrder = 2;")
println(fileout, "Mesh.RecombineAll = 1;")
close(fileout)


############################
# Save case.scr
fileout = open("case.scr", "w")
println(fileout, "spline")
for i=1:num_Point_gear
    println(fileout, Xccc[i],",",Yccc[i],",0")
end
println(fileout, Xccc[1],",",Yccc[1],",0")
println(fileout, "  ")
println(fileout, "circle")
println(fileout, x_0,",",y_0)
println(fileout, d_shaft/2)
close(fileout)


############################
# Save case.sif

fileout = open("case.sif", "w")
println(fileout, "Header")
println(fileout, "  CHECK KEYWORDS Warn")
println(fileout, "  Mesh DB \".\" \".\"")
println(fileout, "  Include Path \"\"")
println(fileout, "  Results Directory \"\"")
println(fileout, "End")
println(fileout, "")
println(fileout, "Simulation")
println(fileout, "  Max Output Level = 5")
println(fileout, "  Coordinate System = Cartesian")
println(fileout, "  Coordinate Mapping(3) = 1 2 3")
println(fileout, "  Simulation Type = Steady state")
println(fileout, "  Steady State Max Iterations = 1")
println(fileout, "  Output Intervals = 1")
println(fileout, "  Timestepping Method = BDF")
println(fileout, "  BDF Order = 1")
println(fileout, "  Solver Input File = case.sif")
println(fileout, "  Post File = case.ep")
println(fileout, "  Coordinate Scaling = Real 0.001")
println(fileout, "End")
println(fileout, "")
println(fileout, "Constants")
println(fileout, "  Gravity(4) = 0 -1 0 9.82")
println(fileout, "  Stefan Boltzmann = 5.67e-08")
println(fileout, "  Permittivity of Vacuum = 8.8542e-12")
println(fileout, "  Boltzmann Constant = 1.3807e-23")
println(fileout, "  Unit Charge = 1.602e-19")
println(fileout, "End")
println(fileout, "")
println(fileout, "Body 1")
println(fileout, "  Target Bodies(1) = 1")
println(fileout, "  Name = \"Body 1\"")
println(fileout, "  Equation = 1")
println(fileout, "  Material = 1")
println(fileout, "End")
println(fileout, "")
println(fileout, "Solver 1")
println(fileout, "  Equation = Result Output")
println(fileout, "  Single Precision = True")
println(fileout, "  Procedure = \"ResultOutputSolve\" \"ResultOutputSolver\"")
println(fileout, "  Save Geometry Ids = True")
println(fileout, "  Output File Name = case_output")
println(fileout, "  Output Format = Gmsh")
println(fileout, "  Exec Solver = After Timestep")
println(fileout, "End")
println(fileout, "")
println(fileout, "Solver 2")
println(fileout, "  Equation = Linear elasticity")
println(fileout, "  Procedure = \"StressSolve\" \"StressSolver\"")
println(fileout, "  Variable = -dofs 2 Displacement")
println(fileout, "  Exec Solver = Always")
println(fileout, "  Stabilize = True")
println(fileout, "  Bubbles = False")
println(fileout, "  Lumped Mass Matrix = False")
println(fileout, "  Optimize Bandwidth = True")
println(fileout, "  Steady State Convergence Tolerance = 1.0e-5")
println(fileout, "  Nonlinear System Convergence Tolerance = 1.0e-7")
println(fileout, "  Nonlinear System Max Iterations = 20")
println(fileout, "  Nonlinear System Newton After Iterations = 3")
println(fileout, "  Nonlinear System Newton After Tolerance = 1.0e-3")
println(fileout, "  Nonlinear System Relaxation Factor = 1")
println(fileout, "  Linear System Solver = Iterative")
println(fileout, "  Linear System Iterative Method = BiCGStab")
println(fileout, "  Linear System Max Iterations = 500")
println(fileout, "  Linear System Convergence Tolerance = 1.0e-10")
println(fileout, "  BiCGstabl polynomial degree = 2")
println(fileout, "  Linear System Preconditioning = Diagonal")
println(fileout, "  Linear System ILUT Tolerance = 1.0e-3")
println(fileout, "  Linear System Abort Not Converged = False")
println(fileout, "  Linear System Residual Output = 1")
println(fileout, "  Linear System Precondition Recompute = 1")
println(fileout, "End")
println(fileout, "")
println(fileout, "Equation 1")
println(fileout, "  Name = \"Equation 1\"")
println(fileout, "  Calculate Stresses = True")
println(fileout, "  Active Solvers(2) = 1 2")
println(fileout, "End")
println(fileout, "")
println(fileout, "Material 1")
println(fileout, "  Name = \"Steel (carbon - generic)\"")
println(fileout, "  Heat Conductivity = 44.8")
println(fileout, "  Youngs modulus = 200.0e9")
println(fileout, "  Mesh Poisson ratio = 0.285")
println(fileout, "  Heat Capacity = 1265.0")
println(fileout, "  Density = 7850.0")
println(fileout, "  Poisson ratio = 0.285")
println(fileout, "  Sound speed = 5100.0")
println(fileout, "  Heat expansion Coefficient = 13.8e-6")
println(fileout, "End")
println(fileout, "")
println(fileout, "Boundary Condition 1")
println(fileout, "  Target Boundaries(47) = 1922 1923 1924 1925 1926 1927 1928 1929 1930 1931 1932 1933 1934 1935 1936 1937 1938 1939 1940 1941 1942 1943 1944 1945 1946 1947 1948 1949 1950 1951 1952 1953 1954 1955 1956 1957 1958 1959 1960 1961 1962 1963 1964 1965 1966 1967 1968 ")
println(fileout, "  Name = \"Fixed\"")
println(fileout, "  Displacement 2 = 0")
println(fileout, "  Displacement 1 = 0")
println(fileout, "End")
println(fileout, "")
println(fileout, "Boundary Condition 2")
println(fileout, "  Target Boundaries(1) = 52 ")
println(fileout, "  Name = \"NormalForce\"")
println(fileout, "  Normal Force = -100")
println(fileout, "End")
println(fileout, "")
close(fileout)


############################
# Save case.bat
fileout = open("case.bat", "w")
println(fileout, "gmsh -2 case.geo")
println(fileout, "ElmerGrid 14 2 case.msh")
println(fileout, "move .\\case\\*.* .\\")
println(fileout, "rmdir case")
println(fileout, "ElmerSolver case.sif")
println(fileout, "gmsh case_output.msh")
close(fileout)


############################
# Save case.sh
fileout = open("case.sh", "w")
println(fileout, "#!/bin/sh")
println(fileout, "gmsh -2 ./case.geo")
println(fileout, "ElmerGrid 14 2 case.msh")
println(fileout, "mv ./case/* ./")
println(fileout, "rmdir case")
println(fileout, "ElmerSolver ./case.sif")
println(fileout, "gmsh ./case_output.msh")
close(fileout)
if OS_NAME == "Linux"
    run(`chmod 751 case.sh`)
end

이것을 실행하려면 우선 Julia가 실행 가능한 환경이어야 한다.
자신의 PC에 Julia가 깔려있지 않고, 또 깔고 싶지 않다면 http://juliabox.org 사이트에서 클라우드 서비스를 제공해 주므로 거기서 소스코드를 복사해서 붙여넣은 후 실행해도 된다.
실행 후 생성된 파일들을 다운로드 받을 수 있기 때문에 전혀 문제가 안된다.

(사실 이것도 불편하기 때문에, 추후에 윈도우용 및 리눅스용 독립 실행파일을 생성해서 제공해 주려고 한다.  아직 Julia v0.3.9 에서는 독립 실행 파일 생성이 불가능하고, 조만간 나올 Julia v0.4 부터는 독립실행파일을 생성해 줄 수 있는 패키지가 안정화되어 나올 예정이기 때문에 시간이 좀 필요하다.)

그리고, 이 프로그램은 기어 파라미터를 설정해 준 별도의 파일을 통해 입력받도록 해 놨기 때문에 입력용 파일을 따로 작성해 둬야 한다.  파일 이름은 input4fgpg.csv 로 해야 한다.
이 파일은 콤마로 필드가 구분된 텍스트파일이다.  엑셀 같은 스프레드 시트에서 읽으면 다루기 쉽다.

내용은 아래와 같다.

Module,m,1
Teeth Number,z,24
Pressure Angle [Deg],alpha_0_deg,20
Offset Factor,x,0.2
Backlash Factor,b,0.025
Addendum Factor,a,1.2
Dedendum Factor,d,1.0
Radius Factor of Edge Round of Hob,c,0.3
Radius Factor of Edge Round of Tooth,e,0.2
Center of Gear,x_0,0
Center of Gear,y_0,0
Segmentation Numbers,seg_circle,360
Segmentation Numbers,seg_involute,20
Segmentation Numbers,seg_edge_r,4
Segmentation Numbers,seg_root_r,10
Segmentation Numbers,seg_center,2
Segmentation Numbers,seg_outter,6
Segmentation Numbers,seg_root,6

기어 사양을 변경하고 싶다면, 위 내용 중에 다른 것은 건드리지 말고 매 행의 맨 마지막 필드의 숫자만 적절히 바꿔주면 된다.
 기어 사양에 관한 자세한 설명은 나중에 버전업 하고 나서 쓸만하게 되면 자세한 도큐먼트를 작성하려고 하므로 일단 생략한다.


아무튼 PC에 Julia 환경이 구축된 경우의 작업 예시를 순서대로 설명해 보겠다.
본 예제에서는 윈도우7 환경이다.

일단 작업하려고 하는 디렉토리를 만들고, 소스코드 및 입력파일을 위와 같이 넣는다.
Juliabox.org 에서 작업할 경우도 역시 마찬가지다.  작업하려는 디렉토리를 만들고 똑같이 넣어두면 된다.

 FGPG_V11.jl 소스코드를 실행하기 위해서는 크게 3가지 방법이 있다.

그냥 cmd 즉 소위 말하는 도스창에서

julia FGPG_V11.jl 

이렇게 커맨드를 직접 때려주는 방법이 있고, 또 다른 방법은 위 그림과 같이 JUNO IDE를 실행시켜서 여기서 실행해 주는 방법이 있다.
마지막으로는 Juliabox.org의 웹브라우저 상에서의 개발환경에서 실행해 줄 수 있다.  Juliabox.org는 기본적으로 IPython Notebook 기반이기 때문에, IJulia Notebook을 자신의 PC에 설치해서 똑같이 해도 된다.

본 예제에서는 위 그림 처럼 JUNO 상에서 실행한다.

 실행 전에 input4fgpg.csv 파일도 열어서 내용 확인하고 적절히 수정해 준 후에, 다시 재 실행하면서 기어 형상 변화를 확인해 갈 수 있을 것이다.

JUNO 에서 소스를 실행하려면, 아래 그림처럼 메뉴 중에서 "Eval All"을 선택해 주면 된다.

참고로, 본 소스코드는 외부 패키지는  PyPlot만 사용한다.
가급적 외부 의존성을 줄여서 간단하게 만들기 위해서이다.
Julia는 JIT 방식이기 때문에, 최초 실행시에는 로딩 및 내부 컴파일 시간이 좀 걸린다.
하지만 두 번째 실행 부터는 상당히 빠르게 수행된다.

실행이 완료되면, 우선 아래와 같은 2개의 그림이 별도 창으로 뜬다.

이것을 보면서 기어의 형상 및 결과로 나오는 주요 치수들을 확인할 수 있다.
기어 사양이 적절치 않다면, input4fgpg.csv 안의 파라미터 값들을 조절해 주면서 반복 실행해 주면 된다.

아울러, 실행이 완료되면 작업 디렉토리 안에 여러 가지 파일들이 만들어져 있음을 볼 수 있다.

우선 case.csv 파일은, 생성된 기어 형상의 원시(Raw) 데이타이다.
이것을 엑셀에서 읽어서 그래프로 그대로 그려줄 수 있다.
또는 Matlab이나 Octave, Scilab 등에서 읽어서 이용할 수도 있다.


case01.svg 및 case02.svg 파일들은 앞서 나타났던 기어 형상 그림들이 파일로 저장된 것이다.
벡터 그림 파일이기 때문에, Inkscape 또는 Adobe Illustrator 같은 것으로 읽어서 활용할 수 있다.
아래의 그림은 Inkscape에서 읽어본 것이다.

case.scr 파일은 AutoCAD Script 파일이다.
AutoCAD 대신 DraftSight에서도 실행 가능하기 때문에, DraftSight에서 실행한 예를 아래 그림과 같이 볼 수 있다.

DraftSight를 실행해서 빈 화면이 나오며느, 화면 아래쪽에 "도면요소 스냅"이 보통 기본적으로 눌러진 상태로 되어 있다.  이것을 마우스로 한 번 눌러줘서 해제 상태로 바꾼다.

"도면요소 스냅"을 해제하는 이유는, 눌러진 상태에서 case.scr 파일을 실행시키면 제대로 그려지지 않고 이상하게 나오기 때문이다.  AutoCAD 에서도 똑같다.

아무튼, 이렇게 실행하면 약간의 수초~십수초 정도의 시간을 소요한 후에 기어 형상이 그려져 있음을 확인할 수 있다.
기어의 치형은 하나의 Spline으로 그려지도록 해 놨다.
또 Addendum 및 Dedendum 값을 비교해서 Internal & External Gear를 자동 판별해서 그려지도록 해 놨다.

이것을 case.dxf로 저장한다.  버전은 R2000~2002 버전으로 적당히 옛날 버전으로 해 두는 습관이 좋겠다.  왜냐면 구버전일수록 여러 소프트웨어를 가리지 않고 잘 열리기 때문이다.

한가지 팁은, R12 버전의 *.dxf 파일로 저장시키면 Spline이 많은 짧은 직선으로 구성된 Polyline으로 바꿔져서 저장된다.  R12 버전에서는 Spline 사양이 없기 때문이다.
이런 특징을 이용해서, Spline 궤적을 생성할 수 없는 공작기계로 가공할 경우에 (구형 와이어 방전 가공기나 밀링머신들은 Spline 궤적 보간 기능이 없고, 오로지 원호 보간 기능만 지원된다.  최신 기계들은 대부분 Spline 보간 기능이 있다.) 대응 가능할 것이다.

한편, case.bat 및 case.sh 파일은 각각 윈도우 및 리눅스용 해석 스크립트이다.
이 파일들을 실행시키면 PC에 이미 Gmsh 및 Elmer가 설치되어 있고 또 path가 잡혀 있는 경우에 한해서 이상없이 FEM 해석이 실시될 것이다.

Batch 파일의 스크립트 내용은 별것 없다.  위 그림과 같다.
순서대로 명령을 실행하는 것이다.
리눅스용 case.sh 파일 역시 같은 동작을 하도록 리눅스용 Bash 쉘에 맞춰서 내용이 들어있다.

gmsh -2 case.geo 명령은, case.geo 형상 파일을 이용해서 gmsh 프로그램으로 2차원 매쉬를 생성하라는 것이다.
매쉬 생성에는 수십초 이상이 소요될 것이다.  커맨트창이 떠서 진행 상황이 확인될 것이다.

gmsh 작업이 끝나면, case.msh 파일이 생성된다.

참고로, case.msh로 만들어지는 매쉬의 특성은, Delaunay 삼각형 자동 생성 알고리즘으로 2nd Order 조건으로 먼저 생성된 후에, Recombine 옵션을 줘 놨기 때문에 삼각형들을 합쳐서 사각형으로 바꿔준다.  따라서 결과물로 나오는 매쉬는 2nd Ordered Quad Mesh가 된다.

매쉬의 조밀함은 input4fgpg.csv 내용 중에서 Segmentation 값들을 올려주면 더 조밀하게 할 수 있다.  대신 해석시간은 기하급수적으로 증가될 것이다.

아무튼 이어서  곧바로 ElmerGrid 명령이 실행된다.  ElmerGrid는 case.msh 파일을 Elmer 전용의 매쉬 포멧으로 바꿔준다.  Elmer용 매쉬 파일은 단일 파일이 아니고, 속성별로 4개의 파일로 저장된다.

따라서 매쉬 생성 결과는, 위 그림과 같이 나타난다.
이어서 곧바로 ElmerSolver가 작동해서 해석을 시작한다.
각종 해석 조건을 설정한 ElmerSolver용 입력파일인 case.sif 파일은 이미 FGPG_V11.jl이 만들어 뒀다.
재질(Carbon Steel), 경계조건의 위치, 힘의 강도(100N)은 상수로 고정시켜 놨기 때문에 소스를 보고 수정해야 바꿀 수 있다.  일단은 잘 되는지 보기 위해 편의상 상수로 뒀다.

아무튼 이에 따라 ElmerSolver가 수행된다.  커맨드창에 나타나는 숫자들은 수렴 오차이다.  수렴오차를 계속 줄여가다가 설정된 오차값 안에 들어오면 해석이 완료되도록 되어 있다.

한참 기다리면 최종적으로 완료된다.  본 예제에서는 해석시간이 406.28초 걸렸다.

이제 어떤 파일이 생겼나 확인해 보면, case_output.msh 파일이 생겨나 있는 것을 확인할 수 있다.  물론 ElmerSolver 결과 파일인 case.ep 역시 함께 생성되어 있다.
이들 파일이 생성된 후에 자동으로 Post Processor로서 gmsh가 실행되어 결과 파일을 띄워준다.

gmsh는 처음 다뤄보는 사람이면 GUI 체계가 좀 생소해서 막막할 수 있는데, 기능들 자체는 매우 단순하다.  이리저리 보면서 건드려 보면 쉽게 습득 가능할 것이다.
일단 본 화면에서는 기본적으로 결과데이타 8가지가 한꺼번에 겹쳐서 보이므로 번잡스럽다.
6번이 VonMises 응력 데이타이므로, 이것을 제외하고 전부 Uncheck 한다.

일단 메뉴에서 Option을 선택해서 화면 상태를 조정해 줄 수 있도록 옵션을 건드려 본다.

Mesh의 Visability 중에서 Surface Edges, Volume Edges를 Uncheck 하면 좀 깔끔해진다.

그리고, Filled iso-values를 선택해서 등고선식으로 보이도록 해 본다.

등고선 단계는 50단계 정도로 해 본다.

그러면 좀 볼 만 하다.  다만 몇 개의 거슬리는 포인트들이 보인다.

거슬리는 포인트는 위 그림처럼 Points를 Uncheck 하면 안보이게 된다.

결과 확인.  그림은 캡춰하거나 다른이름으로 저장할 때 그림포멧으로 지정하면 된다.

이제 7번 Displacement를 본다.

상태 확인하고...

등고선식으로 보이도록 해 보고...

Displacement Factor를 크게 높여서 변형되는 형상을 과장해서 확인해 본다.
이런 식으로 결과 확인이 가능하다.


한 편, 앞서 만들어둔 기어 형상 2D CAD 파일인  case.dxf 파일을 3D CAD에서 읽어들여 모델링해 본다....

일단 FreeCAD에서 잘 되는지 확인.

FreeCAD를 최초로 띄우면 위 그림처럼 나온다.

빈 파일에 case.dxf를 Import 시킨다.

최초에 FreeCAD에는 dxf Reader가 내장되어 있지는 않으나, 우와 같은 창이 뜨면 Yes 해 주면 자동으로 스크립트를 설치해서 읽어들인다.  dxf Reader 부분의 저작권 조건이 달라서 FreeCAD 설치본에 포함하지 못했다고 한다.  (FreeCAD v0.15 기준)

읽어들이는데 성공하면 위 그림과 같다.
간단하게는 이것을 Extrude 시켜서 모델링 하면 된다.

일단은 외곽 라인 부분을 Extrude 해 주는데 일단 Solid가 생성되도록 체크해 준다.

또 기어 치형 부분도 Solid로 Extrude 해 준다.

마지막으로, 생성된 2개의 Extrude Feature를 Boolean Cut 해 준다.
그러면 이상없이 만들어진 모습을 볼 수 있다.
이 파일을 *.step 포멧으로 저장 가능하다.

한편, 프로페셔널하게 쓰는 상용 3D CAD에서도 잘 되는지 확인해 본다...

잘 된다.  모델링 과정 설명은 생략한다.  대충 해 본 거라서..
재밌게 응용해 보기 위해 각도를 살짝 틀어서 Scew를 줘 봤다.  헬리컬기어 느낌이 난다.

이런식으로 응용하면, 상당히 엄밀한 Bevel Gear 및 Bevel Spiral Gear 모델링도 가능할 것이다.




추후에 본 프로그램의 버전업 계획은 아래와 같다.


1. 언더컷 자동 감지 및 겹치는 라인 자동 제거
2. 이끝단 살 소멸 자동 감지 및 겹치는 라인 자동 제거
3. 치형 수정(End Tip Relief) 기능 추가 (Linear Type)
4. 해석 속도 증가를 위해 기어 전체를 해석하지 않고, 3개 정도의 Teeth만 살려서 해석할 수 있도록 수정
5. Theory 및 User Manual 작성


요즘 시간이 점점 쪼들리고 있는데, 신속하게 가능할지는 잘 모르겠다.

이상 끝!

Julia Code : COMBI REDUCER with Planetary & Harmonic Drive v05

MIT에서 만든 Julia Language에 대해서 관심을 계속 가지고 있었는데, 이제 v0.3 이상으로 올라가면서 꽤 쓸만해 진 모양이다.  특히 아마존 웹 서비스에다가 이걸 올려서 클라우드 서비스까지 해 주기 때문에 어렵게 PC에 개발환경 구축하느라 진땀 빼지 않아도 되어서 더 좋다.

나에게 필요한 것은, 컴퓨터공학적인 최적화라던가 대규모 데이타를 효과적으로 다루는 것이 아니고, 높은 수준의 공학용 계산기 내지는 간단한 스크립트 엔진이다.

여러가지 많이 검토해 봤는데 Julia가 역시 최신기술의 집약체이다 보니 모든 조건에 다 부합하는 것 같다.

아래는 어떤 문제를 풀기 위해 Julia를 실무적으로 사용해 본 예이다.
원래는 Matlab m-code로 작성해서 FreeMat, Octave 따위로 돌려보던 것인데, Julia로 테스트해 보기 위해 컨버젼도 해 보고 완전히 새로 작성도 해 보았다.

m-code를 Julia code로 자동 컨버젼해주는 패키지도 Julia 패키지 라이브러리에 있으나, 일단 직접 수동으로 컨버젼 해 보았더니 함수 이름 바꿔준다던가 정도의 문법 수정으로 20분만에 끝났다.

그리고 Julia에 좀 더 부합하도록 몇가지 기능을 추가하고 수식도 새로 적용하고 해서 새로 만들었다.


---------------------------------------------------------------------------------------------------------------


본 Julia Code는 기존의 버전을 Julia 문법에 좀 더 부합하게 고치고, Interactive 기능을 추가한 버전이다.

최외곽 서큘러 스플라인 잇수(Number of Teeth)와 그 안쪽으로 들어가면서 조립되는 각 기어의 잇수의 차이를 기준으로 사양이 자동으로 잡히도록 하였고, Sun Gear의 반지름이 음수로 나오면 색상이 붉은 색으로 변하면서 경고 문구가 표시되도록 하였다.

기타 파라미터들은 관계를 조정하여 기하학적인 모순 상황이 발생하지 않도록 하였다.

타원 변형된 Flex Spline의 원주길이는 다음 공식을 기본으로 하였다.


* Annotations

$a$ : 긴 축의 지름 (Known)

$x$ : 짧은 축의 지름 (Unknown)

$\varepsilon$ : 이심률 (Known)

$L$ : 타원의 원주길이 (Known)


* Notations

$\varepsilon = \sqrt{ 1 - \frac{x^2}{a^2} }$

$L = 2aE(\varepsilon) = \pi a \left[  1 - \left( \frac{1}{2} \right)^2 \varepsilon^2 -  \left( \frac{1 \cdot 3}{2 \cdot 4} \right)^2 \frac{\varepsilon^2}{3} - \left( \frac{1 \cdot 3 \cdot 5}{2 \cdot 4 \cdot 6} \right)^2 \frac{\varepsilon^2}{5} - ... \right]$


* Remarks

$E(\varepsilon)$은 제2종 타원정적분 연산기호이다.  타원정적분은 엄밀하게 풀이하기가 어렵기 때문에, 이를 풀어쓴 우변의 무한급수 식을 이용하기로 한다.  급수의 항이 늘어날 수록 해(Solution)의 개수가 늘어나기 때문에 적합한 해를 판별하는 판단이 더 필요하게 되므로, 간소하게 하기 위해 다음과 같이 첫 번째 항만 살린 근사식을 만들었다.

$L = \pi a \left[  1 - \left( \frac{1}{2} \right)^2 \varepsilon^2 \right]$

해를 신속하게 구하기 위해 SageMath를 이용하였다.
(SageMath 대신 WolframAlpha를 이용해도 되나, 해가 복소수로 나오는 등 형태가 다를 수도 있다.)

a,L = var('a,L')
eq = L == pi*a* (  1-((1/2)^2)*(( sqrt(1-(x^2)/(a^2)) )^2)   )
eq
print solve(eq,x)

SageMath를 이용하여 얻어진 해는 다음과 같다.

L == 1/4*pi*a*(x^2/a^2 + 3)
[
x == -sqrt(-3*a^2 + 4*L*a/pi),
x == sqrt(-3*a^2 + 4*L*a/pi)
]

2개의 해가 얻어졌는데, 이 중에서 양수에 해당하는 두 번째 것을 취한다.  이 식을 Julia Code에 집어넣고 사용한다.


* Julia Coding

PyPlot 패키지는 Python 라이브러리 중에서, Matlab과 거의 비슷하게 Plot을 해 주는 MatPlotLib 기능을 Julia에서 사용할 수 있도록 해 주는 라이브러리이다.  따라서 이것은 내부적으로는 Python MatPlotLib가 그대로 동작하며, Julia Code에서 좀 더 사용하기 쉽도록 Matlab 문법과 거의 유사하게 간소화되어 있다.  기계공학도에게 Julia가 Python보다 접근하기 쉬운 점은, 굳이 객체지향 방법으로 코딩하지 않아도 된다는 점이다.  Matlab을 다루듯이 편하게 쓰면 된다.

Interact 패키지는 매우 간단하게 기본적인 GUI를 만들 수 있도록 해 준다.  여기 포함된 @manipulate 매크로를 이용하여, 자유롭게 조절하고 싶은 파라미터들을 설정해 주기만 하면 된다.  내부적으로는 무한루프를 돌리는 간단한 형태다.

#######################################################
# COMBI REDUCER with Planetary & Harmonic Drive
# V05
# 2015.05.18
# Dymaxion.kim
#######################################################
 
using PyPlot
using Interact
 
# Basic Parameter
alpha=[0:0.02:2*pi]
max_n_cs = 300
max_dn_fs = 20
max_dn_ig = 20
max_dn_sg = max_n_cs - max_dn_fs - max_dn_ig
# Making New Figure
f=figure(figsize=(8,8))
# Interact Macro
@manipulate for m = 0.1:0.1:2, # Module of The Gear
    n_cs = 50:1:max_n_cs, # Number of Teeth for Circular Spline
    dn_fs = 1:1:max_dn_fs, # Different Number of Teeth for Flex Spline
    dn_ig = 1:1:max_dn_ig, # Different Number of Teeth for Internal Gear
    dn_sg = 1:1:max_dn_sg # Different Number of Teeth for Sun Gear
    
    withfig(f) do
        grid("on")
        axis("auto")
        title("COMBI REDUCER v05")
        
        # Draw Circular Spline
        r_cs = m*n_cs/2
        plot( r_cs*sin(alpha), r_cs*cos(alpha), linewidth=2, color="Magenta" )
        annotate(["n_cs",int(n_cs)],xy=(0,r_cs),ha="center",va="bottom", color="magenta")
        
        # Draw Flex Spline in Idle Status
        r_fs = m*(n_cs-dn_fs)/2
        plot( r_fs*sin(alpha), r_fs*cos(alpha), linestyle="--", color="cyan" )

        # Draw Internal Gear in Idle Status
        r_ig = m*(n_cs-dn_fs-dn_ig)/2
        plot( r_ig*sin(alpha), r_ig*cos(alpha), linestyle="--", color="cyan" )
        
        # Draw Sun Gear
        r_sg = m*(n_cs-dn_fs-dn_ig-dn_sg)/2
        n_sg = 2*r_sg/m
        if r_sg > 0
            plot( r_sg*sin(alpha), r_sg*cos(alpha), linewidth=2, color="blue" )
            annotate(["n_sg",int(n_sg)],xy=(0,0),ha="center",va="center", color="blue")
        else
            plot( r_sg*sin(alpha), r_sg*cos(alpha), linewidth=4, linestyle="--", color="red" )
            annotate(["n_sg is in minus",int(n_sg)],xy=(0,0),ha="center",va="center", color="red")
        end
        
        # Draw Planet Gears
        t_fs = r_fs-r_ig # Thickness of Flex Spline's Plastic
        r_pg = (r_cs-t_fs-r_sg)/2
        plot( r_pg*sin(alpha), r_pg*cos(alpha)+(r_sg+r_pg), linewidth=2, color="magenta" )
        plot( r_pg*sin(alpha), r_pg*cos(alpha)-(r_sg+r_pg), linewidth=2, color="magenta" )
        n_pg = 2*r_pg/m
        annotate(["n_pg",int(n_pg)],xy=(0,r_sg+r_pg),ha="center",va="center", color="magenta")
        annotate(["Thickness",t_fs],xy=(0,r_sg+2*r_pg),ha="center",va="top", color="blue")
        
        # Draw Flex Spline in assembled
        L_fs = 2*pi*r_fs
        Dl_fs = 2*r_cs
        Ds_fs = sqrt(-3*Dl_fs^2 + 4*L_fs*Dl_fs/pi)
        plot( (Ds_fs/2)*sin(alpha), (Dl_fs/2)*cos(alpha), linewidth=2, color="blue" )
        n_fs = 2*r_fs/m
        annotate(["n_fs",int(n_fs)],xy=(-0.9*r_fs,0),ha="left",va="bottom", color="blue")

        # Draw Internal Gear in assembled
        L_ig = 2*pi*r_ig
        Dl_ig = Dl_fs - 2*t_fs
        #Ds_ig = sqrt(-3*Dl_ig^2 + 4*L_ig*Dl_ig/pi)
        Ds_ig = Ds_fs - 2*t_fs
        plot( (Ds_ig/2)*sin(alpha), (Dl_ig/2)*cos(alpha), linewidth=2, color="blue" )
        n_ig = 2*r_ig/m
        annotate(["n_ig",int(n_ig)],xy=(-0.9*r_ig,0),ha="left",va="top",color="blue")
        
        # Reduction Ratio
        RRp = (n_sg+n_ig)/n_sg
        RRh = n_fs/(n_fs-n_cs)
        RRt = RRp * RRh
        annotate(["Total Reduction Ratio",RRt],xy=(0,-r_cs),ha="center",va="top", color="black")

    end
end

아래 그림은 http://juliabox.org에서 구동한 결과이다.

검증결과 동일한 조건의 파라미터를 입력했을 때 기존 버전과 동일한 결과를 얻을 수 있었다.