2016년 6월 29일 수요일
2016년 4월 9일 토요일
STEP AP242 Standard
STEP AP242 Standard
3D 설계데이타 교환을 위한 국제표준 중립 포멧으로 STEP이 제정된지 수십년이 흘렀는데, 아직 CAD 분야에서는 STEP Native Modeler는 개발되지 않았다. 그리고 현재 현업에서는 대부분 AP203 및 AP214를 사용한다.우선 AP203에 관해서.
IGES를 대체하기 위해 STEP 표준이 논의되면서, 제일 먼저 나온 사양이 AP203이다. AP 명세를 구현하는 실체는, 텍스트로 표현된 Express 스키마 문법이다. 이때는 아직 XML이 나오기 전이므로, 전용의 구문 문법을 개발해서 사용한 것이다. AP203을 주도적으로 제정한 곳은 미국 국방분야의 항공회사들로 보통 알려져 있다. 즉 보잉 같은 회사를 말한다. 초창기 표준이므로 가벼운 대신 담고 있는 정보가 좀 부족하다. 예를 들어, 내부에 구성되어 있는 엔티티들의 색상을 정의하는 구문이 존재하지 않는다. 따라서 초창기 AP203 표준에 따라 구현된 STEP 파일 생성기로 만들어낸 STEP 파일을 CAD에서 읽어들이면, 색상정보가 다 빠져버려 탈색된 빨래 같은 것이 뜨게 된다.물론 시대가 흐르면서 AP203도 계속 업그레이드 되었다. 따라서 요즘 사용하는 버전의 3D CAD 툴에서 생성한 AP203 STEP 파일들은 색상 정보를 담아낼 수 있다. 그리고, 어셈블리 구조를 모듈화하여 표현한 업그레이드 버전인 AP203-ed2가 거의 최종 버전인 것 같은데… 이 표준으로 생성할 경우에는 여러 CAD 시스템 간에 데이타 교환시 실패 확률이 매우 높아진다. 서로 제대로 구현하지 못했다는 이야기다.
좋은 예로, CREO 2에서 AP203-ed2로 STEP 파일을 생성한 후에, 그 파일을 CATIA V5R20에서 읽어들이면 형상 재구성에 실패한다. 때문에 STEP 버전을 낮춰서, CREO 2에서 구버전의 AP203으로 생성한 후에, 그 파일을 CATIA V5R20에서 읽어들이면 성공적인 결과를 볼 수 있다. 이런 식의 미묘함 때문에, STEP을 이용한 이기종 CAD 간의 데이타 교환시에는 가능한 구버전의 STEP 표준으로 파일을 생성하는 것이 유리하다. 물론 그 대신, STEP 파일 안에 담기는 정보는 더 적어지게 되고, 실질적으로는 IGES 수준의 형상 정보 정도만 교환하는데 만족할 수 밖에 없다.
사실 AP203 표준에는 더 많은 정보를 담아내서 교환하는데 써야 하는게 맞는데, CAD 시스템에서 제대로 지원해 주지를 못한다. 예를 들어, 어떤 부품이 Made 인지 Bought 인지를 구분해 주는 아규먼트가 AP203의 Express 스키마 구문 중에 분명히 있지만(PRODUCT 관련 구문), CREO 이든 CATIA이든 간에 이 항목을 제대로 매치시켜 STEP으로 생성해 주지 못한다. 그냥 없는 항목처럼 무시하는 것이다.
따라서 이렇게 생성된 STEP 파일을 PDM 시스템에 업로드 해서 해석하여 데이타베이스화하고자 해도, 형상 정보 이외의 이런 부가 정보들은 현실적으로 활용되지 못하게 된다. 각 CAD 시스템이 STEP 표준을 이렇게 대충 지원하는 이유는, 일반적으로 자사 전용 포멧을 고객에게 강요함으로써 시장지배력을 올리고 또 고객을 자사 솔루션에 종속시키기 위한 것으로 생각된다.
그리고 AP214에 관해서.
한편, AP214는 독일쪽 자동차 업계에서 주도하여 2000년 경에 확정되었다. 미국에 정면 대항하려는 목적도 있지 않았을까. 당연히 AP203보다 더 발전된 사양으로, 방대한 명세를 포함한 거대한 포멧이 되었다. 그러나 이기종 CAD간 호환성은 확실하게 AP203보다 확연히 떨어진다. AP214로 생성한 파일을 협력업체에 보냈을 때, 높은 확률로 "이 파일 읽어지지가 않아요"라는 전화를 받곤 한다.내 경우엔, 이전 직장에서는 시제품 제작을 위해 목업 업체에 보낼 데이타는 AP214를 사용했다. 함께 서로 포멧을 조정한 기간을 거친 후에, 안정적으로 데이타 교환을 했기 때무에 문제가 없었다. 문제는 직장을 옮기고 다른 협력업체와 일을 하게 되자 그 새로운 회사에서 사용하는 소프트웨어가 달라졌기 때문에 자료 교환에 문제가 생긴 것이다. 결국 울며 겨자먹기로 AP203으로 포멧을 바꿀 수 밖에 없었다.
아무튼 AP214는 보다 더 많은 정보를 담아낼 수 있어 편리성은 높아졌지만, 이기종 CAD 간에 살짝 미묘하게 호환성이 떨어졌기 때문에 범용적으로 사용하기는 곤란한 점이 있었다.
AP242
마침내 AP203과 AP214를 합쳐서 AP242가 제정되었다. 그리고 그게 다시 업그레이드 되어 AP242-ed2 확정이 다 된 거나 마찬가지다. ed2로 올라가면서, 온갖 정보들을 더 많이 담게 되었다. PDM 관련한 각종 부가정보들, 부품간의 조립 구속조건, 파라메트릭 구속조건, 제조 관련 정보들, 전장 하네스 설계 정보, 히스토리 정보, 프리젠테이션을 위한 정보, 3차원 기구학 애니메이션을 위한 정보 등등. 이제까지의 전례를 봤을 때, 각 CAD 개발사들이 과연 AP242를 얼마나 적극적으로 준수해서 제대로 된 STEP 포멧 생성기를 만들어낼지 약간 의문이긴 하다. CAD 개발사들이 적극적으로 지원하지 않으면 사실상 죽은 표준이 될 공산도 크고…아무튼 PTC사의 경우에는, 이 자료를 신뢰한다면, CREO 4에 AP242 관련 기능들이 다 구현되어 제공될 것으로 보인다. 경쟁관계에 있는 CATIA 및 UG 쪽도 앞서거니 뒷서거니 하면서 비슷한 시기에 나오지 않을까 싶다.
한가지 가정을 해 보자.
CAD 시스템에서 올바르게 AP242 표준에 맞도록 데이타를 만들어 낸다면?
그러면 오랫동안 사용자들이 원해왔던, 그리고 CAD 개발사들이 그렇게 오랫동안 방해해왔던 이기종 CAD 간의 그나마 쓸만한 자료교환이 가능해 질 것이다.
그렇게 될 경우, 특정 CAD의 Native Format을 지원하는 모듈을 별도 개발할 필요 없이, AP242만 지원하는 오픈소스 모듈만 가지고도 필요한 PDM/PLM, CAD 및 기타 관련 소프트웨어들의 구현이 될 것이고, 이것만으로도 프로페셔널한 수준의 진정한 산업표준 데이타 교환이 가능해져 제조업 효율이 좋아질 것이다.
개인적으로, 특정 CAD 솔루션에 종속된 설계환경은 벗어나야 한다고 본다. 생각해 보자.
CREO를 능숙하게 사용하는 구조설계자가 있을 경우, 이 사람은 CATIA를 사용하는 업체에 취업할 수 없다. 이 얼마나 바보같은 일인가?
CAD 소프트웨어가 본질이 아니고, 공학적 설계능력 자체가 핵심인데 그전에 엉뚱한 허들이 형성되어 있는 것이다.
또, CATIA를 쓰는 업체는 다른 CAD로 전환하거나 다종의 CAD 시스템을 운용하기가 어렵다. 인력 운용면에서도 그렇고 시스템 측면에서도 그렇다. (물론 대규모 자동차 업계 같은 데서는 이기종 CAD간의 조합이 이루어지기는 하지만, 악영향이 최소화되도록 모듈별로 철저하게 구분된 방식이라고 생각된다.)
따라서 CATIA에 종속된 제조업체는, 높은 비용을 고정적으로 강요받게 되고, 자금운용에서 경직적인 부분이 발생한다. 이것은 경영상의 위협요인이 아닐 수 없다.
대한민국의 현실
아마 정년퇴임하셨을 것으로 추측되는 이상헌 교수라는 분이 80~90년대에 걸쳐 초창기 CAD 분야에서 업적을 꽤 쌓았다고 알고 있다. 당시 이 분야에서의 화두 중의 하나였던 비다양체 위상 기하학을 다룰 수 있는 구조체를 구현한 CAD커널이 한창 개발되고 실용화 단계에 접어들 때의 이야기다. 현재는 비다양체 CAD커널은 그냥 아주 당연한 기본 특성이 되어 있다.이상헌 교수의 경우, 자신의 연구에 펀딩만 좀 받았더라면 독자적인 CAD커널로 제대로 발전시킬 수 있었을 것이다. 왜냐면 같은 시기에 다른 나라들의 수준과 별 차이가 없었기 때문이다. 그러나 그런 장기적인 연구를 지원하는 시스템은 한국에 없었다. 사업가들의 눈에 띄지도 못했다. 이상헌 교수에게는 글로벌 경쟁에 나설 기회가 아예 주어지지 못했던 것이다.
그래서 지금의 한국은 그냥 CAD 및 PLM 분야에서 종속 국가가 되었다.
한국이 제조업 기반으로 흥한 나라라고들 하는데, 현재 상태로는 '제조업 2.0 (?)'이라던가 하는 말은 과대망상 내지는 허언증이다. 기초가 부실하다고 하는 좋은 예가 될 수 있을지 모르겠다.
2015년 6월 17일 수요일
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 이기 때문에, 매우 손쉽게 고수준 해석에 도전할 수 있는 환경이 되는 것이다.
2015년 6월 14일 일요일
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
총 시간은 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 결과를 좀 효율이 떨어지더라도 얻어낼 수 있다는 것.
2015년 6월 11일 목요일
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를 한 예가 거의 없는 것 같다. 아무리 검색해도 안 나온다.
마지막으로, 이상 검토한 조건을 기반으로, 충분히 더 잘게 쪼개어서 해상도를 높여서 돌려본 결과도 올려둔다..
2015년 6월 8일 월요일
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만불 목표금액을 채우는데 성공했나 보다.
이름하여 Vortex Bladeless.
별 건 아니고, 원형 단면 기둥을 하나 세워서 와류가 형성되면 기둥이 흔들거릴 테니깐, 그걸로 유도기전력을 발생시켜 전력을 생산한다는 거다.
당연히 효율은 좀 떨어지겠지만, 기존의 블레이드 방식의 경우 고가의 고정밀 동력전달기구 및 회전형 발전기가 들어가기 때문에 가격이 엄청 높아진다는 점을 공략했다.
고속 회전기구가 안 들어가므로 매우 싸게 만들수 있기 때문이다.
재료도 가볍고 싼 수지 계열로 했다.
또 회전 날개가 없으니깐 새들이 죽어나가는 일도 없다는 '명분'도 만들었다.
칼만 와류 모르는 기계공학도는 없겠지만, 이런 식으로 명분 만들어서 사업화 하는건 코페르니쿠스적 발상이라고 생각된다. 또 MIT라는 브랜드를 등에 업었으니 더 쉽게 되겠지.
이 제품의 디자인을 보면, 위로 갈수록 단면 지름이 커지도록 해 놨는데, 이는 지표면에서 위로 올라갈수록 점성에 의해 바람의 속도가 증가하게 되므로, 이를 고려하여 구배값을 결정한 것 같다.
또한, 세워진 원통 막대를 지지하는 카본 목 부위의 탄성 등을 감안하여, 막대의 고유진동수를 칼만 와류에 의해 가진되는 진동수와 잘 맞춰주도록 설계함으로써, 진동에 의해 발생하는 운동에너지를 극대화하려고 노력했을 것이다.
그리고, 유량계의 경우 사이즈가 작아서 진동이 초음파대로 나오도록 설계되어 있지만, 이 녀석의 경우에는 충분히 직경이 크기 때문에 인디고고에서 이 사람들이 밝힌 바에 따르면 음파는 20Hz 미만의 낮은 주파수대로 나온다고 한다. 따라서 새나 기타 동물, 인간에게 초음파에 의한 악영향이 없을 거라고 밝히고 있다.
2015년 6월 3일 수요일
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
자신의 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 작성
요즘 시간이 점점 쪼들리고 있는데, 신속하게 가능할지는 잘 모르겠다.
이상 끝!
2015년 5월 21일 목요일
GPG > DraftSight > FreeCAD > GMSH > ELMER > PARAVIEW 작업 흐름 강좌 (개정판)
초보자에게 간명하게 전체 작업 흐름을 파악할 수 있도록 하기 위해서 만든 자료이다.
작업 목표는 기어 이빨을 하나 설계해서, 그것의 굽힘응력을 보는 것이다.
기존의 예제와 달리 다음 몇가지 조건을 준수하였다.
1. Freeware 및 OpenSource 소프트웨어만을 사용한다. 유료 소프트웨어는 전혀 사용하지 않는다.
2. 기어 이빨은 실제 사용 가능한 수준의 엄밀성을 가지도록 설계한다.
(실제 인볼류트 곡선 함수 및 트로코이드 필렛을 적용한 모델을 설계 사용)
3. 매쉬 작업은 충분히 세밀하게 하고, 특히 관심부위를 조밀하게 생성한다.
4. 여러개의 프레임을 얻어내어 동영상을 만든다.
기어 이빨의 작도를 위해서는 FreeCAD의 옵션으로 들어있거나 제공되는 매크로를 사용해도 되지만, 검토해 보니까 정밀도가 낮고 제대로 된 인볼류트 함수를 사용한 것이 아니라서 제외하였다.
GPG, DraftSight로 기어 이빨 형상 설계
FreeCAD로 3D 모델링
GMSH로 매쉬 작업
ELMER로 해석
PARAVIEW로 결과확인 및 동영상 만들기
최종 결과물
100 프레임으로 다시 퀄리티를 높여 만든 동영상 결과물
피드 구독하기:
글 (Atom)