Isogeometric Analysis for Two-dimensional Multipatch Model

Min-Geun Kim; Bonyong Koo

doi:10.7734/COSEIK.2017.30.6.515

Preview

Journal of the Computational Structural Engineering Institute of Korea. 2017. 515-522
https://doi.org/10.7734/COSEIK.2017.30.6.515

Isogeometric Analysis for Two-dimensional Multipatch Model

2차원 멀티패치 모델의 아이소-지오메트릭 해석

Min-Geun Kim¹

Bonyong Koo²^†

김 민 근¹

구 본 용²^†

¹Korea Institute of Machinery & Materials(KIMM), Daejeon, 34103, Korea

²School of Mechanical Convergence System Engineering, Kunsan National Univ., Kunsan, 54150, Korea

¹한국기계연구원 스마트기계연구실

²군산대학교 기계융합시스템공학부

^†Corresponding author : +82-63-469-4713; bykoo@kunsan.ac.kr

License (open-access, http://creativecommons.org/licenses/by-nc/3.0/):

This is an Open-Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/3.0) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.

ABSTRACT

In this paper, an isogeometric analysis for multipatch problem is investigated, in which two or more geometries are connected at the interface in a conforming or non-conforming conditions. To express higher continuity at the patch interface, two approaches such as Nitsche based method and master-slave method are formulated for the linear elasticity problem and discretized using the isogeometric approach using NURBS basis functions. A short comparison between two approaches in formulations reveals the pros and cons of them with the applicability in the isogeometric multipatch problem. In addition, a NURBS based stress recovery is adopted to express a better stress continuity through the post-processing. Numerical examples indicate the effectiveness of Nitsche method in the non-conforming patch, following the exact solution well. For the stress concentration problem with the conforming patch, introduced two methodologies show comparative results, meanwhile the NURBS based stress recovery presents an improved smooth stress contour in the whole domain including the patch interface.

Keywords

isogeometric analysis

multipatch

nitsche method

stress continuity

본 논문에서는 2개 이상의 기하형상이 순응 또는 비순응 경계면에서 접합된 멀티패치 문제에 대한 아이소-지오메트릭 해 석에 대해서 연구하였다. 패치 경계면에서 응력의 연속성을 표현하는 방법으로 Nitsche 방법론과 마스터-슬레이브 방법에 기반한 방법론에 대해서 지배방정식을 유도하고 아이소-지오메트릭 이산화를 수행하였다. 멀티패치 문제에 대해서 두 방법 론의 차이점을 간단하게 비교하였으며, 후처리 과정에서 사용되는 NURBS 곡면 기반의 응력 복원법에 대해서 기술하였다. 수치예제에서 비순응 경계면을 가지는 멀티패치 빔 문제를 통해 Nitshce 방법론을 검증하였으며, 응력집중을 가지는 문제에 서 소개된 두 방법론이 유사한 결과를 보이는 것을 확인하였다. 소개된 NURBS 곡면 기반의 응력 복원법을 후처리에서 도 입할 경우 멀티패치 문제의 경계면에서 개선된 연속적인 응력을 보임을 알 수 있다.

키워드

아이소-지오메트릭 해석

멀티패치

Nitsche 방법론

응력의 연속적 표현

MAIN

1. 서 론

CAD 모델은 곡면의 기하학적 형상을 수치적으로 표현하는 방법으로 20세기 후반컴퓨터의 발전과 함께활성화되어 자동차, 항공, 조선을 포함한 전 산업분야에서 널리 사용되고 있으며, 반면 경계조건을 가지는 문제에 대해서 지배 미분방정식을 이산화(discretization)하여 수치해를 구하는 유한요소해석은 CAD 모델과 독립적으로 개발되어 구조해석, 열전달과 같은 문제에 대해 표준 해석방법으로 사용되고 있다. CAD 모델은 정확한 경계를 표현하기 위해 고차원의 기저함수를 사용하며, 유한요소해석의 경우 1차 기저함수를 사용하여 경계와 내부 요소를 표현하는 것이 일반적이다.

경계의 기하학적 형상을 표현하는 CAD 모델과 경계조건에 지배되는 지배방정식의 해를 구하는 유한요소해석의 공통점 으로 인해 CAD 모델과 해석을 통합하는 다양한 연구들이 수행되어 왔으며, CAD 모델의 형상을 표현하는 NURBS (Non-Uniform Rational B-Spline)를 활용하여 해석공간의 기저함수를 구성하는 “아이소-지오메트릭 해석(IGA)”이 제안 되었다(Hughes et al., 2005). 아이소-지오메트릭 해석은 CAD 모델 구성과 동일한 이산화가 가능하여, 경계에서 법선 벡터, 곡률 등의 기하 정보를 해석공간에서 직접 고려할 수 있는 장점을 가진다(Ha et al., 2012). 아이소-지오메트릭 해석 모델의 최소 단위는 1개의 곡면을 표현하는 패치(patch)로서 패치 단위로 갤러킨(galerkin) 정식화와 이산화를 수행하게 된다.

일반적으로 곡선 또는 곡면의 접합은 CAD에서 두 개 이상의 패치로 구성된 멀티패치(multipatch)로 이루어지며, 패치 접합면(interface)의 표현 파라메트릭 변수가 동일하지 않으 므로 연속성이 감소하게 된다. 이러한 멀티패치 곡선 또는 곡면의 연속성을 나타내기 위해서 파라메트릭 연속, Cⁿ과 기하 연속, Gⁿ의 개념을 도입할 수 있다. 파라메트릭 변수로 정의된 곡선 또는 곡면이 변수에 대해 단순 연속인 경우 C⁰, 변수의 1차 미분이 연속인 경우 C¹으로 표현할 수 있으며, 반면 기하 연속의 경우 기하학적으로 곡선 또는 곡면이 단순 연결된 경우 G⁰,접선또는 접하는면을 공유할경우G¹으로 정의할수 있다. 일반적으로 Cⁿ인 곡면은 Gⁿ임이 알려져 있으며, CAD 모델 에서는 곡면의 부드러운 연결을 위해서 C¹ 이상의 고차 기저 함수를 사용하게 된다.

멀티패치 곡면을 가지는 CAD 모델에 대해서 아이소-지오 메트릭 이산화를 수행할 경우 패치 간 변위에 대해서 연속인 C⁰인 해를 얻게 된다. 그러나C¹ 이상의 연속성을 가지는 고차 기저함수를 사용하였으므로, 멀티패치 간 변위의 1차 미분값의 조합, 즉 응력의차이를 최소화하여 표현할수 있게 되며, 이러한 멀티패치를 가지는 아이소-지오메트릭 해석에 대한 여러 연구가 수행되어 왔다. 패치 간 NURBS 곡면 조정점이 일치하는 경우 (conforming patch)와 일치하는 않는 경우(non-conforming patch)에대해서패치간마스터-슬레이브성질(Cottrell et al., 2007), 라그랑지 곱수(Brivadis et al., 2015), Nitsche 방법론(Nguyen et al., 2015) 등을 적용하여 멀티패치 사이 에서 해의 연속성을 표현하고 있다.

본 논문에서는 이러한 멀티패치 간 연속성을 표현하기 위한 방법론 중 하나인 Nitsche 방법론에 대한 멀티패치 지배 방정 식의 유도와 이산화 및 마스터-슬레이브 방법론을 소개하고, 수치예제로 2차원 멀티패치 문제에 적용하여 Nitsche 방법론 과 마스터-슬레이브 성질을 이용한 방법론의 결과를 비교 검토 하려고 한다. 추가적으로 NURBS 근사에 기반한 응력 복원 (Hassni et al., 2012) 방법을 도입하여 후처리 단계에서 멀티패치 간 응력의 차이를 최소화하여 표현하는 연구를 수행 하고자 한다.

2. 본 론

이 장에서는 NURBS 기저함수를 소개한 다음, 멀티패치 문제에 대해 Nitsche 방법론과 마스터-슬레이브 방법론을 적용한 지배방정식의 전개와 아이소-지오메트릭 이산화를 수행하고자 한다.

2.1. NURBS 기저함수

CAD 모델과 아이소 지오메트릭 해석의 해 공간 기저함 수로 사용되는 NURBS를 정의하기 위해서 재귀적으로 구성된 B-Spline을 먼저정의한다. B-Spline 곡선은 1차원 B-Spline $ϕ_{i, p} (ξ)$ 에 조정점(control point) B_i을 곱한 형태로 표현되며, p차의 기저함수가 주어졌을 때 아래와 같이 표현된다.(2)

(1)

C (ξ) = \sum_{i = 1}^{n} ϕ_{i, p} (ξ) B_{i}

여기서,(2)(3)

(2)

ϕ_{i, 0} (ξ) = {\begin{array}{l} 1 & if ξ_{i} \leq ξ \leq ξ_{i + 1} \\ 0 & o t h e r w i s e \end{array} (p = 0)

(3)

\begin{array}{l} ϕ_{i, p} (ξ) = \frac{ξ - ξ_{i}}{ξ_{i + p} - ξ_{i}} ϕ_{i, p - 1} (ξ) + \frac{ξ_{i + p + 1} - ξ}{ξ_{i + p + 1} - ξ_{i + 1}} ϕ_{i + 1, p - 1} (ξ) \\ p = 1, 2, 3, \dots, n \end{array}

B-spline 곡면과 3차원의 입체는 곡선과 동일한 형태로 아래와 같이 표현된다.(4)(5)

(4)

S (ξ, η) = \sum_{i = 1}^{n} \sum_{j = 1}^{m} ϕ_{i, p} (ξ) ϕ_{j, q} (η) B_{i, j}

(5)

V (ξ, η, ζ) = \sum_{i = 1}^{n} \sum_{j = 1}^{m} \sum_{k = 1}^{l} ϕ_{i, p} (ξ) ϕ_{j, q} (η) ϕ_{k, r} (ζ) B_{i, j, k}

B-spline 함수에 가중치(weight)를 곱하여 분수 형태로 나타낸 기저함수를 NURBS 라고 정의하며 아래와 같이 표현 한다.(6)

(6)

N_{i, p} (ξ) = \frac{ϕ_{i, p} (ξ) w_{i}}{\sum_{i = 1}^{p} ϕ_{i, p} (ξ) w_{i}}

B-spline 곡선과 마찬가지로 n개의 p차 기저함수 N_i,p와 대응되는 조정점 B_i가 주어졌을 때, NURBS 곡선은 다음과 같이 정의된다.(7)

(7)

C (ξ) = \sum_{i = 1}^{n} N_{i, p} (ξ) B_{i}

또한 NURBS 곡선의 텐서 곱으로 NURBS 평면을 표현할 수 있으며, 이 때 NURBS 평면의 파라메트릭 노트벡터는 $ξ = {ξ_{1}, ξ_{2}, \dots, ξ_{n + p + 1}}, η = {η_{1}, η_{2}, \dots, η_{l + q + 1}}$ 이다.(8)

(8)

S (ξ, η) = \sum_{i = 1}^{n} \sum_{j = 1}^{m} N_{i, p} (ξ) N_{j, p} (η) B_{i, j}

2.2. 선형 탄성 문제의 멀티패치 지배방정식

선형 탄성문제에 대한 멀티패치의 정의와 경계면(interface) 을 가지는 문제에 대해 Nitsche 방법론을 도입한 지배 방정식을 기술하고자 한다. Fig. 1에서 경계(boundary) $Γ \equiv \partial Ω$ 를 가지는 영역(domain), Ω에 대해 서로 겹치지 않는 영역 Ω^m, m=1,2를 정의하고 경계에 주어지는 디리클레(dirichlet) 경 계, $Γ_{d}^{m}$ 와 노이만(neumann) 경계, $Γ_{t}^{m}$ 그리고 패치 간 경계 면(interface) $Γ_{i}$ 를 정의하였다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2017-030-06/07COSEIK-30-515/images/COSEIK-30-515_F1.jpg

Fig. 1

Linear elasticity problem with the interface

Fig. 1에서 응력 $σ_{i, j, j}^{m}$ , 체적력 $b_{i}^{m}$ , 변위 $z_{i}^{m}$ , 디리클레 경계 조건 $g_{i}^{m}$ , 노이만 경계조건 $h_{i}^{m}$ 을 정의하면 주어진 선형 탄성 문제의 지배방정식은 다음과 같다.(9a)(9b)(9c)(9d)(9e)

(9a)

\begin{array}{l} - σ_{i j, j}^{m} = b_{i}^{m} & on Ω^{m} \end{array}

(9b)

\begin{array}{l} z_{i}^{m} = g_{i}^{m} & on Γ_{d}^{m} \end{array}

(9c)

\begin{array}{l} σ_{i j}^{m} n_{j}^{m} = h^{m} & on Γ_{t}^{m} \end{array}

(9d)

\begin{array}{l} z_{i}^{1} = z_{i}^{2} & on Γ_{I} \end{array}

(9e)

\begin{array}{l} σ_{i j}^{1} n^{1} + σ_{i j}^{2} n^{2} = 0 & on Γ_{I} \end{array}

H¹ 소볼레프 공간(sobolev space)에 속하는 해공간 (solution space) Z와 시험함수 공간(trial function space) Z 를 정의하고 가상일의 원리를 적용하면 약형(weak fom)으로 표현된 지배방정식 (10a)과 부분적분을 수행한 식 (10b)를 얻게 된다.

(10a)

\int_{Ω} \bar{z_{i}} (σ_{i j, j} (z) + b_{i}) d Ω = 0

(10b)

\int_{Ω} \bar{z_{i}} σ_{i j, j} (z) = \int_{Ω} {(\bar{z_{i}} σ_{i j} (z))}_{, j} d Ω - \int_{Ω} {(\bar{z_{i, j}} σ_{i j} (z))}_{, j} d Ω

식 (10b)의 우변 첫 번째 항을 경계면에서 변위는 일치하고 수직 벡터는 반대이며 합력이 0인 것을 이용하여 다음과 같이 나타낼 수 있다.

(11)

\begin{array}{l} \int_{Ω} {({\bar{z}}_{i} σ_{i j} (z))}_{, j} d Ω = \\ \int_{Γ_{t}} \bar{z_{i}} t_{i} d Γ + \int_{Γ_{I}} ‖ \bar{z_{i}} σ_{i j} (z) ‖ n_{j}^{1} d Γ \end{array}

여기서, 경계면 차이에 대한 표현식 $‖ \cdot ‖ = {(\cdot)}^{1} - {(\cdot)}^{2}$ 과, 평균에 대한 표현식 $〈 \cdot 〉 = 1 / 2 {{(\cdot)}^{1} - {(\cdot)}^{2}}$ 을 이용하여 식 (11)의 우변 두 번째 항은 다음과 같이 표현 가능하다.

(12)

‖ \bar{z_{i}} σ_{i j} (z) ‖ = ‖ \bar{z_{i}} ‖ 〈 σ_{i j} (z) 〉 + 〈 \bar{z_{i}} 〉 ‖ σ_{i j} (z) ‖

경계면에서 $‖ {\bar{z}}_{i} ‖ = 0$ 것을활용하여정리된식(11), (12)을 식 (10b)에 대입하면 다음과 같은 패치 간 경계의 응력 조건을 포함하는 약형(weak form)의 지배방정식을 얻게 된다.(13)

(13)

\begin{array}{l} \int_{Ω} {\bar{z}}_{i, j} σ_{i j} (z) d Ω - \int_{Γ_{I}} 〈 {\bar{z}}_{i} 〉 ‖ σ_{i j} (z) ‖ n_{j}^{1} d Γ = \\ \int_{Ω} {\bar{z}}_{i} b_{i} d Ω + \int_{Γ_{t}} \bar{z_{i}} h_{i} d Γ \end{array}

또한 경계면에서 적합성 조건 $z_{i}^{1} = z_{i}^{2}$ 과 Nitsche 방법론의 안정화 계수(stability parameter) α와 대칭성을 부여하기 위해 적합성 조건 $‖ z_{i} ‖ = 0$ 을 다시 고려하면 다음과 같은 최종 형태의 적분 지배방정식으로 표현 가능하다.

(14)

\begin{array}{l} \int_{Ω} {\bar{z}}_{i, j} σ_{i j} (z) d Ω - \int_{Γ_{I}} 〈 {\bar{z}}_{i} 〉 ‖ σ_{i j} (z) ‖ n_{j}^{1} d Γ = \\ - \int_{Γ_{t}} ‖ \bar{z_{i}} ‖ 〈 σ_{i j} (\bar{z_{i}}) 〉 n_{j}^{1} d Γ + \int_{Γ_{I}} α ‖ \bar{z_{i}} ‖ ‖ z_{i} ‖ d Γ \\ = \int_{Ω} {\bar{z}}_{i} b_{i} d Ω + \int_{Γ_{I}} \bar{z_{i}} h_{i} d Γ \end{array}

여기서, 안정화계수 α는해의정확도와수치적안정성을위해서 아래 식의 C₁의 제곱보다 크도록 정해야 하며, 일반적으로 $α = 2 C_{1}^{2}$ 으로 정할 수 있다(Dolbow and Harari, 2009).(15)

(15)

{‖ \nabla \bar{z} \cdot n ‖}_{L^{2}, Γ_{I}} \leq C_{1} {‖ \bar{z} ‖}_{E, Ω}

이러한 멀티패치 Nitsche 방법론의 상세한 유도 과정은 관련 참고문헌(Sanders et al., 2011)에서 찾을 수 있다.

2.3. Nitsche 방법론 기반 아이소-지오메트릭 이산화와 NURBS 근사기반 응력 복원

탄성 문제의 해석에 사용되는 형상 변수, x^m과 변위, z^m를 NURBS 기저 함수 N_k와 형상 변수 벡터, $D_{k}^{m}$ , 변위 조정점 벡터, $y_{k}^{m}$ 로 이산화하여 나타내면 다음과 같다.(16a)(16b)

(16a)

x^{m} (Ξ) = \sum_{k} N_{k}^{m} (Ξ) D_{k}^{m}

(16b)

z^{m} (Ξ) = \sum_{k} N_{k}^{m} (Ξ) y_{k}^{m}

이산화된 형상 변수와 변위를 멀티패치 문제의 지배방정식 에 대입하면 식 (14)의 좌변 항은 아래와 같이 표현된다.

(17a)

\begin{array}{l} \int_{Ω} {\bar{z}}_{i, j} σ_{i j} (z) d Ω = \\ \sum_{K, J} {\begin{matrix} y_{K}^{1} \\ y_{K}^{2} \end{matrix}}^{T} [\begin{matrix} \int_{Ω} B_{K}^{1T} C^{1} B_{J}^{1} d Ω \int_{Ω} B_{K}^{1T} C^{2} B_{J}^{2} d Ω \\ \int_{Ω} B_{K}^{2T} C^{1} B_{J}^{1} d Ω \int_{Ω} B_{K}^{2T} C^{2} B_{J}^{2} d Ω \end{matrix}] {\begin{matrix} y_{J}^{1} \\ y_{J}^{2} \end{matrix}} \end{array}

(17b)

\begin{array}{l} \int_{Γ_{I}} ‖ {\bar{z}}_{i} ‖ 〈 σ_{i j} (z) 〉 n_{j}^{1} d Ω = \\ \frac{1}{2} \sum_{K, J} {\begin{matrix} y_{K}^{1} \\ y_{K}^{2} \end{matrix}}^{T} [\begin{matrix} \int_{Γ_{I}} N_{K}^{1T} {nC}^{1} B_{J}^{1} d Γ \int_{Γ_{I}} N_{K}^{1T} {nC}^{2} B_{J}^{2} d Γ \\ - \int_{Γ_{I}} N_{K}^{2T} {nC}^{1} B_{J}^{1} d Γ \int_{Γ_{I}} N_{K}^{2T} {nC}^{2} B_{J}^{2} d Γ \end{matrix}] {\begin{matrix} y_{J}^{1} \\ y_{J}^{2} \end{matrix}} \end{array}

(17c)

\begin{array}{l} \int_{Γ_{I}} ‖ {\bar{z}}_{i} ‖ 〈 σ_{i j} (z) 〉 n_{j}^{1} d Ω = \\ \frac{1}{2} \sum_{K, J} {\begin{matrix} y_{K}^{1} \\ y_{K}^{2} \end{matrix}}^{T} [\begin{matrix} \int_{Γ_{I}} B_{K}^{1T} C^{1 T} n^{T} N_{J}^{1} d Γ - \int_{Γ_{I}} B_{K}^{1T} C^{2T} n^{T} N_{J}^{2} d Γ \\ \int_{Γ_{I}} B_{K}^{2T} C^{1 T} n^{T} N_{J}^{1} d Γ - \int_{Γ_{I}} B_{K}^{2T} C^{2T} n^{T} N_{J}^{2} d Γ \end{matrix}] {\begin{matrix} y_{J}^{1} \\ y_{J}^{2} \end{matrix}} \end{array}

(17d)

\begin{array}{l} \int_{Ω} α ‖ {\bar{z}}_{i} ‖ ‖ {\bar{z}}_{i} ‖ d Γ = \\ \sum_{K, J} {\begin{matrix} y_{K}^{1} \\ y_{K}^{2} \end{matrix}}^{T} [\begin{matrix} \int_{Γ_{I}} α N_{K}^{1T} N_{J}^{1} d Ω - \int_{Γ_{I}} α N_{K}^{1T} N_{J}^{1} d Ω \\ - \int_{Γ_{I}} α N_{K}^{1T} N_{J}^{1} d Ω - \int_{Γ_{I}} α N_{K}^{1T} N_{J}^{1} d Ω \end{matrix}] {\begin{matrix} y_{J}^{1} \\ y_{J}^{2} \end{matrix}} \end{array}

여기서 C는 탄성행렬, B는 변위-변형률 행렬이며 n은 패치간 경계면에서의 법선 벡터를 나타낸다. 위의 이산화 형태 의식을 지배방정식에 대입하여 정리하면 다음과 같은 최종적인 이산화된 행렬 방정식을 얻게 된다.

(18)

\sum_{J} [K_{K J}^{b} + K_{K J}^{n} + K_{K J}^{n^{T}} + K_{K J}^{s}] {\begin{matrix} y_{J}^{1} \\ y_{J}^{2} \end{matrix}} = {\begin{matrix} F_{K}^{1} \\ F_{K}^{2} \end{matrix}}

여기서, 행렬 $K_{K J}^{b}, K_{K J}^{n}, K_{K J}^{n^{T}}, K_{K J}^{s}$ 는 각각 식 (17a), (17b), (17c), (17d)의 기저함수로부터 정의된 행렬을 나타낸다. 우변의 하중 벡터는 일반적인 아이소-지오메트릭 이산화와 동일하며 다음과 같이 정의된다.(19)

(19)

F_{K}^{m} = \int_{Ω^{m}} N_{K}^{mT} b^{m} d Ω + \int_{Γ_{t}^{m}} N_{K}^{mT} h^{m} d Γ

식 (18)로부터 변위에 대한 조정점의 해 y를 얻게 되며 응력은 응력과 변위의 관계식 σ = C B y에서 계산되게 된다. 일반적으로응력값은가우스적분점(gauss quadrature point) 에서 가장 정확한 것으로 알려져 있으며 NURBS 기반의 응력 복원 곡면을 다음과 같이 정의할 수 있다.

(20)

σ^{*} (Ξ) = \sum_{k} N_{k} (Ξ) P_{k}

여기서, P_k는 복원 응력 곡면을 표현하는 조정점으로서 식 (20)와 각 가우스 적분점에서 얻어진 응력의 차이를 최소화하는 식 (21)에서 결정 가능하다.

(21)

\frac{\partial (\sum_{k} {(σ^{*} (Ξ) - σ)}^{2})}{\partial P_{k}} = 0

2.4. 마스터-슬레이브 방법론을 이용한 멀티패치 아이소 -지오메트릭 이산화

아이소-지오메트릭 해석에서 멀티패치 사이의 연결을 처음 으로 다룬 마스터-슬레이브 방법론(Cottrell et al., 2007)을 간단하게 설명하고자 한다. 패치 1과 패치 2에서의 변위 조정점 y^m을 내부 변위 조정점 $y_{n}^{m}$ 과 경계 변위 조정점 $y_{i}^{m}$ 로 분리하여 표현하면 식 (22)와 같다.

(22)

y^{1} = [\begin{matrix} y_{n}^{1} \\ y_{i}^{1} \end{matrix}], y^{2} = [\begin{matrix} y_{n}^{2} \\ y_{i}^{2} \end{matrix}]

세분화(refinement) 전 경계면에서 패치간 동일한 변위 조정점을 적용하므로 $y_{i}^{1} {= y}_{i}^{2}$ 이며 h-세분화를 위한 Oslo 알고 리즘을 적용하여 세분화된 패치 2의 변위 조정점 $y_{h}^{2}$ 를 표현하면 다음과 같다.

(23)

y_{h}^{2} = [\begin{matrix} T_{n} \\ T_{i} \end{matrix}] [\begin{matrix} y_{n}^{2} \\ y_{i}^{2} \end{matrix}]

기존 선형탄성 문제를 기술하는 식 (10)에서 이산화를 수행 할 경우 멀티패치 문제에 대해서 식 (24)로 표현 가능하며 식 (23)을 대입하면 식 (25)로 표현된다. 식 (25)로 사용하여 이산화된 멀티패치의 해를 구하는 것을 마스터-슬레이브 방법 이라고 하며, 초기에 멀티패치 사이의 변위 조건을 표현하는 $y_{i}^{1} {= y}_{i}^{2}$ 식이 필요하므로 이러한 조건을 만족하는 순응 경계면 (conforming interface)에 대해 적용할 수 있다.

(24)

[\begin{matrix} \begin{matrix} K_{n n}^{1} K_{n i}^{1} \\ K_{n i}^{1} K_{i i}^{1} \end{matrix} \\ \begin{matrix} K_{n i}^{2} K_{n i}^{2} \\ K_{n i}^{2} K_{n i}^{2} \end{matrix} \end{matrix}] [\begin{matrix} y_{n}^{1} \\ y_{i}^{2} \\ y_{n}^{2} \\ y_{i}^{2} \end{matrix}] = [\begin{matrix} f_{n}^{1} \\ f_{i}^{2} \\ f_{n}^{2} \\ f_{i}^{2} \end{matrix}]

(25)

[\begin{matrix} \begin{matrix} K_{n n}^{1} & K_{n i}^{1} \end{matrix} \\ K_{i n}^{1} (K_{i i}^{1} + T_{i}^{T} + K_{i i}^{2}) T_{i}^{T} K_{i n}^{2} \\ \begin{matrix} K_{n i}^{2} + T & K_{n n}^{2} \end{matrix} \end{matrix}] [\begin{matrix} y_{n}^{1} \\ y_{i}^{1} \\ y_{n}^{2} \end{matrix}] = [\begin{matrix} f_{n}^{1} \\ f_{n}^{1} + T_{i}^{T} + f_{i} 2 \\ f_{n}^{2} \end{matrix}]

앞서 설명한 Nitsche 방법론 기반의 멀티패치 아이소-지오 메트릭 해석과 마스터-슬레이브 방법론 기반의 해석법을 경계 면의 조건, 강성 행렬의 대칭성 및 구현 난이도 측면에서 간단 하게 비교하였다.Table 1

Table 1

Comparison of multipatch isogeometric analysis

Nitsche based method	Master-slave method
Non-conforming interface	Conforming interface
Symmetric stiffness matrix	Non-symmetric stiffness matri
Complicated implementation	Simple implementation

2.5. 수치예제 결과 및 고찰

Nitsche 방법론에 기반한 멀티패치 문제의 아이소-지오메 트릭 해석의 정확성을 검증하기 위해서 정해(exact solution) 를 가지는 2차원 멀티패치 평면 빔 문제에 대해서 아이소-지오 메트릭 해석을 수행하고자 한다. Fig. 2에 나타낸 것처럼 원점 에서 고정되고 하중 P를 받는 빔에 대해서, 길이 방향 가운데 면을 분리하여 2개의 패치로 구성하였으며 패치간 곡면 조정 점이 일치하지 않는 비순응(non-conforming) 문제로 패치를 구성하였다. 수치예제에서 3차의 NURBS 기저함수를 사용 하였고 탄성계수 E = 3×10⁷, 푸아송 비 ν = 0.3, 길이 L = 48, 높이 D = 6, 두께 t = 1로 가정하였다. 아이소-지오메트릭 해 석의 표준 격자 세분화(refinement) 방법인 h-세분화를 적용 하여 Coarse 격자, Fine 격자 두 가지 경우를 해석하였고, Coase 격자의 경우 패치 1의 자유도(degree of freedom)은 266, 패치 2의 자유도는 50 이며, Fine 격자의 경우 패치 1의 자유도는 585, 패치 2의 자유도는 50 이다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2017-030-06/07COSEIK-30-515/images/COSEIK-30-515_F2.jpg

Fig. 2

2-dimensional multipatch beam

두 경우의 해석에서 얻어진 $σ_{x x}, σ_{x y}$ 응력을 Fig. 3과 Fig. 4에 나타내었다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2017-030-06/07COSEIK-30-515/images/COSEIK-30-515_F3.jpg

Fig. 3

Stress contour of beam(coarse mesh)

https://cdn.apub.kr/journalsite/sites/jcoseik/2017-030-06/07COSEIK-30-515/images/COSEIK-30-515_F4.jpg

Fig. 4

Stress contour of beam(fine mesh)

Fig. 2의 수치예제에서 응력에 대한 정해는 다음 식 (26)와 같이 주어진다. y = 1.5 지점에서 응력을 계산하여 아래의 그림 Fig. 5, Fig. 6에서 수치해석 결과와 정해를 비교하였다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2017-030-06/07COSEIK-30-515/images/COSEIK-30-515_F5.jpg

Fig. 5

Stress comparison between exact solution and isogeometric analysis(coarse mesh)

https://cdn.apub.kr/journalsite/sites/jcoseik/2017-030-06/07COSEIK-30-515/images/COSEIK-30-515_F6.jpg

Fig. 6

Stress comparison between exact solution and isogeometric analysis(fine mesh)

(26)

σ_{x x} (x, y) = \frac{P (L - x) y}{I}, σ_{x y} (x, y) = - \frac{P}{2 I} (\frac{D^{2}}{4} - y^{2})

경계조건의 영향을 받게 되는 x = 0 근처를 제외하고는 수직 응력과전단응력이정해와잘일치하는경향을확인할수있으며, 특히 패치간 경계면의 영향을 확인하기 위해서 패치의 경계인 x = 24에서 전단응력을 계산하여 Fig. 7에 비교하였고 정해와 일치하는 것을 확인 가능하다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2017-030-06/07COSEIK-30-515/images/COSEIK-30-515_F7.jpg

Fig. 7

Stress comparison between exact solution and isogeometric analysis at the interface(fine mesh)

멀티패치 아이소-지오메트릭 Nitsche 방법론을 검증하기 위한 두 번째 수치예제로 Fig. 8의 인장하중 T_x를 받는 평판의 구멍에서의 응력집중 문제를 다루기로 한다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2017-030-06/07COSEIK-30-515/images/COSEIK-30-515_F8.jpg

Fig. 8

Hole stress concentration in plate

위의 수치예제에서 응력에 대한 정해는 다음 식 (27)과 같 이 주어진다.

(27a)

σ_{r r} (r, θ) = \frac{T_{x}}{2} (1 - \frac{R^{2}}{r^{2}}) + \frac{T_{x}}{2} (1 - 4 \frac{R^{2}}{r^{2}} + 3 \frac{R^{4}}{r^{4}}) cos 2 θ

(27b)

σ_{θ θ} (r, θ) = \frac{T_{x}}{2} (1 + \frac{R^{2}}{r^{2}}) - \frac{T_{x}}{2} (1 + 3 \frac{R^{4}}{r^{4}}) cos 2 θ

문제에서 전체 크기 8×8, 구멍의 반지름 1, 두께 1인 문제 에 대해서 탄성계수 E = 3×10⁷, 푸아송 비 ν = 0.3로 가정하 였다. 대칭 경계조건을 활용하여 1/4 문제에 해석을 수행하고 2차의 NURBS 기저함수를 사용하였으며, 패치 1에서 원주방향 자유도는 20, 패치 2에서 원주방향 자유도는 36 이다. 경계면의 영향을 확인하기 위해서 반지름 2, 1.5, 1.25인 면을 패치 간의 경계면으로 정의하였고 각 경우에 대해서 Nitsche 방법론과 마스터-슬레이브 방법론을 사용하여 아이소-지오메트릭 해석을 수행한 후 수직응력σxx의 값을 Fig. 9, Fig. 10, Fig. 11에 나타내었다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2017-030-06/07COSEIK-30-515/images/COSEIK-30-515_F9.jpg

Fig. 9

Multipatch with the interface at R=2 and stress contour comparison

https://cdn.apub.kr/journalsite/sites/jcoseik/2017-030-06/07COSEIK-30-515/images/COSEIK-30-515_F10.jpg

Fig. 10

Multipatch with the interface at R=1.5 and stress contour comparison

https://cdn.apub.kr/journalsite/sites/jcoseik/2017-030-06/07COSEIK-30-515/images/COSEIK-30-515_F11.jpg

Fig. 11

Multipatch with the interface at R=1.25 and stress contour comparison

두 가지 방법론 모두 정 해를 잘 묘사하는 것을 확인할 수 있으며, 멀티패치에서 두 가지 방법론의 차이를 명확하게 파학 하기 위해서 경계면에서 다음 식과 같이 패치 간 응력차이를 정의하였다.

(28)

Δ (r, θ) = \frac{| σ_{r r}^{2} (r, θ) - σ_{r r}^{1} (r, θ) |}{σ_{r r}^{e x a c t} (r, θ)}

이브 방법, Nitsche 방법 그리고 Nitsche 방법에 2.3절에서 설명한 NURBS 곡면 기반의 응력 복원법을 적용하여 응력값을 계산하고 식 (24)에 정의된 응력차이를 통해 상대적인 오차를 비교하였다.

Fig. 12, Fig. 13, Fig. 14에서 경계면을 따라 즉 아래면을 θ = 0°으로 정의하고 시계방향을 따라 θ = 90°까지 증가하면서 응력의 상대적 오차의 크기를 나타내었다. θ = 0° 근처에서는 정해로 주어지는 분모의 응력값이 작아서 오차가 큰 것으로 보이나 전체적으로는 경계면이 R=2인 경우 2%, R=1.5인 경우 5%, R=1.25인 경우 13% 이내의 패치간 응력 오차를 가지는것을확인할수있다. 경계면의반지름이작아질때응력 오차가 커지는 이유는 구멍 주위의 응력의 변화가 큰 부분이 Coarse 격자로 포함되는 이유인 것으로 판단된다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2017-030-06/07COSEIK-30-515/images/COSEIK-30-515_F12.jpg

Fig. 12

Relative stress difference at the interface(R=2)

https://cdn.apub.kr/journalsite/sites/jcoseik/2017-030-06/07COSEIK-30-515/images/COSEIK-30-515_F13.jpg

Fig. 13

Relative stress difference at the interface(R=1.5)

https://cdn.apub.kr/journalsite/sites/jcoseik/2017-030-06/07COSEIK-30-515/images/COSEIK-30-515_F14.jpg

Fig. 14

Relative stress difference at the interface(R=1.25)

Nitsche 방법론의 결과와 마스터-슬레이브 방법론의 결과를 비교하면경계면을 따라 거의 동일한 응력 오차값을가지는것을 확인할 수 있으며, NURBS 곡면 기반의 응력 복원법을 적용 하여 응력을 계산하였을 경우, 빨간 점선으로 표기한 결과에서 패치간 경계면에서 응력차이가 현저하게 줄어드는 것을 확인 할 수 있다. 마지막으로 주어진 문제의 3차원 응력분포를 Fig. 15에 도시하여 세가지 방법론을 사용한 멀티패치 해석 결과를 비교하였고, NURBS 곡면 기반의 응력 복원법을 적용할 경우 패치 내부 및 경계면에서 연속적인 응력 모사가 가능한 것을 알 수 있다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2017-030-06/07COSEIK-30-515/images/COSEIK-30-515_F15.jpg

Fig. 15

Stress contour of multipatch analysis

3. 결 론

본 논문에서는 멀티패치를 가지는 아이소-지오메트릭 해석 문제에대해서경계면의응력차를 최소화 하는방법으로Nitsche 방법론과 조정점의 마스터-슬레이브 관계를 이용한 방법론을 구현하고 그 결과를 비교하는 연구를 수행하였다. 마스터- 슬레이브 관계를 이용한 방법론은 상대적으로 간단하게 구현 가능하나 순응 경계면에서만 적용가능한 단점을 가지나, Nitsche 방법론을 도입할 경우 비순응 경계면에서 적용 가능하고 대칭의 강성행렬을 얻게 되는 장점을 가진다.

수치 예제를 통해서 정해를 가지는 문제에 대한 Nitsche 방법론 기반의 멀티패치 아이소-지오메트릭 해석을 수행하여 경계면에서 적확한 결과가 나오는 것을 확인하였으며, 응력 집중을 가지는 문제에 대해서 Nitsche 방법론의 결과와 마스터 -슬레이브 관계를 이용한 방법이 유사한 결과를 가지는 것을 알 수 있었다. 특히 후처리 과정에서 NURBS 곡면 기반의 응력 복원법을 적용할 경우 패치 내부 및 경계면에서 연속적인 응력 모사가 가능한 것을 알 수 있다.

본 연구의결과는여러물체가 접합된비순응경계면을가지는 멀티패치아이소-지오메트릭해석에직접적용 가능하며, 해석의 후처리과정으로 본 연구에서검토한NURBS 곡면 기반의 응력 복원법을 사용할 수 있다.

감사의 글

본 연구는 2017년도 정부(과학기술정보통신부)의 재원으로 한국기계연구원 개인시드창의형과제(NK207D)와 한국에너지 기술평가원(KETEP)의 ‘네트워크 기반 유체기기 고효율화 고급트랙’ 인력양성 사업(No. 20174010201350)의 지원을 받아 수행되었습니다.

References

Brivadis, E., Buffa, A., Wohlmuth, B. and Wunderlich, L., 2015. Isogeometric Mortar Methods.. Comput. Methods Appl. Mech. Eng., 284, pp.292-319.

10.1016/j.cma.2014.09.012

Cottrell, J., Hughes, T. and Reali, A., 2007. Studies of Refinement and Continuity in Isogeometric Structural Analysis.. Comput. Methods Appl. Mech. Eng., 196, pp.4160-4183.

10.1016/j.cma.2007.04.007

Dolbow, J. and Harari, I., 2009. An Efficient Finite Element Method for Embedded Interface Problems.. Int. J. Numer. Methods Eng., 78, pp.229-252.

10.1002/nme.2486

Ha, Y.D., Yoon, M.H. and Cho, S., 2012. Isogeometric Shape Sensitivity Analysis in Generalized Curvilinear Coordinate Systemsn Isogeometrical approach to Error Estimation and Stress Recovery.. J. Comput. Struct. Eng., 25, pp.497-504.

10.7734/COSEIK.2012.25.6.497

Hassani, B., Ganjali, A. and Tavakkoli, M., 2012. An Isogeometrical approach to Error Estimation and Stress Recovery.. Eur. J. Mech. A, Solids, 31, pp.101-109.

10.1016/j.euromechsol.2011.08.001

Hughes, T., Cottrell, J. and Bazilevs, Y., 2005. Isogeometric Analysis: CAD, Finite Elements, Exact Geometry and Mesh Refinement.. Comput. Methods Appl. Mech. Eng., 194, pp.4135-4195.

10.1016/j.cma.2004.10.008

Nguyen, V., Kerfridenu, P., Brino, M., Bordas, S. and Bonisoli, E., 2014. Nitsche ?(tm)s Method for Two and Three Dimensional NURBS Patch Coupling.. Comput. Mech., 53, pp.1163-1182.

10.1007/s00466-013-0955-3

Sanders, J.D., Laursen, T. and Puso, M.A., 2011. A Nitsche Embedded Mesh Mothod.. Comput. Mech., 49, pp.243-257.

10.1007/s00466-011-0641-2

Journal of the Computational Structural Engineering Institute of Korea ISSN:1229-3059(Print) 2287-2302(Online) 한국전산구조공학회 논문집

Preview

Isogeometric Analysis for Two-dimensional Multipatch Model

ABSTRACT

MAIN

(1)

(2)

(3)

(4)

(5)

(6)

(7)

(8)

Fig. 1

(9a)

(9b)

(9c)

(9d)

(9e)

(10a)

(10b)

(11)

(12)

(13)

(14)

(15)

(16a)

(16b)

(17a)

(17b)

(17c)

(17d)

(18)

(19)

(20)

(21)

(22)

(23)

(24)

(25)

Table 1

Fig. 2

Fig. 3

Fig. 4

Fig. 5

Fig. 6

(26)

Fig. 7

Fig. 8

(27a)

(27b)

Fig. 9

Fig. 10

Fig. 11

(28)

Fig. 12

Fig. 13

Fig. 14

Fig. 15

감사의 글

References