곡면의 특이적분을 위한 가상 매핑 방법

Pseudo Mapping Method for Singular Integral of Curved Panels

Article information

J. Ocean Eng. Technol. 2019;33(1):17-25
Publication date (electronic) : 2019 February 27
doi : https://doi.org/10.26748/KSOE.2019.001
*Department of Naval Architecture and Ocean Engineering, PNU, Busan, Korea
이익재,*orcid_icon, 권순홍*orcid_icon
*부산대학교 조선해양공학과
Corresponding author Ik-Jae Lee: +82-51-510-2753, ick7140@pusan.ac.kr
It is noted that this paper is revised edition based on proceedings of the 7th Green’s Function Seminar in Shanghai, China.
Received 2019 January 9; Revised 2019 January 29; Accepted 2019 February 22.

Trans Abstract

A numerical method is suggested for evaluating the singular integral of curved panels in the higher-order boundary element method. Two-step mapping procedures that are significantly related to the physical properties of singular behaviors were developed and illustrated. As a result, the singular behaviors were significantly alleviated, and the efficiency and robustness of the present method for tangentially and axially deformed elements were proven. However, inaccuracies and numerical instabilities of twisted elements were discovered as a result of nonlinearities.

1. 서 론

과학 및 공학 분야에서 복잡한 기하학적 경계를 포함하는 경계치 문제는 경계적분방정식을 통해 모델링할 수 있으며, 경계요소법은 이 경계적분방정식을 효율적으로 풀 수 있는 중요한 수치 계산 방법 중 하나로 자리매김해 왔다. 경계요소법은 해양유체, 전자기, 음향, 탄성재료 등의 다양한 영역의 해석을 위해 적용될 수 있으며, 특히 관성이 지배적인 해양 환경에 놓인 물체의 운동 해석에서는 포텐셜 이론에 기반한 유체력 계수 등을 계산할 때 주로 사용되고 있다.

경계요소법은 라플라스(Laplace) 또는 헬름홀츠(Helmholtz) 방정식을 만족하는 연속체 영역에서 복잡한 기하 형상에서 적용되는 경계치 문제를 표현한 경계적분방정식을 이산화시켜 푸는 가장 대표적인 방법으로, 풀고자 하는 전 영역이 아닌 이 영역들을 둘러싸고 있는 경계들에 적용되는 경계조건만으로 문제를 풀 수 있다는 장점을 가지고 있다. 하지만 경계적분방정식은 무한 원방 경계를 가지는 라플라스 방정식의 기본해(Fundamental solution)인 랜킨 소스(Rankine source) 항과 미분항인 다이폴(Dipole)을 포함하고 있으며, 이들은 공간에 대한 델타 함수(Delta function)의 충격 응답 함수(Impulse response function)로서 충격이 가해진 위치에서 특이점(Singular point)을 가지는 문제점을 안고 있다.

한편, 경계요소법은 선형 요소 기반의 일정요소법(Constant panel method)과 고차 다항식 근사를 이용한 고차경계요소법(Higher-order boundary element method)으로 나눠질 수 있는데 일정요소법과 고차경계요소법의 수치계산 정확도 및 효율은 Liu et al.(1991) 등에 의해서 비교된 바 있으며, 고차경계요소법의 넓은 적용범위와 높은 정확도는 잘 알려져 있다. 경계적분방정식에서는 앞서 언급한 특이점 위치에서 특이적분을 수행해야 한다는 어려움이 있는데, 고차경계요소법에서는 복잡한 곡면의 경계 요소들이 고차 다항식으로 근사되어 높은 정확도의 수치 적분을 위해서는 평면 요소에 비해 상대적으로 많은 계산 시간이 요구된다.

곡면에서의 특이점 적분법은 오래 전부터 많은 연구자들에 의해서 개발되어 왔는데, 이 중 Guiggiani and Gigante(1990)Guiggiani et al.(1992)는 경계적분방정식이 가질 수 있는 특이성(Singularity)을 약한 특이(Weakly singular), 강한 특이(Strongly singular), 초특이(Hypersingular)의 세 가지 경우로 분리하여 테일러 급수 전개(Taylor’s series expansion)를 이용한 준해석적(Semi-analytical) 방법을 적용한 바 있다. 최근 Rong et al.(2014)Guiggiani and Gigante(1990)Guiggiani et al.(1992)의 테일러 급수 전개 접근 방식에서 특이점 부근의 각-방향(Angular-direction) 특이 거동을 완화시킬 수 있는 특정조건들을 발견하여 Conformal polar coordinate transformation(CPCT)라고 이름한 매핑(Mapping) 방법을 고안하였고, 삼각 요소에 대한 CPCT법의 효능을 검증한 바 있다. 그리고 Lv et al.(2015)Rong et al.(2014)의 삼각 요소 문제를 사각 요소에서 적용할 수 있는 CPCT법으로 확장하여 정확도 및 효율에 대해 검증하였다.

일반적으로 고차경계요소법에서는 3차원 영역에서의 곡면 요소를 2차원 국부 요소(Local element)로 매핑한 후 가우스-르장드르 구적법(Gauss-Legendre quadrature) 등의 수치 적분방법을 사용하고 있는데, 특이점을 포함하는 요소의 경우는 특이점을 기준으로 극좌표계를 만들고, 이 특이점에서 국부 요소의 꼭지점들까지 선을 그어 여러 개의 삼각형 요소들로 나눠서 각 삼각형에 대해 수치적분을 수행하는 방법이 주로 사용되고 있다.

앞서 언급한 CPCT법은 3차원 공간에서 국부 영역(Local domain) 또는 주영역(Master domain)으로의 일차적인 매핑 이후에 적용되며, 특정 조건들을 만족하는 또 다른 국부 영역으로의 추가적인 매핑을 의미하는데, 특정 조건들은 랜킨 소스를 국부 영역 특이점 기준으로 테일러 급수 전개할 때 나타나는 여러 항들 중 가장 지배적인 일차항의 각의존성(Angular dependency)을 제거하여 상수로 만드는 조건을 의미한다. 즉, CPCT법을 적용한 매핑 후의 새롭게 정의된 국부 영역에서 극좌표계로 변환하여 적분을 수행하게 된다.

앞서 언급했듯이, Rong et al.(2014)Lv et al.(2015)에 의해 CPCT법의 정확도 및 효율은 검증이 되었지만, 두 저자가 제시한 매핑 방법들에서는 각방향 특이성 완화 메커니즘에 대한 물리적 의미들의 해석이 어려우며, 새롭게 매핑되는 영역으로부터 역변환(Inverse transformation)을 취해야하는 다소 부자연스럽고 복잡한 절차를 요구하고 있다. 반면 본 연구에서는 CPCT법의 물리적인 의미를 명확하게 보여주면서도 수식이 간단하고 Rong et al.(2014)Lv et al.(2015)가 제시한 방법들의 정확성과 효율성을 동시에 가지는 가상 매핑 기법을 개발하였다. 본 연구에서 제시한 방법은 2장에서 자세히 언급하였고, 3장에서는 이산화된 경계적분식에서 특이점을 가지는 국부 요소의 극좌표 변환 및 테일러 급수 전개를 통해 CPCT에서 제시한 두 가지 조건을 유도하였다. 4장에서는 위 두 조건에 해당되는 가상 매핑 방법 두 단계의 절차를 수식을 통해 정식화 하였으며, 그림을 통해 시각적으로 도식화하였다. 5장에서는 가상 매핑 방법을 적용하기 전, 후의 특이 요소의 수치 적분 결과를 나타내었으며, 6장에서 수치 결과와 가상 매핑 방법에 대한 고찰 및 결론으로 마무리하였다.

2. 가상 매핑 방법

가상 매핑은 CPCT법에서 주어진 두 가지 조건과 특이점 부근의 물리량 왜곡(Distortion) 현상을 기반으로 하는데, 두 조건은 3차원 곡면이 국부 영역으로 옮겨가는 과정에서 발생할 수 있는 다양한 왜곡 현상 중 두 가지 선형 왜곡 현상을 의미한다. 첫 번째 조건은 마주하는 변(삼각 요소의 경우 한 꼭지점과 마주하는 변)들의 접선 방향 왜곡에 해당하고, 다른 한 조건은 수직 방향 압축(또는 인장)에 해당한다. 이에 가상 매핑 방법 또한 두 단계로 나누어 각 단계가 두 조건에 대응될 수 있도록 했으며, 각 매핑을 접선 변형(Tangential deformation)과 수직 변형(Normal deformation)으로 이름하였다. 가상 매핑을 이용하면 일차적인 매핑에서 발생한 접선 및 수직 방향의 선형 변형 현상을 2차원 공간 내에서 반대편으로 되돌려 줌으로서 특이점 주위를 둘러싸는 특이 거동을 상당히 완화시킬 수 있다. 이 두 변형들은 곡면뿐 아니라 평면에서도 발생할 수 있지만, 평면 요소의 경우는 수직벡터를 기준으로 새로운 국부좌표계를 형성하면 3차원 공간의 요소를 아무 변형없이 2차원으로 옮겨올 수 있으므로 가상 매핑을 사용할 필요가 없게 된다. 이에 본 연구에서는 곡면 요소에 대해서만 고려하였다. 실제 평면 요소에서 선형 보간함수를 이용해 매핑을 수행한 경우는 특이점 주위 왜곡 현상들의 대부분을 제거할 수 있다.

한편, Rong et al.(2014)Lv et al.(2015)은 새롭게 정의된 매핑 영역에서 반경-방향(Radial-direction) 적분에 대해서는 해석적인 방법을 사용하였는데, 저자의 경험에 따르면 배가 전진하며 운동하는 특수한 문제의 경우는 경계적분방정식에 m-항과 같은 복잡한 경계조건을 필요로 하며, Rong et al.(2014)Lv et al.(2015)에서처럼 해석적 적분을 수행하기 위해서는 수치 미분을 포함하는 m-항을 테일러 전개하는 꽤나 복잡한 계산을 필요로 하게 된다. 이에 본 연구에서는 Rong et al.(2014)Lv et al.(2015)의 방법들과 달리 CPCT법을 적용한 뒤 가우스 구적법을 이용한 수치 적분만을 수행하였으며, CPCT법의 적용이 더 복잡한 사각 요소에 대해서만 수행하였다. 실제로 CPCT법을 통한 특이점 적분의 정확도와 효율은 위 두 저자의 연구에 의해 검증되었고, 동일한 두 조건을 만족한다면 수치 정확도의 차이는 발생하지 않는다. 이에 본 논문에서는 CPCT법에서 제시된 두 조건들이 가지는 좀 더 넓은 물리적 의미들을 근거하여 각 조건들에 대응하는 두 가지 매핑 방법의 수식 유도과정과 수치 결과에서 나타나는 비선형적 왜곡현상들의 문제점을 언급하는데 초점을 두었다. 이후부터는 본 연구에서 제시한 방법을 가상 매핑 방법(Pseudo mapping method, PMM)으로 명명한다.

3. 특이 적분

3.1 경계적분방정식

라플라스 방정식을 지배방정식으로 하는 속도 포텐셜 ϕ와 그린 함수 G를 통해 그린 제2 항등식으로부터 다음과 같은 경계적분방정식을 유도할 수 있다.

(1) cϕ(x)+sϕ(x˜)G(x,x˜)n(x˜)dS(x˜)=sϕ(x˜)n(x˜)G(x,x˜)dS(x˜)

여기서, c=0,xS, Ωμ,xS4π,xΩ 이며, xSsmooth의 매끄러운 면에서는 입체각(Solid angle) μ는 2π를 가진다.

또한, S, Ω는 연속체를 둘러싸는 경계와 경계의 내부 영역을 의미하며, 해양 유체 영역에서의 S는 자유 수면, 바닥 경계, 원방 경계 등을 포함하는 모든 유체 경계를 의미한다. X~, X는 실제 유체 영역 필드점(Field point)의 위치 벡터와 소스점(Source point)에서의 위치 벡터를 의미하며 다음과 같이 쓸 수 있다.

(2) x˜=(x˜1,x˜2,x˜3)
(3) x=(x1,x2,x3)

그리고 그린 함수 G는 경계조건이 적용되지 않은 소스 1/r이며, r=x~-x2는 필드점과 소스점 사이의 거리를 의미한다. 소스점 x가 유체의 경계에 위치할 때, 식 (1)은 이산화 과정을 통해 다음과 같이 특이점이 포함된 요소를 분리해 낼 수 있다.

(4) α(x˜)+β(x˜,x)+β*(x˜,x)=γ(x˜,x)+γ*(x˜,x),x˜S
(5) α(x˜)=μϕ(x˜)β(x˜,x)=e(ee*)Seϕ(x)G(x˜,x)n(x)dS(x),β*(x˜,x)=P.V.Se*ϕ(x)G(x˜,x)n(x)dS(x)γ(x˜,x)=e(ee*)Seϕ(x)n(x)G(x˜,x)dS(x),γ*(x˜,x)=Se*ϕ(x)n(x)G(x˜,x)dS(x)

여기서, β*, γ*는 각각 강한 특이성과 약한 특이성을 가지는 특이점 적분 항들이며, P.V.는 코시 주요값을 의미한다.

3.2. 특이점을 포함하는 적분

Fig. 1은 3차원 사각 곡면 형상 및 물리량을 9개 점으로 근사해 국부 요소로의 매핑 과정과 동시에 특이점 주위의 특이성 분포 왜곡 현상을 나타낸다.

Fig. 1

Mapping from 3-D curved panel to local element and distortion of singular behavior near the singularity

각 절점에 위한 위치 벡터와 물리량들은 라그랑지 보간함수(Lagrange interpolating function)을 이용해 다음과 같이 근사된다.

(6) xi=k=19xikNk(ξ1,ξ2),(i=1,2,3)
(7) ϕ=k=19ϕkNk(ξ1,ξ2)
(8) ϕn=k=19(ϕn)kNk(ξ1,ξ2)

여기서 라그랑지 보간함수 Nk(ξ1, ξ2)는 다음과 같다.

(9) Nk=14ξ1kξ1(1+ξ1kξ1)ξ2kξ2(1+ξ2kξ2)k=1,3,5,7Nk=12(1ξ12)ξ2kξ2(1+ξ2kξ2)k=2,6Nk=ξ1kξ1(1+ξ1kξ1)12(1ξ22)k=4,8Nk=(1ξ12)(1ξ22)k=9

여기서, kFig. 1에 부여된 절점 번호를 의미한다.

식 (4)의 특이점을 포함하는 항들을 식 (6)-(9)를 이용하여 국부 요소로 매핑하고 이산화하면 다음 결과를 얻을 수 있다.

(10) γ*=k=19(ϕn)k1111Nk(ξ)r1(ξ˜,ξ)J(ξ)2dξ1ξ2
(11) β*=k=19ϕk1111Nk(ξ)x(r1(ξ˜,ξ))n(ξ)J(ξ)2dξ1dξ2

여기서, J2는 자코비안 놈(Jacobian norm), n은 유체 경계 바깥 방향으로의 법선 벡터이며 다음과 같이 얻을 수 있다.

(12) J(ξ)2=x(ξ)ξ1×x(ξ)ξ22
(13) n(ξ)=J(ξ)J(ξ)2

식 (6)-(8)과 같이 고차경계요소법에서는 3차원 영역에서 국부영역으로의 매핑 과정에서 위치 벡터, 속도 포텐셜 등이 고차다항식으로 근사하게 되는데, 본 연구에서는 수치 결과 비교를 위해 특이점 주위의 기하학적 변형에 따른 변화만을 고려하여 속도 포텐셜과 그 미분 값은 상수로 가정하였다. 실제로 라그랑지 보간함수는 경우 특이성(Singularity)을 전혀 가지지 않는 매끄러운 함수이며, 본 연구에서 살펴볼 특이점 적분에는 거의 영향을 끼치지 못한다. 따라서, 식 (7)의 노드 포텐셜과 법선방향 미분 값들을 모두 1인 간단한 문제로 가정하면 라그랑지 보간함수들의 완전성 조건(Completeness condition)에 의해 위 식들은 다음과 같이 된다.

(14) γ*=1111r1(ξ˜,ξ)J(ξ)2dξ1dξ2
(15) β*=1111x(r1(ξ˜,ξ))n(ξ)J(ξ)2dξ1dξ2

위 두 가지 경우처럼 일반적으로 사용하는 경계적분방정식에서 얻을 수 있는 특이 적분은 약한 특이, 강한 특이의 두 종류이지만 음향 문제에 사용되는 Burton-Miller 공식 등에서는 다음과 같은 초특이 적분을 포함하는 항들도 다루어지고 있다.

(16) Ihyper=1111r3(ξ˜,ξ)J(ξ)2dξ1dξ2

4. 가상 매핑 방법

4.1. 거리 r의 테일러 급수 전개

먼저, 소스점과 필드점의 i 방향에 대한 위치 벡터 성분 차는 테일러 급수 전개로부터 근사적으로 다음과 같이 쓸 수 있다.

(17) xix˜i=ρ(u1icosθ+u2isinθ)+ρ22(u11icos2θ+u22isin2θ+u12icosθsinθ)+O(ρ3)=ρAi(θ)+ρ2Bi(θ)+O(ρ3)

여기서 사용된 극좌표계는 (ξ1, ξ2)=(ξ~1+ρ cosθ, ξ~2+ρ sinθ)의 관계를 가진다. 그리고 uαi, uαβixi/ξαξ=ξ~2xi/ξαξβξ=ξ~로 각각 계산할 수 있으며, 전자는 공변 기저벡터(Covariant basis vector)이고, 후자는 공변 벡터 방향의 곡률(Curvature)로 생각할 수 있다(α, β = 1, 2).

위 결과로부터 거리의 제곱 r2은 다음과 같이 구할 수 있다.

(18) r2=i=13(xix˜i)2=ρ2(|u1|2cos2θ+|u2|2sin2θ+2u1u2cosθsinθ)+ρ3{u1u11cos3θ+u2u22cos3θ+(u1u22+2u1u12)cos2θsinθ+(u2u11+2u2u12)cosθsin2θ}+O(ρ4)=ρ2AiAi+ρ3AiBi+O(ρ4)

여기서, uα=x/ξαξ=ξ~, uαβ=2x/ξαξβξ=ξ~이고, 반복된 인덱스(Index)는 아인슈타인의 합 규약(Einstein summation convention)을 따른다(AiAi=A12+A22+A32).

4.2. 각-방향 의존성을 완화시키는 두 가지 조건

식 (18)에서 AiBi는 각의 변화에 따라 바뀌는 값으로 각-의존성을 가지고 있다. 앞서 언급한 CPCT법의 두 가지 조건은 각-의존적인 첫 번째 항의 AiAi를 상수로 만들어주는 역할을 하며 다음과 같이 쓸 수 있다.

(19) u1u2=0|u1|=|u2|

위 두 조건들이 각각 가지는 물리적 의미들은 다음과 같다.

(1) 곡면(또는 뒤틀린 평면) 특이점의 공변 기저벡터들은 서로 직교(Orthogonal)한다.

(2) 두 방향의 공변 기저벡터들의 크기는 같다.

하지만 일반적인 ξ-영역에서는 위 조건들을 만족시키지 못한다(즉, u1·u20, u1u2). 따라서, 위 두 조건을 만족시키는 ζ-영역을 찾는 것이 본 연구의 목표이다. 먼저, 첫 번째 조건을 만족시키는 η-영역과 두 조건을 모두 만족시키는 ζ-영역에서의 공변 기저벡터들은 연쇄법칙(Chain rule)에 의해 다음과 같이 정의할 수 있다.

(20) vα=xξ1|ξ=ξ˜ξ1ηα|η=η˜+xξ2|ξ=ξ˜ξ2ηα|η=η˜=uβξβηα|η=η˜
(21) wα=xξ1|ξ=ξ˜(ξ1η1|η=η˜η1ζα|ζ=ζ˜+ξ1η2|η=η˜η2ζα|ζ=ζ˜)+xξ2|ξ=ξ˜(ξ2η1|η=η˜η1ζα|ζ=ζ˜+ξ2η2|η=η˜η2ζα|ζ=ζ˜)=uγξγηβ|η=η˜ηβζα|ζ=ζ˜

매핑을 통해 최종적으로 ζ-영역에서는 다음과 같은 두 조건을 만족시키도록 한다.

(22) w1w2=0|w1|=|w2|

여기서, ξ=ξ(η), η=η(ζ)이며, 다음과 같은 관계를 갖는다.

(23) ξα=kξαkNk(η1,η2),(α=1,2)ηα=kηαkNk(ζ1,ζ2)

또한, 여기서 사용된 보간 함수는 식 (6)처럼 반드시 고차 다항식일 필요가 없고, 식 (22) 조건은 위치 벡터들의 1차 미분만을 사용하므로 선형 보간함수를 사용하는 것으로 충분히 만족시킬 수 있다.

4.3. 가상 매핑 방법

가상 매핑 방법은 위 두 조건에 따라 두 단계로 분리할 수 있다. 첫 번째 매핑된 영역에서는 식 (20)η-영역 공변 기저벡터들이 서로 직교하도록 만들어주며, Fig. 2와 같이 특이점 부근에서 접선 방향의 왜곡 현상을 완화시켜주는 역할을 한다. 두 번째 매핑된 영역에서는 식 (21)ζ-영역 공변 기저벡터들이 첫 번째 조건을 만족시키면서 동시에 서로 크기가 같도록 만들어준다. Fig. 3에서 볼 수 있듯이, 특이점에서의 서로 직교한 방향의 서로 다른 크기들이 같도록 해서 최종적으로 각방향 의존적인 문제를 완화시키는 역할을 한다.

Fig. 2

Pseudo mapping method: Tangential deformation

Fig. 3

Pseudo mapping method: Normal deformation

먼저 Fig. 2의 가상의 사각형은 ξ-영역에서 η-영역으로 매핑되며 식 (23)과 같이 라그랑지 보간함수로 쓸 수 있다. 이 때, ξ-영역에서의 노드 값들은 다음과 같다.

(24) (ξ11,ξ21)=(1+a,1)(ξ12,ξ22)=(1+a,1)(ξ13,ξ23)=(1a,1)(ξ14,ξ24)=(1a,1)

이로부터 η-영역에서 정의된 공변 기저벡터는 연쇄법칙에 의해 다음과 같이 얻을 수 있다.

(25) v1=u1ξ1,η1|η=η˜+u2ξ2,η1|η=η˜=u1v2=u1ξ1,η2|η=η˜+u2ξ2,η2|η=η˜=au1+u2

첫 번째 조건, v1 ․ v2 = 0을 통해 다음과 같이 얻을 수 있다.

(26) a=u1u2|u1|2

마찬가지 Fig. 3에서 η-영역의 가상의 사각형 노드 값들은 다음과 같다.

(27) (η11,η21)=(1,b)(η12,η22)=(1,b)(η13,η23)=(1,b)(η14,η24)=(1,b)

이로부터 ζ-영역에서 정의된 공변 기저벡터는 연쇄법칙에 의해 다음과 같이 얻을 수 있다.

(28) w1=v1η1,ζ1|ζ=ζ˜+v2η2,ζ1|ζ=ζ˜=v1w2=v1η1,ζ2|ζ=ζ˜+v2η2,ζ2|ζ=ζ˜=bv2

그리고 두 조건은 다음과 같이 만족시킬 수 있다.

(29) w1w2=0|w1|=|w2|bv1v2=0|v1|=b|v2|

이에 b는 다음과 같이 구할 수 있다.

(30) b=|v1||v2|

따라서, ξ-영역과 ζ-영역사이의 관계를 식 (23)을 통해 얻을 수 있고, 이로부터 ζ-영역에서 옮겨진 특이점 위치와 꼭지점 노드 값들은 다음과 같다.

(31) (ζ˜1,ζ˜1)=(ξ˜1+aξ˜2,ξ˜2/b)
(32) (ζ11,ζ21)=(1a,1/b)(ζ12,ζ22)=(1a,1/b)(ζ13,ζ23)=(1+a,1/b)(ζ14,ζ24)=(1+a,1/b)

식 (24)(27)은 가상의 사각형 노드 값들인 반면, 식 (32)는 실제 계산되는 영역의 노드 값들을 가리킨다. 결과적으로 ζ-영역으로 매핑된 적분방정식은 다음과 같이 쓸 수 있다.

(33) γ*=G(ζ˜,ζ)J(ζ)2bdζ1dζ2β*=Gn(ζ˜,ζ)J(ζ)2bdζ1dζ2

여기서 ξ-영역에서 ζ-영역으로의 매핑 자코비안은 b이며, 면 문제에서 자코비안은 매핑 전, 후 영역의 면적비로 정의될 수 있고 접선 방향 변화는 면적에 영향을 주지 않는다. 이에 수직방향 압축(또는 인장)에 의한 면적 변화비인 b가 자코비안 값이 됨을 알 수 있다.

최종적으로 ζ-영역의 특이점 적분을 수행하기 위해 특이점 중심의 극좌표계를 도입하면 4개의 삼각형 요소로 나눌 수 있다. 이 때, 적분방정식은 다음과 같이 쓸 수 있다.

(34) γ*=Δ=14ΔG(ρ,θ)J(ρ,θ)2bρdρdθβ*=Δ=14ΔGn"(ρ,θ)J(ρ,θ)2bρdρdθ

여기서, G''(ρ,θ)=G(ξ(ζ(ρ,θ))), J''(ρ,θ)=J(ξ(ζ(ρ,θ)))이며 J''=b-1(x/ζ1×x/ζ2)와 같은 관계를 가진다.

5. 수치 계산 결과

수치 결과 비교를 위해 ξ-영역에서 극좌표계로 변환하는 기존의 방법과 가상 매핑을 적용한 ζ-영역에서 극좌표계로 변환한 방법을 비교하였으며, 전자는 PCT(Polar coordinate transformation), 후자는 PMM이라고 표기하였다. 그리고 서로 다른 기하학적 특성을 가지는 두 개의 고차요소에 대해 수치 계산을 수행하였다. 첫 번째 요소는 높은 가로세로 비(Aspect ratio)와 마주하는 변들에 대해 큰 전단방향 변형을 가지고 있고 한 방향에 대해 높은 곡률을 가지고 있고, 두 번째는 KCS(KRISO container ship) 선형 중 상대적으로 큰 뒤틀림(Twisting) 변형을 가지는 고차요소를 선정하였다(Table 1). 또한, 수치 계산은 특이점이 ξ-영역에서 특이점이 다양한 위치들에 놓여진 경우들에 대해 수행하였으며, 특이점 위치들은 Table 2에 표기하였다.

Position vectors of test elements

Locations of singular points

PCT와 PMM모두 반경-방향 구적점은 5개로 고정시켰으며, 각-방향 구적점들은 2개에서 20개까지 1개씩 변화시켜가며 수치계산을 수행하였고, PCT에서 반경방향과 각방향 모두 구적점을 300개씩 사용한 결과를 기준으로 잡아 상대 오차(Relative errors)를 계산하였다(Figs. 4-5). 상대오차는 다음과 같이 계산할 수 있다.

Fig. 4

Numerical results for 1st element: weakly and strongly singular cases

Fig. 5

Numerical results for 2nd element: weakly and strongly singular cases

(35) Relativeerror=|I300×300I5×NGPθI300×300|

여기서 NGPθ는 각방향 가우스 구적점 수를 의미한다.

두 방법의 비교에서 가장 두드러진 점은 각-방향 구적점의 변화에 따라 PMM의 결과치들은 PCT에 비해 상대적으로 점진적 감소(Monotonically decreasing) 경향을 가진다는 사실이다. 이는 수치 계산 결과의 수렴성과 강인성 및 구적점 사용 시의 적절한 구적점 수를 판단하는데 있어 매우 중요한 요소이다. 먼저, 첫 번째 요소의 경우 두 번째 요소에 비해 가상 매핑의 영향이 더 크게 나타났는데, 이는 앞서 언급한 사실과 같이 두 가상 매핑 방법들이 가지는 물리적 의미들과 첫 번째로 선정한 요소의 기하학적 특성이 상당히 잘 부합한다는 데에 있다. 첫 번째 요소는 큰 전단 변형과 높은 가로세로비를 가지고 있고 이는 두 가상 매핑 방법에 각각 대응한다. 그리고 약한 특이성보다는 강한 특이성 문제에서 PMM의 효능이 훨씬 더 크게 나타났다는 점이 주목할 만한데, 이는 약한 특이성과 강한 특이성 문제가 서로 다른 특성들을 가지는 데에서 찾을 수 있다. 약한 특이성의 경우는 근사해가 아닌 엄밀한 적분해가 존재하지만 강한 특이성 문제는 코시-주요값(Cauchy-principal value)문제로, 수치 적분을 통해 점근적인 근사값을 얻을 수 있으며, 특이점 기준으로 반경방향으로 얼마나 균일하게 떨어진 위치들에 구적점들이 놓이는지가 수치 계산 정확도와 수렴율을 결정할 수 있기 때문에 PMM을 적용하면 더 좋은 결과를 얻을 수가 있다. 실제 PMM이 적용된 이후의 매핑 영역에서는 반경방향 구적점들이 각방향으로 둘렀을 때 특이점 부근의 구적점들 배치가 원 모양에 가깝게 적용되므로, 코시-주요값문제의 경우 높은 수렴율을 보여줄 수가 있게 된다.

두 번째 요소는 큰 뒤틀림이 있는 요소이며, 첫 번째 요소의 결과와 같이 PMM의 상대오차가 점진적으로 감소하는 특성을 보여주지만, 전체적인 상대오차 감소율은 PCT와 비슷한 수준이다. 이는 두 번째 요소가 갖고 있는 뒤틀림 현상은 일차적인 매핑에서 발생하는 현상들 중 선형 왜곡이 아닌 비선형적 왜곡현상이기 때문이며 곡면 요소의 경우는 이러한 비선형적 왜곡현상이 크게 나타날 수 있다. 물론, 선형 보간함수를 사용한 평면 사각 요소 문제에서는 2x/ξ1ξ1와 같은 한 방향에 대한 2차 공변 미분항은 항상 0이지만, 2x/ξ1ξ2와 같은 이중선형(Bilinear)항은 0이 아닐 수 있으므로 역시 비선형 왜곡 현상을 가질 수 있다. 이러한 비선형적 왜곡 현상들이 상대적으로 큰 요소들에서는 가상 매핑된 영역에서의 공변 기저벡터들의 변화가 비선형적 왜곡 특성들을 더 악화시킬 수 있기에 좋지 않은 결과를 초래할 수 있다. 하지만 두 번째 요소와 같은 기하학적으로 복잡한 요소들은 격자 생성 시 잘 쓰이지 않고, 이런 요소들을 곧바로 사용하면 특이점 적분 뿐 아니라 위치 벡터들이나 속도 포텐셜의 다항식 근사에서도 큰 오차를 유발할 수 있기 때문에 격자 생성 시 피하는 것이 좋다. 반대로, 첫 번째 요소는 일반적인 격자에서 흔하게 발생할 수 있는 선형 왜곡 현상을 극대화시킨 모델로 PMM의 선형 왜곡 완화 메커니즘만으로도 대부분 고차요소에서는 높은 정확도와 수렴율을 가지는 수치적분이 가능하다는 것을 보여주며, 격자 생성 시 가로세로비의 크기와 전단 방향 등 선형 왜곡 특성들을 고려하지 않아도 되므로 시간 차분에 따라 격자가 변하는 시간 영역 수치 시뮬레이션 문제에서 보다 큰 이점을 가질 것으로 예상된다.

6. 결론 및 고찰

본 연구에서는 연속체의 경계치 문제의 수치적 해법으로 주로 사용되고 있는 고차경계요소법의 특이점 적분을 위한 가상 매핑 방법이 제안되었다. 가상 매핑 방법은 경계요소가 3차원 영역에서 국부 요소로 매핑될 때 발생하는 물리량 왜곡 현상의 특성을 반영한 방법으로, 이 중 선형 왜곡 거동을 완화시켜주는 역할을 한다. 첫 번째 접선 변형에서는 공변 기저벡터들이 직교하도록 만들어주고, 두 번째 수직 변형에서는 서로 다른 공변 기저벡터들의 크기가 같도록 만들어주어 랜킨 소스의 테일러 급수 전개로부터 얻을 수 있는 두 조건들을 만족시킬 수 있도록 해준다. 수치 계산은 선형 왜곡 특성이 큰 첫 번째 요소와 비선형 왜곡 특성이 큰 두 번째 요소에 대해 수행되었고, 일반적인 고차요소를 대변하는 첫 번째 요소에서는 높은 정확도와 수렴율을 보여주었다. 반면, 기하학적으로 복잡한 두 번째 요소에서는 가상 매핑의 영향이 크지 않았는데, 이는 비선형 왜곡 현상에서 기인한 것으로 판단되며 비선형적 특이 거동을 완화시킬 수 있는 추가적인 가상 매핑 방법의 개발이 요구된다.

Acknowledgements

이 논문은 부산대학교 기본연구지원사업(2년)에 의하여 연구되었음.

References

Guiggiani M, Gigante A. 1990;A General Algorithm for Multidimensional Cauchy Principal Value Integrals in the Boundary Element Method. Journal of Applied Mechanics 57(4):906–915. http://www.doi.org/10.1115/1.2897660.
Guiggiani M, Krishnasamy G, Rudolphi TJ, Rizzo FJ. 1992;A General Algorithm for the Numerical Solution of Hypersingular Boundary Integral Equations. Journal of Applied Mechanics 59(3):604–614. http://www.doi.org/10.1115/1.2893766.
Liu YH, Kim CH, Lu XS. 1991;Comparison of Higher-order Boundary Element and Constant Panel Methods for Hydrodynamic Loadings. International Journal of Offshore and Polar Engineering 1(1):8–17.
Lv JH, Feng XT, Yan F, Jiang Q. 2015;Efficient Evaluation of Integrals with Kernel 1/rx for Quadrilateral Elements with Rrregular Shape. Engineering Analysis with Boundary Elements 61:33–41. https://doi.org/10.1016/j.enganabound.2015.06.010.
Rong J, Wen L, Xiao J. 2014;Efficiency Improvement of the Polar Coordinate Transformation for Evaluating BEM Singular Integrals on Curved Elements. Engineering Analysis with Boundary Elements 38:83–93. https://doi.org/10.1016/j.enganabound.2013.10.014.

Article information Continued

Fig. 1

Mapping from 3-D curved panel to local element and distortion of singular behavior near the singularity

Fig. 2

Pseudo mapping method: Tangential deformation

Fig. 3

Pseudo mapping method: Normal deformation

Fig. 4

Numerical results for 1st element: weakly and strongly singular cases

Fig. 5

Numerical results for 2nd element: weakly and strongly singular cases

Table 1

Position vectors of test elements

Element type 1st element 2nd element
Node number 1 (0, −2, 0) (9.60383570, 0, −2.30673860)
2 (1, 0, 0) (10.6937260, −0.32298802, 2.17538260)
3 (1, 2, 0) (12.016290, −0.57100546, −2.34276250)
4 (0.92387953, 3, 0.38268343) (11.3972740, −1.04141720, −1.84616940)
5 (0.70710678, 4, 0.70710678) (10.6447760, −1.54068710, −1.55561950)
6 (0.70710678, 2, 0.70710678) (9.40466430, −0.91600309, −1.51056070)
7 (0.70710678, 0, 0.70710678) (8.51099060, 0, −1.55062450)
8 (0.92387953, −1, 0.38268343) (9.08797830, 0, −1.88489470)
9 (0.92387953, 1, 0.38268343) (10.0882940, −0.65084745, −1.73879260)

Table 2

Locations of singular points

Locations Singular points (ξ̃1, ξ̃2)
Inner (0, 0)
Edges (0, 1), (1, 0)
Vertices (1, 1), (1, −1)