개인 위키를 어디다 만들까?

개인 위키를 어디다 만들까?


블로그에다가 각종 테스트한 것들을 계속 올려보니 문서들의 조직화가 안되기 때문에 한계성을 느끼게 되었다.  또 html을 직접 건드리면서 편집하는것도 이제는 염증이 생긴다.
그래서 위키 쪽으로 생각을 해 보았는데, 집에다가 서버 돌릴 생각하니 너무 끔찍해서 그냥 공짜로 위키를 쓸 수 있는 방안을 모색해 보기로 했다.

* 첫번째 시도

일단 맨 처음에는 호스팅 업체 쪽에서 제공하는 무료 플랫폼을 알아봤는데, 고르고 골라 본 게 호스팅어(Hostinger) 였다.
무료계정 신청하고 들어가 보니, 온갖 패키지들을 원클릭으로 설치해서 사용할 수 있도록 엄청나게 잘 되어 있었다.  미디어위키 역시 그냥 누르고 설정만 몇 개 해 주니까 그냥 쨔잔 하고 되는게 아닌가!!!
다만 서버가 영국에 있다고 나오는지라 반응 속도는 조금 느린 편이지만 못 쓸 수준은 아니고 상당히 좋았다.  제공 용량도 충분하고...
하지만 문제점은, 처음 셋팅하고 이것저것 건드려보다 보니 갑자기 접속이 끊어져 버리는 현상이 나타났다.  CPU 점유율이 올라가면 서버쪽에서 트리거링 되어서 자동으로 접속 끊어버린다고 하던데 이것도 그런 경우인가 싶기도 한데, 이런 현상이 자주 나타난다면 신뢰성 측면에서 도저히 안심이 되질 않는다.  데이타가 통째로 날아가 버리는 사태가 발생할 확률이 아무래도 무시 못할 듯 하다.
또 호스팅 형식으로 한다면 패키지 버전 관리다 뭐다 해서 골치아픈 일이 많을 것 같다.
결국 보류해 두기로 했다.

* 두번째 시도

꼭 개인 계정으로 꽁꽁 숨겨둘 필요는 없으므로, 다른 관점에서 접근해 보기로 했다.
바로 위키독스(WikiDocs)로 온라인 책을 써 보는 거다.
위키독스 만든 분도 꽤 괜챦은 개념인 같아 보여서 믿음도 가고, 또 내용이 풍부해지면 공개해서 기여도 할 수 있을 것 같아서 좋아 보인다.
게다가 마크다운 태그 기반으로 먹일 수 있어서 더 좋다.  세심하게 개발되어 있어서, 그림을 업로드해서 붙여주는게 거의 원클릭으로 된다.
단점은 위키독스 같은 서비스는 거의 전부 소프트웨어 개발자들 위주의 것들이다보니 기계공학 쪽 주제로 하면 왠지 혼자만 튈 것 같아 좀 쪽팔린다.
또 너무 퍼블릭한 서비스라서 완성도가 낮은 결과물일 경우에도 쪽팔림이 있을까 두렵네...  책을 쓴다는 느낌은 아무래도 부담스럽다.
또하나 치명적인 단점은, 데이타를 내 쪽으로 통째로 백업받는 기능이 안 보인다.  그냥 전부 온라인에서 작업하는 것만 가능한 것 같다.  이렇게 하는게 간단하긴 한데 그래도 뭔가 안되는게 있으니 찜찜한 느낌이다.

* 세번째 시도

깃허브(GitHub)에서도 위키가 있으니 그걸 써 볼까 싶은 생각이 들었다.  실제로 그런식으로 쓰는 사람도 있다고 하니...
이걸 쓰면 부수적인 효과로 Git 공부도 되고 좋겠다 싶다.
그림 올리는 건 Issues 기능에다가 그림 올린 다음에 그걸 퍼다 쓰면 되는 꼼수가 있는 것도 알았다.  또 자체 홈페이지 운영도 가능하므로 웹을 잘 다룰 수 있게 된다면 그것도 좋을 듯 하다.  일단 최고로 간단하게 설정해서 만들어 봤는데, 나름대로 봐줄만 하긴 하다.  문제는, 나는 웹개발에 별 관심이 없다는 거다.  앵귤라, 부트스트랩, 지킬 등 온갖 자바스크립트 프레임웍들을 이용해서 뭔가를 만들거나 다뤄야 하는데 영 자신도 없고 의지도 없다.  이거 공부해서 팔 거면 진작에 웹 개발자 했겠지... ㅋㅋ
그냥 있는거 갖다 써서 운영하는 것도 가능하므로 욕심만 안 내면 되긴 되겠다.
확실한 건 깃허브가 요즘 대세이고, 그게 아주 오래 갈 것 같다.
대세에 편승한다는 점에서는 장점이 있을지도 모른다.


이상 세가지 경우에 대해서 시험삼아 기본적인 셋팅은 전부 다 해 봤다.
전부 다 운영하는건 미친 짓이므로 하나를 골라야 될 텐데, 어느쪽이 좋을지 아직 확실히 잘 모르겠지만 현재까지는 깃허브 쪽으로 마음이 기우는 듯 하다.

마크업 태그는 사용해 본지 얼마되진 않았는데, 사내 자료 축적용으로 운영하는 미디어위키에다가 자료들을 작성해 보니 꽤 좋다는 느낌이 들었다.
위키독스 및 깃허브 쪽은 미디어위키의 것보다 좀 더 발전된 마크다운을 사용하도록 되어 있는데, 써보니까 좀 더 간결해서 더 좋다는 느낌이 들었다.

아무튼 점수를 내 보니까, 첫번째 호스팅 서비스 쓰는것은, 확장성이나 자유도 면에서는 우월하지만, 유지관리에 들어가는 리소스가 너무 많을 것이므로 아무래도 피하는게 맞는 것 같다.
위키독스는 내용작성 그 자체에만 집중하는데는 좋겠지만 부수적 효과는 기대하기 어렵다.
깃허브는 이왕 하는 김에 Git도 좀 다뤄보고 하면 일석이조가 될 것 같다.  신뢰성도 가장 월등하고 용량제한도 없고...  또 clone 해서 오프라인으로 작업해도 무방하고...
어차피 독자를 염두에 두고 기록을 남기는 것도 아니므로, 깃허브로 가는게 맞지 않나 하는 생각이 든다.

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 기본 환경 구축 방법 간단 요약

* 윈도우에서 Julia 기본 환경 구축 방법

1. Python 환경 구축을 먼저 한다.

  1.1. 정신건강을 위해 그냥 일반적인 Python 설치본을 사용하지 말고, Anaconda 설치본을 사용하도록 한다.  Anaconda 설치본에는 Julia에서 주로 갖다쓰게 될 과학계산용 패키지들이 전부 통합되어 있기 때문에 한 방에 끝난다.  Python 깔고 Matplotlib 깔고 뭐 깔고 하는 과정 생략 가능.
  1.2. Python 버전은 반드시 2.7로 한다.  3.X 버전은 사용하지 말 것.
  1.3. Anaconda 다운로드 받는 곳 ::: http://continuum.io/downloads
  1.4. Anaconda를 깔 때, 반드시 기존에 이미 깔려있는 Python이 있다면 가급적 말끔히 지워주고 까는 것이 좋겠다.
  1.5. 설치 과정 중에 몇가지 옵션 체크가 기본으로 되어 있는데 건드리지 말고 그냥 디폴트로 깐다.
  1.6. 설치 장소는 c:\Anaconda 로 단순하게 고정해 두는 것이 좋겠다.
  1.7. 설치후 시스템 설정에서 Path를 잡아주도록 한다.  순서는 다음과 같다.
  윈도우 제어판 → 시스템 및 보안 → 시스템 → 설정변경 → 고급(탭) → 환경 변수 까지 찾아들어간 다음, 사용자 변수에 다음 내용을 추가해 준다.

PYTHONHOME = c:\Anaconda
PYTHONPATH = c:\Anaconda\Lib

  1.8. 또 Anaconda Launcher를 실행해서 패키지를 업데이트 해 주는 것이 좋겠다.

2. Julia를 깐다.

  1.1. Julialang.org 가서 받아다 깔면 된다.
  1.2. 설치 위치는 c:\julia 으로 고정한다.
  1.3. 설치 직후 터미널에서 julia 쳐 보고 작동되는지 확인.
  1.4. Julia 프롬프트상에서 Pkg.add("PyPlot") 해서 이상없이 패키지 추가되는지 확인.  PyPlot은 Python 및 Matplotlib가 작동되어야 이상없이 되기 때문에 이상여부 확인하기에 좋음.
  1.5. 설치후 시스템 설정에서 Path를 잡아주도록 한다.  순서는 다음과 같다.
  윈도우 제어판 → 시스템 및 보안 → 시스템 → 설정변경 → 고급(탭) → 환경 변수 까지 찾아들어간 다음, 시스템 변수 중에서 Path 항목에 다음 내용을 추가해 준다.

;c:\julia;c:\julia\bin

  새롭게 바뀐 Path를 적용하기 위해서는 윈도우 시스템을 재부팅해야 하는 것 같다.
  이후에

echo %path%




를 쳐서 잘 등록되어 먹어들어갔는지 확인해 보고, 아무 디렉토리에서 julia를 쳐서 실행 되는지 확인해 본다.

3. Juno를 깐다.  (Julia 전용 에디터)

  3.1. Julialang.org 가서 받아다 깔면 된다.
  3.2. 압축을 그냥 풀면 된다.  위치는 c:\Juno 로 한다.
  3.3. 이제 실행해서 쓰면 된다.

* 우분투/데비안 리눅스에서 Julia 기본 환경 구축 방법

1. Python 환경 구축을 한다.

  1.1. 웬만하면 python은 뭐 이미 깔려 있을테니 안해도 될 거고...
  1.2. MatPlotLib 랑 IPython을 깔아준다.

sudo apt-get install ipython-notebook python-matplotlib

2. Julia를 깐다.

  2.1. Julia의 레포지토리를 등록한다.

sudo add-apt-repository ppa:staticfloat/julianightlies  (최신버전)
sudo add-apt-repository ppa:staticfloat/juliareleases   (안정화버전)
sudo add-apt-repository ppa:staticfloat/julia-deps

  2.2. Julia를 깔고 확실히 조져준다.

sudo apt-get update
sudo apt-get install julia
sudo apt-get update
sudo apt-get upgrade

  2.3.  깔고 나서 julia 쳐서 작동 잘 되는지 확인해 본다.

3. Juno를 깐다.  (Julia 전용 에디터)

  3.1. Julialang.org 가서 받아다 깔면 된다.
  3.2. 압축을 그냥 풀면 된다.  위치는 내 경우에는 보통 /home/아이디/APP/juno 이런식으로 한다.
  3.3. 아무데서나 실행 가능하도록 실행파일의 심볼릭 링크를 걸어준다.

sudo ln -s /home/아이디/APP/juno/Juno /usr/bin/Juno

  3.4. 여기서 실행명령 이름이 Juno인데, 첫 글자가 대문자이므로 유의한다.
  3.5. 가능하면 /usr/share/applications 디렉토리 안에 Juno.desktop 파일도 하나 만들어 주면 메뉴에 등록되므로 편리하겠다.


끝.


Juno를 전용 에디터로 추천하는데, 다른 일반적인 에디터에 비해서 나은 점은 딱 하나다.
별다른 설정 없이도 그냥 실행해서 쓰면 된다는 거다.
Julia 코드 실행해 주는 건 다 기본 설정 되어 있고, 또 주요 예약어들의 자동완성 기능 같은것도 된다.
나머지는 허접하다.

Juno가 싫으면 IPython의 notebook Jupiter 위에 올라타서 구현된 Julia notebook도 나쁘지 않겠지만, 사실 이건 Juliabox.org 가서 써 보면 체험할 수 있으니 굳이....

윈도우에서 Notepad++ 에디터 사용자를 위한 설정 파일도 함께 Julia에 동봉되어 있는 것 같은데 확실히 확인은 안 해 봤다. (개발중인 Julia 0.4에는 있는데, 현재 안정화 배포버전인 Julia 0.3.9에 포함되어 있지는 않은 듯)


이상 환경이 구성되면, 이제 자유롭게 막 쓰면 된다.
Julia는 헤프고 쉬운 여자다.  딱 내 타입이다.