1. 서 론
유한요소법(Bathe et al., 1984)은 전산 상에 이산화(Discretization) 된 유한요소모델(Finite element model)을 통해 실제 구조물의 응답을 예측하며, 복잡한 구조시스템의 합리적인 근사치를 얻기 위해서는 일반적으로 많은 자유도(Degrees of freedom) 가 필요하다. 특히 선박 및 해양 구조물 해석은 고려해야 할 하중 Case가 상당히 많으며, 설계 최적화를 위한 재해석 과 정, 비선형 해석을 위한 반복 근사 과정까지 고려한다면 유 한요소 모델의 특성과 목적은 유지하되 자유도 크기를 줄여 한정된 전산 비용을 효율적으로 사용할 필요가 있다(Boo, 2019).
한편 실제 설계에서는 특정 부위에 대한 변위 및 응력 이 력 평가가 주를 이룬다는 점을 상기한다면, 전체 자유도의 결과를 도출하기보다는 관심 부(Local)의 결과만을 도출하는 모델 축소기법(Model reduction method)을 이용하는 것이 효율 적이다. 모델 축소기법은 모델의 전체 자유도를 각각 생략 된 부분(Omitted part)과 대상 부분(Target part)으로 행렬 순열 (Matrix permutation) 및 Partitioning 한 후 간단한 행렬 연산을 통해 자유도를 축소, 국부 해석한다(Boo and Oh, 2017). 이를 통해 전체 모델을 해석하지 않고도 관심 영역의 해석 결과 를 얻을 수 있으며 자유도의 개수를 줄임으로써 가장 직관 적이고 확실하게 전산 비용을 줄이는 방법이다. 정적 해석 을 위한 자유도 기반 모델 축소기법은 강성응축기법(Wilson, 1974)이 대표적이다. 강성응축기법은 Guyan(1965)의 질량, 강 성응축기법에서 관성 효과를 무시하고 강성만 응축함으로 써 정확한 정적 해석 솔루션을 얻을 수 있으며 Improved reduced system(IRS) 기법(O 'Callahan, 1989), Iterative 기법(Friswell et al., 1995;Xia and Lin, 2004), The component mode synthesis (CMS) 기법(Craig and Bampton, 1968;MacNeal, 1971)과 같은 여러 동적 모델축소기법의 초석으로도 사용되었다.
또한, 반복적인 유한요소모델 해석과정에서 설계 변경마 다 전체 모델의 강성 행렬과 힘 벡터를 구성하고 방정식을 풀어 변위를 구하는 과정을 되풀이하는 것은 비효율적이라 는 점을 개선하기 위해, 축소 모델을 사용한 구조 재해석 기 법(Boo, 2019)이 개발되었다. 축소 모델 재해석 기법은 유한 요소 모델을 설계 변경이 수시로 발생하는 부분(Target part) 과 축소되어 제거될 부분(Omitted part)으로 구분하여 자유도 를 축소, 설계 변경 발생 시 발생 부분만을 갱신(Updating)하 여 재해석을 수행한다.
본 연구는 유한요소법 기반의 정적 축소기법과 축소 모델 재해석 기법을 이용하여 반복 근사 과정으로 인해 많은 시 간과 전산 비용이 필요한 비선형 해석을 비선형 영역(Target part)만 갱신하여 효율적으로 수행하는 핵심 아이디어를 소 개하고, 제안된 기법을 수치 예제에 적용하여 축소 모델 비 선형 해석기법의 성능을 보여 주었다.
2. 선형 정적 해석
본 장에서는 유한요소법(Finite element method)을 이용한 선형 정적 해석 방정식의 유도과정을 서술한다. 유한요소 이론은 무한한 자유도(DOFs, Degrees of freedom)를 가진 연속 체 구조물을 이산화(Discretization) 과정을 통해 유한개의 자 유도를 갖는 유한요소(Finite element)로 모델링 하여 해를 구 하는 방법이다.
유한요소법의 변위(Displacement)와 변형률(Strain)은 연속 체 역학(Continuum mechanics)을 기반으로 정의되며, 응력 (Stress)은 외력이 작용할 때 구조물의 표면 또는 내부에서 힘의 평형(Equilibrium of force)을 만족하기 위해 작용하는 힘 으로, 물체의 Momentum balance를 만족한다. 이를 방정식으 로 나타내면 다음과 같이 표면과 내부의 힘의 평형에 관련 된 2개의 식으로 나타낼 수 있다.
여기서, 아래 첨자 i, j 는 텐서 형식(Tensor form)의 1st order와 2nd order를 의미한다. τ 는 응력 요소를 나타내며, fS 와 fB 는 각각 구조물의 표면 Sf 과 부피 V 에 적용된 하중 벡터(Force vector), n 은 접선 벡터(Tangent vector)를 의 미한다. 식(1)의 방정식들은 각각 탄성체 모델의 내부와 외 부 힘의 평형을 나타낸 방정식이다.
유한요소해석은 강체 운동을 피하고자 경계 조건 (Boundary condition)이 필요하며 이에 관한 식으로 물체 표면 변위 호환성(Compatibility) 식을 얻을 수 있다.
여기서, ui 는 변위 요소, Su 는 경계 조건 표면적, 는 경계 조건에서의 변위를 의미하며, 호환성 방정식의 경계 조건 변위는 항상 0이다.
마지막으로 Stress-Strain law에 의해 응력과 변형률 관계를 다음과 같이 나타낼 수 있다.
여기서, C 는 강성 텐서(Stiffness tensor), τ 는 응력 벡터 (Stress vector), ε 는 변형률 벡터(Strain vector)이다.
앞서 식(1), (2), (3)에서 얻은 방정식들은 유한요소해석의 해를 구하는 데 필요한 방정식들이며, 이를 강형(Strong form) 이라고 한다. 강형을 이용해 해를 구하면 자유도에 따라 수 많은 미지수를 갖는 연립방정식을 풀어야 하지만 가상 일 법칙(Principal of virtual work)을 이용한다면 이를 하나의 방정 식으로 풀 수 있다. 이는 다음과 같은 식으로 표현된다.
여기서, δu 는 가상변위(Virtual displacement), δε 는 가상 변형률(Virtual strain)이다. 식(4)의 좌 항은 내부 가상일, 우 항은 외부 가상 일을 나타내며, 내 외부 힘의 평형을 나타낸 다.
식(2)의 경계 조건을 통해 지배방정식을 만족하는 가상변 위와 실제 변위(Real displacement)를 도출하기 위한 내부 요 소 m 의 내부 m 변위 장(Internal displacement field)과 내부 변형률 장(Internal strain field)은 다음과 같이 표현할 수 있다.
여기서, H 와 B 는 각각 변위 보간 행렬(Displacement interpolation matrix)과 변형률 보간 행렬(Strain interpolation matrix), U 는 노드 변위 벡터(Nodal displacement vector)이다.
가상 일 법칙에서 응력은 식(3)의 Hooke’s law를 이용하여 강성 텐서와 변형률의 곱으로 나타내며 가상변위 벡터와 가 상 변형률 벡터 또한 식(5)에서 사용된 보간 행렬과 같은 행 렬을 사용하여 정의된다. 이를 이용해 식(4)의 좌 항과 우 항 은 각각 식(6)과 (7)로 표현할 수 있다.
여기서, 와 는 각각 내력 벡터(Internal force vector), 외력 벡터(External force vector)를 의미한다.
이때 내부 가상변위 장(Internal virtual displacement field)이 1의 물리량을 가진다고 가정하면 다음과 같이 정적평형 방 정식을 정의할 수 있다.
식(8)에서 K 는 강성 행렬을 의미하고, U 와 F 는 각각 변위 벡터와 힘 벡터를 의미한다.
3. 비선형 정적 해석
유한요소 모델에 외력이 가해졌을 때 발생하는 응력이 재 료의 항복점(Yield point)을 넘지 않을 경우, 응력과 변형률의 관계를 나타내는 그래프는 선형적으로 증가하며 이때 함수 의 기울기를 탄성계수(Young’s modulus)라고 한다. 반면에, 항 복점보다 강한 응력에 의해 영구변형이 발생하여 응력-변형 률 곡선의 기울기가 탄성 영역과 달라지는 구간을 소성 영 역이라고 한다.
본 연구에서는 비선형 유한요소 방정식을 풀기 위해 가 장 널리 사용되는 반복 근사 기법인 뉴턴-랩슨 기법 (Newton-Raphson method)을 사용한다(Bathe et al., 1984;Dvorkin and Bathe, 1984). 뉴턴 랩슨 반복은 소성 응력이 발생하였을 때, 내외부 힘이 평형할 때까지 내력과 변위를 반복하여 근 사한다(Friswell et al., 1995).
요소 응력에 대응하는 노드 힘 벡터(Nodal force vector)가 외력 벡터와 일치하는 유한요소 평형은 다음과 같이 나타낼 수 있다.
여기서, R*, F 는 각각 외력 벡터(External force vector), 내력 벡터(Internal force vector)이며 비선형 해석에서 내력 벡 터는 노드 변위(Node displacement)에 비선형적으로 의존하기 때문에 다음과 같은 반복 근사 과정이 필요하다.
여기서, 윗 첨자 i 는 반복 횟수를 뜻하며 식(10)을 변위– 힘 그래프에 대입하면 다음과 같은 변위증분(Displacement increment) 방정식을 얻을 수 있다.
여기서, K, ΔU 는 각각 접선 강성 행렬(Tangential stiffness matrix)과 증가 변위 벡터(Incremental displacement vector)를 의미하며 접선 강성 행렬은 힘 - 변위 그래프를 변 위에 대해 미분하여 아래와 같이 얻을 수 있다.
또한, 증가 변위 벡터는 반복 횟수가 증가함에 따라 계속 해서 갱신되어 변위 벡터에 더해지며
반복 근사의 극한에서는 식(9)의 유한요소 평형을 만족하 는 강성 행렬과 변위 벡터가 계산된다. 앞서 서술한 비선형 정적 해석의 뉴턴 랩슨 반복 과정은 아래 Fig. 1과 같다.
탄소성(Elastoplastic) 해석에서, 식(3)의 응력-변형률 관계는 다음과 같이 증분 형식으로 새롭게 나타낼 수 있다.
여기서, dεtotal 와 dεp 은 각각 총 변형률 증분(Total strain increment)과 소성 변형률 증분(Plastic strain increment)이며, dτ 는 응력 증분(Stress increment)을 의미한다.
4. 강성응축 기법을 활용한 비선형 해석
유한요소모델은 엔지니어의 의사결정에 의해 해석 대상 영역(Target part)과 축소되어 제거될 부분(Omitted part)으로 구분할 수 있으며, Node 번호를 이용한 간단한 행렬 순열 (Matrix permutation) 기법과 부분행렬 화(Matrix partitioning) 기 법을 이용하면 식(8)의 정적평형 방정식을 다음과 같이 새롭 게 표현할 수 있다.
여기서, 아래 첨자 o 와 t 는 각각 Omitted part와 Target part를 의미하며 ot 는 두 항의 결합 된 항을 의미한다.
이때 식(15)의 1, 2번째 행을 전개하면 다음과 같으며
여기서, 식(16)의 Ψ 는 Constraint mode matrix(Craig and Bampton, 1968;Craig and Hale, 1988)를 나타낸다.
이때 식(16), (17)를 식(15)에 대입하면 간단한 행렬 계산을 통해 다음과 같은 축소된 행렬과 벡터로 이루어진 방정식으 로 나타낼 수 있으며
마지막으로, 방정식의 0 행렬을 제거하면 아래와 같은 축소 정적평형 방정식을 얻을 수 있다.
앞서 강성응축기법을 통해 축소된 축소 정적평형 방정식 의 강성 행렬과 힘 벡터는 각각 다음과 같이 간소화하여 나 타낼 수 있다.
여기서, the over-bar ― 는 Target part를 제외한 나머지 축소 행렬 또는 벡터 항을 의미한다.
축소 모델 재해석 기법(Boo, 2019)에서는 축소 모델 행렬 이 Target part와 Omitted part의 간단한 덧셈으로 표현된다는 점을 이용한다. Fig. 2와 같이 Omitted part는 축소 후 경계 영 역(Boundary term)에 적용하여 재사용(Re-using)하며, 국부 해 석 영역은 갱신(Updating)되어 설계 변경이 적용된 새로운 축 소 모델을 구성한다.
따라서, 설계 변경이 적용된 축소 강성 행렬과 힘 벡터는 아래와 같이 표현할 수 있으며
식(19)의 축소 방정식은 아래와 같이 재정의된다.
여기서, 위 첨자 New 는 설계 변경이 적용되었음을 나타 낸다.
본 연구에서는, 앞서 설명한 축소 모델 재해석 절차가 축 소된 국부 모델의 강성과 힘을 갱신할 수 있다는 점을 이용 하여, 축소 모델의 강성 행렬과 힘 벡터를 갱신함으로써 반 복 근사 비선형 정적 해석을 수행한다.
정적축소된 유한요소 모델 해석에서 힘의 평형은 다음과 같은 방정식으로 나타낼 수 있다.
여기서, 는 각각 축소된 외력 벡터(Reduced external force vector)와 Target part의 내력 벡터(Internal force vector)이며, 식(11)의 변위증분 방정식은 다음과 같이 Target part의 변위증분 방정식으로 새롭게 표현할 수 있다.
접선 강성 행렬과 변위 벡터 또한 아래와 같이 Target part 에 대한 방정식으로 다시 얻을 수 있다.
앞서 3장과 마찬가지로, 축소 모델 반복 근사의 극한에서 는 식(23)의 유한요소 평형을 만족하는 강성 행렬과 변위 벡 터를 얻을 수 있으며, 강성응축 과정을 포함한 비선형 해석 절차는 Fig. 3과 같다.
본 연구에서는 선형 영역에 해당하는 모델을 축소하였기 때문에 모든 비선형 해석은 Target part에서만 발생한다. 따라 서, 축소 모델의 내부 요소 m 의 내부 변위 장과 내부 변형 률 장은 다음과 같이 표현할 수 있으며
여기서, 증분 형식의 응력-변형률 관계는 다음과 같다.
5. 수치 예제
제안된 방법의 성능을 검증하기 위한 수치 예제의 강성 행렬과 힘 벡터 모델링 및 구조 해석은 상용 FE 해석 software인 ABAQUS를 이용하였다. 유한요소 모델을 축소하 여 축소 모델을 구성하는 과정은 MATLAB 기반의 Static condensation code를 사용하였으며 계산된 축소 모델은 ABAQUS의 MATRIX INPUT 기능을 통해 모델의 경계 영역 에 적용하여 축소 모델 해석을 수행하였다.
본 연구에서는 해석 결과의 정확도 비교를 위해, 전체 모 델 해석 결과 변위 벡터와 축소 모델 해석 결과 변위 벡터의 L2-norm(Euclid norm)을 비교하였다. L2-norm의 비교는 행렬 크기 비교를 위해 널리 쓰이는 방법으로, 아래와 같이 행렬 벡터 내적 값을 비교함으로써 행렬 크기를 비교한다.
또한, 예제의 전체 모델 해석과 축소 모델 해석 시간을 측 정하여 비교함으로써 제안된 기법의 전산 시간효율을 검증 하였으며, 수치 예제의 축소비율을 제시하기 위해 수치 예 제로 사용한 Hatch corner plate 모델과 Barge 모델의 전체 DOFs 수와 축소된 모델의 DOFs 수를 Table 1에 나타내었다.
5.1 Hatch corner plate problem
본 예제는 아래 Fig. 4와 같이 B = 57.6 mm, H = 25.4 mm, 두 께 t = 2 mm인 평판에 반지름 3 mm의 아치와 길이 6.6 mm의 돌출부를 가지는 자유도 13050개의 Hatch corner plate를 Shell element로 모델링 하여 구조적 불연속 부에 국부 압력 하중 이 적용되는 정적 해석을 수행하였다.
또한, 유한요소 모델에 사용된 Shell 요소 타입은 ABAQUS 의 S4, 재질은 연철을 사용하였으며 물성값으로 탄성계수 E 는 206000 MPa, Poisson's Ratio는 0.3, 소성해석을 위한 Yield stress는 373.5 MPa로 설정하였다.
압력 하중은 돌출부에 선형 탄성 변형이 발생하는 10 MPa 과 소성 영역 변형이 발생하는 50 MPa을 각각 선형 하중 조 건, 비선형 하중 조건으로 구분하여 적용하였다. 하중이 적 용된 돌출부 전체를 Target part로 지정하였으며 경계 조건은 평판의 위아래 끝단에 적용된다. 선형 해석과 비선형 해석 의 결과 응력을 Figs. 5, 6에 전체 모델과 축소 모델을 함께 제시하였다.
본 예제는 선형 해석과 비선형 해석 모두 같은 모델을 사 용하였기 때문에 최초 모델축소 과정을 통해 제작된 축소 모델을 반복적으로 재사용하여 해석을 수행하였다. 전체 모 델과 축소 모델의 국부 영역(Target part) 해석 결과 변위를 X, Y, Z 자유도 별로 Figs. 7, 8에서 비교 하였으며 L2-norm 값 을 Table 2에서 검토하였다(Boo et al., 2016).
또한, 제안된 기법의 전산 효율을 검증하기 위해 전체 모 델의 해석 시간과 축소 모델의 해석 시간을 해석 조건별로 Table 3을 통해 각각 비교하였다.
5.2 Large barge problem
제안된 기법의 대형구조물 해석 효율을 검토하기 위해 아 래 Fig. 9와 같이 길이 L = 12.45 m, 폭 B = 3 m, 높이 H = 0.99 m 의 자유도는 362783개를 가지는 Barge 구조물을 Shell element 로 모델링 하여 정중앙 Deck에 국부 압력 하중이 적용되는 정적 해석을 수행한다.
또한, 유한요소 모델에 사용된 Shell 요소 타입은 ABAQUS 의 S4, 재질은 판 부재 알루미늄 5083H321, 보강재 알루미늄 6061T6을 사용하였으며 물성값 탄성계수 E = 68900 MPa, 70300 MPa을 적용하였으며 밀도 ρ 는 2660 kg/m3, 2700 kg/m3, Poisson's Ratio는 0.3을 각각 적용하였다. 또한, 소성해석을 위한 Yield stress는 241 MPa로 설정하였다. 압력 하중은 Barge의 Deck 중 앙 부에 0.2 MPa을 적용하였으며 경계 조건은 Barge의 양 끝 단에 적용하였다. 비선형 해석 결과 응력을 Fig. 10에 전체 모델과 축소 모델을 함께 제시하였다.
앞선 예제와 마찬가지로 전체 모델과 축소 모델의 국부 영역(Target part) 해석 결과 변위를 비교하였다. Fig. 11에서 X, Y, Z 자유도 별로 비교하였으며 L2-norm 값은 Table 4에서 검토하였다.
또한, 제안된 기법의 전산 시간효율을 검증하기 위해 전 체 모델의 비선형 정적 해석 시간과 축소 모델의 비선형 정 적 해석 시간을 Table 5를 통해 비교하였다.
6. 결 론
본 연구에서는 강성응축기법을 활용하여 비선형 해석을 효율적으로 수행하는 해석기법을 제안하였다. 제안된 기법 의 효율성을 검증하기 위해 Hatch corner plate 구조물과 대형 Barge 구조물을 정적 축소하여 축소 모델 선형, 비선형 해석 을 수행하고 전체 모델 해석 결과와 비교하여 검토하였다.
축소 모델 국부 정적 선형 해석의 경우 전체 모델 해석보 다 2배가량, 비선형 해석의 경우 최대 23배가량 빠르게 계산 할 수 있었다. 특히 비선형 해석은 해석에 필요한 정량적인 시간이 훨씬 크기 때문에 더 효율적이며 축소 과정에 근사 과정이 없어 해석의 정확도가 보장된다는 점도 중요하다.
이 연구에서 제안한 축소 모델 해석기법을 바탕으로 기존 Static condensation 기법을 이용한 축소를 대형구조물의 정적 축소에 특화된 축소기법인 Algebraic multi-level substructuring 기반의 Automated static condensation 기법으로 개선하여 축소 시간을 단축하고, 사용성 개선 및 축소 시간 단축을 위해 FORTRAN으로 구성된 ABAQUS의 Subroutine 모듈화 관련 연 구가 필요하다.
마지막으로, 엔지니어링 실무에서의 제안된 기법 성능 검 증을 위해 수백만 개 이상의 자유도를 가지는 실제 구조물 을 축소 모델 해석하여 결과를 비교할 필요가 있다.