OpenSees를 다루기 시작하면서[코드 수정]

 오픈시스 (OpenSees: the Open System for Earthquake Engineering Simulation)

 이름에서도 알 수 있듯이 지진 응답 등을 시뮬레이션 할 수 있는 오픈소스이다. 캘리포니아 소재의 대학들에서 개발되어 PEER (Pacific Earthquake Engineer Center)에서 주도적으로 개발을 진행하였고 현재는 미국에서 지진연구에 사용하는 가장 일반적인 프로그램 중 하나가 되었다. 당장 Caltech, UC Berkeley, UCLA, Stanford University 등 캘리포니아 명문대학 출신의 엔지니어들과 학자는 이것을 활용하여 새로운 해석 툴 등의 자료를 공유하기에 최신의 연구 내용에 관심이 있다면 익혀둘 필요가 있겠다. 한국에서도 관련 전문가들 사이에서 점점 활용 빈도가 높아지고 있다.

 내가 해당 프로그램을 익혀야겠다 결심을하는 이유 중 하나는 국내에 압도적인 점유율을 가지고 있는 구조해석 상용프로그램인 MIDAS GEN의 경우 계산 모듈과 계산과정이 가려진채 입력값과 결과값 만을 알 수 있다보니 특수한 상황에서 내가 의도한대로 계산을 한 것인지 확인이 힘들다. 가령 구조적으로 강성차이가 매우 큰 부재가 연속된 위치의 경우 순간적인 힘에 대해서는 전단 지연등을 포함하여 힘의 전달이 용이하지 않기에 실제 강성보다 매우 작거나 혹은 무한대에 가까운 강성을 넣는 등의 작업을 하여 비교하기도 하는데, 컴퓨터의 경우 실수 계산시 부동소수점을 사용하므로 매우 큰수와 매우 작은 수를 같이 활용할 때 수학적 접근이 올바르지 않으면 오차가 생각 이상으로 커지기도 한다. 그러나 오픈시스를 사용하면 해석을 진행할때 수렴시키는 system 알고리즘을 자신이 직접 선택하기 때문에, 자신이 관련 참조 문헌과 논문을 읽고 방식의 차이를 이해 할 수 있다면 가장 정확한 방법을 선택가능하다. (대형건물 해석이나 수치해석을 진행한다면 알고리즘 취사선택의 영향력은 수십배의 성능 체감이 느껴질만큼 매우 강력하며 제대로 다루기는 매우 어렵다.)

 따라서 일반 실무자라도 보다 전문성을 가지려면 이러한 부분에 있어서 제대로 검토할 수학적 지식과 계산 과정을 직접 확인하거나 신뢰할 수 있는 해석 모듈을 다룰 수 있어야 하는데 그런 부분적 검토에서도 오픈시스는 매우 유용하다. 더욱이 해당 프로그램은 구조물의 지진해석을 위해 고안되었지만 범용성이 매우 뛰어나 국소적인 부재해석이 모두 가능하다. 커뮤니티에서 자료 등을 찾다보면 메쉬모델로 알류미늄 캔을 만들어 찌그러트리는 시뮬레이션을 진행한 것도 있는데 구조적 지식이 충분하다면 응용성이 상당하다.  (그렇다해도 지식의 깊이가 없다면 아예 아무것도 할 수 없는 프로그램이기에 연구자 출신이 아니라면 다루기 힘들긴하다... 나도 애를 먹고 있는 중이고...) 

 기존 상용해석 프로그램들과 달리 GUI를 활용한 pre,post 기능이 없으며 CLI (command line interface)를 활용하여 코딩만으로 구조물의 해석을 진행하게 된다. 물론 써드파티 프로그램으로 필요에 따라 관련 기능을 보완할 수 있지만, pre processing을 해당 프로그램으로 진행하면 노드 번호와 element 번호 등을 규칙에 맞추어 배정하여 문자 편집기를 통해 일괄적으로 구조물을 편집할 수 있는 코딩의 장점이 사라질 수 있다. post processing과 관련해서는 현재의 상용프로그램 중에서는 NextFEM 정도가 거의 유일하게 유용한듯 싶다. 

문자 편집기를 활용하여 오픈시스 사용

오픈시스 인터프리터 자체 실행시의 모습. 꽤나 업데이트가 자주 있으며 2020.08.28 현재는 2020.07.03에 배포된 3.2.2버젼이 최신이다.--- 2021년 4월 23일 까지도 버젼 유지되고 있는 상태인 것으로 보아 정식 알고리즘 추가 채택이 한동안은 없을지도 모르겠다.

 오픈시스는 TCL 언어를 활용한 구조해석용 인터프린터라 TCL 언어를 새로 익혀야 사용할 수 있는데 해당 언어는 C언어 다음으로 오랜된 프로그래밍 언어이면서 다른 언어와 문법적으로 차이가 많아 처음 익힐 때 다소 어색함이 있다. 그러나 언어 자체는 크게 어렵지 않은 편이라 금방 익힐 수 있었다. 다만 국내 프로그래머들은 거의 활용하지 않는 언어이다보니(활용 순위로보면 30~50위 사이에 있던 것으로 기억한다.) 실질적으로 미국 자료만으로 독학해야하는 불편함은 있다.

 (이와 관련해서는 파이썬용으로 새로 라이브러리가 2년 전에 개설된 것을 뒤늦게 확인하였는데, 기본 매뉴얼이 tcl기반이고 opensees wiki 등에서 제공되는 예제와 수십년의 레퍼런스 문헌들이 전부 tcl코드 기반으로 제공되는 프로그램이기 때문에 tcl 언어로 코딩하는게 여전히 편하게 느껴진다. 사실 매트랩으로 코드를 짜고 최종 입력 변수만 tcl언어로 변환시키는 과정을 스스로 선택하였기에 코딩언어 장벽은 크게 문제되지 않았다.)

 

 아래는 오픈시스를 처음 익히면서 점심시간에 기본 예제 중 하나를 다뤄본 것인데 철근 보강된 콘크리트에 힘을 작용할 시 곡률을 계산하는 모듈이다. 맨 아래 출력결과를 보면 철근의 항복점을 기반으로 한 추정 곡률은 계산이 되는데 실제 단면과 작용하중을 고려한 계산 모듈은 에러가 떠서 출력되지 않았다. 아무래도 내가 수학적으로 잘못 접근을 했거나 오픈시스 관련 문법에 미숙하여 생긴 문제이거나 둘 중 하나일텐데 지금은 오픈시스를 다루는 것 자체가 미숙하여 검토하는데 오래 걸릴듯하니 추후 천천히 확인을 해봐야겠다. [저녁식사 후 해당 내용 수정하여 아래에 첨부하였다.]

 
#proc MomentCurvature
#Arguments  
#    secTag-- tag identifying section to be analyzed
#    axialLoad--axial load applied to section (negative is compression)
#    maxK--maximum curvature reached during analysis
#    numlncr -- number of increments used to reach maxK (default 100)
# Sets up a recorder which writes moment-curvature results to file
#sectio$setTag.out... the moments is i coulumn 1, and curvature in column 2
proc MomentCurvature {secTag axialLoad maxK {numlncr 100} } {
    #Define two nodes at (0,0)
    node 1 0.0 0.0
    node 2 0.0 0.0
    #fix all degrees of freedom except axial and bending at node 2
    fix 1 1 1 1
    fix 2 0 1 0
    # Define element 
    #        tag ndi ndJ secTag
    element zeroLengthSection 1 1 2 $secTag
    #Create recorder
    recorder Node -file section$secTag.out -time -node 2 -dof 3 disp
    #Define constant axial load
    pattern Plain 1 "Constant" {
        load 2 $axialLoad 0.0 0.0
    }
    #Define analysis parameters
    integrator LoadControl 0 1 0 0
    system SparseGeneral -piv;
    test NormUnbalance 1.0e-9 10
    numberer Plain
    constraints Plain
    algorithm Newton
    analysis Static
    #Do one analysis for constant axial load
    analyze 1
    #Define reference moment 
    pattern Plain 2 "Linear" {
        load 2 0.0 0.0 1.0
    }
    #compute curvature increment
    set dK [expr $maxK/$numlncr]
    #Use displacement control at node 2 for section analysis
    #integrator DisplacementControl $node $dof $incr <$numIter $<math>\Delta U \text{min} </math> $<math>\Delta U \text{max}</math>>
    integrator DisplacementControl    2     3   $dK     1             $dK                               $dK
    #Do the section analysis
    analyze $numlncr
}

#example2.1 Moment Curvature analysis of a reinforced Concrete Section
#for the zero length element, a section discretized by concrete and steeel is created to represent the resultant behavior.
 
#units: kips,in,sec
#Define model builder
#--------------------------------
wipe
model BasicBuilder -ndm 2 -ndf 3
#Define materials for nonlinear coulumns
#--------------------------------
#Concrte             taG    f'c     ec0 f'cu ecu
#Core concrete (confined)
uniaxialMaterial Concrete01 1 -6.0 -0.004 -5.0 -0.014
#Cover concrete(unconfined)
uniaxialMaterial Concrete01 2 -5.0 -0.002 0.0 -0.006
#STEEL
#Reinforcing steel
set fy 60.0; #Yield stess
set E 30000.0; #Young's modulus
#        tag fy E0 b
uniaxialMaterial Steel01 3 $fy $E 0.01
#Define cross-section for nonlinear columns
#--------------------------------
#set some parameters
set colWidth 15
set colDepth 24
set cover 1.5
set As 0.60;    #area of no. 7bars
#some variables derived from the parameters
set y1 [expr $colDepth/2.0]
set z1 [expr $colWidth/2.0]
section Fiber 1 {
    #create the concrete core fibers
    patch rect 1 10 1 [expr $cover-$y1] [expr $cover-$z1] \
         [expr $y1-$cover] [expr $z1-$cover]
    #Create the concrete cover fibers(patch rect $mattag $numsubdivy $numsubdivZ $yi $zi $yj $zj)
    patch rect 2 10 1 [expr -$y1] [expr $z1-$cover]\
                       $y1 $z1
    patch rect 2 10 1 [expr -$y1] [expr $z1] $y1 [expr $cover-$z1]
    patch rect 2 2 1 [expr -$y1] [expr $z1-$cover] [expr $cover-$y1] [expr $z1-$cover]
    patch rect 2 2 1 [expr $y1-$cover] [expr $cover-$z1] $y1 [expr $z1-$cover]
    #Create the reinforcing fibers(left, middle, right)
    layer straight 3 3 $As [expr $y1-$cover] [expr $z1-$cover] [expr $y1-$cover] [expr $cover-$z1]
    layer straight 3 2 $As 0.0 [expr $z1-$cover] 0.0 [expr $cover-$z1]
    layer straight 3 3 $As [expr $cover-$y1]                [expr $z1-$cover]            [expr $cover-$y1]        [expr $cover-$z1]
}
#estimate yield curvate
#(Assuming no axial load and only top and bottom steel)
set d [expr $colDepth-$cover]    ;    #d--from cover to rebar
set epsy [expr $fy/$E]    ;#steel yield strain
set Ky [expr $epsy/(0.7*$d)]
#Print estimate to standard output
puts "Estimated yield curvature: $Ky"
#Set axial load
set P -180
set mu 15; #Target ductility for analysis
set numlncr 100; #number of analysis increments
#call the section analysis procedure
source MomentCurvature.tcl
MomentCurvature 1 $P [expr $Ky*$mu] $numlncr

결과값을 볼때 해석모듈을 잘못 작성한듯 싶지만 당장은 수정이 힘드니 추후에 확인을 해 보아야겠다.

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

앞서 코드에서 잘못된 부분을 내가 해당 언어를 잘 몰라서 잘못 입력했거니 했는데 저녁시간에 확인해보니 기하학적 좌표입력에서 부호를 잘못 고려하여 사각 튜브를 잘못 설정했기 때문이었다.

위의 코드 중

 

patch rect 2 10 1 [expr -$y1] [expr $z1] $y1 [expr $cover-$z1]

patch rect 2 2 1 [expr -$y1] [expr $z1-$cover] [expr $cover-$y1] [expr $z1-$cover]

 

부분이 다음과 같이 수정되면 정상적으로 계산이 됨을 확인하였다.

 

patch rect 2 10 1 [expr -$y1] [expr -$z1] $y1 [expr $cover-$z1]

patch rect 2 2 1 [expr -$y1] [expr $cover-$z1] [expr $cover-$y1] [expr $z1-$cover]

정상실행시 결과

실행시 output 값으로 생성되는 최종 해석값 [여기서 좌측값은 linear하게 하중을 넣는 것과 관련된 time 값이며 우측값이 최족적으로 나오는 곡률이다.

코드를 검토해준 Michael H. Scott, Ph.D. 에게 감사인사를 남긴다.

https://portwooddigital.com/

Portwood Digital