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 이기 때문에, 매우 손쉽게 고수준 해석에 도전할 수 있는 환경이 되는 것이다.