Journal Search Engine
Search Advanced Search Adode Reader(link)
Download PDF Export Citaion korean bibliography PMC previewer
ISSN : 1229-3059(Print)
ISSN : 2287-2302(Online)
Journal of the Computational Structural Engineering Institute of Korea
Vol.32 No.4 pp.215-222

DOI : https://doi.org/10.7734/COSEIK.2019.32.4.215

Analysis of Interface Problem using the MLS Difference Method with Interface Condition Embedment

Young-Cheol Yoon1†
1Department of Civil Engineering, Myongji College, Seoul, 03656, Korea
Corresponding author: Tel: +82-2-300-1135; E-mail: ycyoon@mjc.ac.kr
April 22, 2019 July 17, 2019 July 18, 2019

Abstract


The heat conduction problem with discontinuous material coefficients generally consists of the conservative equation, boundary condition, and interface condition, which should be additionally satisfied in the solution procedure. This feature often makes the development of new numerical schemes difficult as it induces a layered singularity in the solution fields; thus, a special approximation is required to capture the singular behavior. In addition to the approximation, the construction of a total system of equations is challenging. In this study, a wedge function is devised for enriching the approximation, and the interface condition itself is embedded in the moving least squares(MLS) derivative approximation to consistently satisfy the interface condition. The heat conduction problem is then discretized in a strong form using the developed derivative approximation, which is named as the interface immersed MLS difference method. This method is able to efficiently provide a numerical solution for such interface problems avoiding both numerical quadrature as well as extra difference equations related to the interface condition enforcement. Numerical experiments proved that the developed numerical method was highly accurate and computationally efficient at solving the heat conduction problem with interfacial jump as well as the problem with a geometrically induced interfacial singularity.



계면경계조건이 매입된 이동최소제곱 차분법을 이용한 계면경계문제 해석

윤 영 철1†
1명지전문대학 토목공학과

초록


복합재료의 열전달 문제는 일반적으로 만족시켜야 하는 보존방정식과 경계조건 외에 추가적으로 만족시켜야 하는 계면경 계조건의 존재로 인해 새로운 수치기법의 개발에 어려움이 있다. 계면경계조건이 미분방정식의 해에 불연속성을 유발시키기 때문에 이것을 적절하게 처리할 수 있는 특수한 함수의 도입이 필요하며, 이산화를 통한 계 방정식의 구성도 쉽지 않다. 본 논문에서는 계면경계의 불연속성을 모사하는 특수함수를 포함하면서 계면경계조건을 항상 만족시킬 수 있도록 계면경계식 자체를 매입한 미분근사식을 제안하고, 불연속 재료상수를 갖는 열전달 문제를 무요소 강형식으로 이산화한 이동최소제곱 차분법을 제시한다. 개발된 수치기법은 기존의 수치기법들과 달리 수치적분과 계면경계조건을 만족시키기 위한 별도의 구속 방정식이 필요없으며, 빠르고 정확하게 이종재료 열전달 문제의 수치해를 구해준다. 개발된 수치기법으로 다양한 복합재료 열전달 문제를 해석하고 오차의 수렴률을 조사한 결과, 높은 정확성과 계산 효율성을 갖는다는 것을 확인할 수 있었으며, 특 히, 계면경계가 기하학적 특이성을 나타내는 문제에서도 우수한 성능을 발휘하는 것을 보였다.



    National Research Foundation of Korea
    NRF-2017R1A1A1A0 5001196

    1. 서 론

    이종재료의 열전달 문제는 계면경계(interface)를 따라 해의 미분이 불연속한 특징을 나타내며, 보존방정식과 경계조건 외에 열유속(heat flux)의 점프(jump)에 대한 계면경계조건을 추가로 만족시켜 주어야 한다. 따라서 일반적인 수치해석 방법 으로는 정밀한 해석이 어렵기 때문에 불연속성을 처리할 수 있는 근사식을 갖는 확장된 유한요소법(XFEM; Moes et al., 1999)이나, 계면경계의 기하학적 특성을 그리드 구성 시 반영 하여 차분식을 수정한 유한차분법 기반의 IIM(Immesed Interface Method; LeVeque and Li, 1994)과 같이 특수 하게 고안된 수치기법의 적용이 요구된다. 특히, 유한차분법의 경우, Osher와 Sethian(1988)이 선단추적기법(Front Tracking Method)이라는 이름으로 움직이는 경계를 추적할 수 있는 기법을 제안한 이후로 Tu와 Peskin(1992)은 지배 미분방정 식의 열원(heat source)항을 변형시켜 계면경계조건을 만족 시킬 수 있는 Immersed Boundary Method(IBM)을 제시 하였고, 전술한 바와 같이 IIM(LeVeque and Li, 1994)이 여러 가지 이동경계문제를 수치 모사하였으며, Hou와 Liu (2005)는 불연속 재료상수를 갖는 다양한 타원형 미분방정식을 해석할 수 있는 기법을 제시하였다. 무요소법의 일종으로 약정 식화(weak formulation)를 기반으로 한 Element-Free Galerkin(EFG)법과 최종 평형방정식에 추가적인 구속조건을 부여할 수 있는 penalty 방법을 조합하여 이종재료의 응력해석 문제를 풀 수 있는 수치기법도 제시된 바 있다(Cordes and Moran, 1996).

    그러나 이러한 수치기법들도 임의로 위치하는 불연속면으로 인해 근사함수의 구성이 상당한 제약을 받고 계면경계 모델링 으로 인해 최종적인 계방정식의 크기가 증가되는 단점이 지적 되어 왔다. 한편, 요소망없이 절점만 사용하는 무요소법에서 빈번하게 사용되는 이동최소제곱법은 근사함수의 구성에 상당한 유연성을 갖고 있어서 추가적인 구속조건을 만족시키는 것이 용이하고, 특히, 강정식화(strong formulation)를 통해 절점 단위에서 강형식(strong form)을 직접 이산화하고 풀면 이러 한 특이성을 갖는 문제에 유리하다는 것이 알려져 있다(Yoon et al., 2007;Yoon and Song, 2014a;2014b). 강정식화 된 무요소법의 장점을 잘 활용하면, 기존의 유한요소법이나 유한차분법에서 나타난 요소망이나 그리드에 의한 제약으로 인한 단점을 효과적으로 극복할 수 있다. Yoon 등(2007)은 규칙적으로 분포된 절점을 바탕으로 이동최소제곱법을 이용하여 불연속 함수를 포함하는 차분식을 구성하고 이종재료의 열전 달문제를 정밀하게 해석했지만, 계면경계의 모델링으로 인해 계 방정식의 크기가 증가되어 효율성이 감소되는 문제점이 있었다. 본 연구에서는 계면경계조건을 근사함수에 직접 매입하여 미분근사식이 계면경계조건을 자동으로 만족하도록 함으로써 계 방정식 구성 시 자유도의 증가가 필요없는 수치기법을 제시하고 다양한 이종재료의 열전달 문제를 해석하여 계산 효율성과 정확 성을 검증한다.

    2. 복합재료 열전도 방정식 해석을 위한 계면경계조건이 매입된 근사함수

    열전달 방정식에서 시간에 대한 항을 제외하고 계면경계 Γ에 대한 계면경계조건이 추가된 이종재료의 열전달 문제의 지배 방정식은 다음과 같다.

    · ( κ u ) = f in Ω / Γ
    (1)

    u = u ¯ on Ω
    (2)

    [ u ] Γ = 0 on Γ
    (3)

    [ κ u n ] Γ = h ¯ on Γ
    (4)

    여기서, f는 열원이고, κ는 열전달계수인데 이종재료에 대해 서로 다른 값을 갖는다. 계면경계조건은 식 (3)의 해의 연속 조건과 식 (4)의 계면경계를 통과하는 κ를 포함한 수직방향 열유속의 점프에 대한 조건식으로 주어진다.

    연속함수에 대한 이동최소제곱 차분법의 근사식은 Taylor 전개에 근거하여 다음과 같이 주어진다(Yoon et al., 2007).

    u R ( x , x ¯ ) = p m T ( x , x ¯ ) a ( x ¯ )
    (5)

    여기서, p m T ( x , x ¯ ) m차까지의 다항식 항들로 구성된 기저함 수이고, a(x)는 Taylor 전개식의 미분계수를 포함하는 미지의 벡터이다. 즉, 식 (5)는 x 근방에서 전개한 Taylor 다항식 으로서 임의의 위치 x에서의 u(x) 에 대한 근사값을 제공한다. 이종재료의 열전달 문제를 해석하기 위한 계면경계 근방에서의 미분불연속을 갖는 확장된 근사식은 다음과 같이 쓸 수 있다.

    u S ( x , x ¯ ) = p m T ( x , x ¯ ) a ( x ¯ ) + b s ( x ¯ ) b Γ ( x , x ¯ )
    (6)

    여기서, b Γ ( x , x ¯ ) 는 미분불연속을 갖는 쐐기함수이고, bs (x) 는 쐐기의 각도 또는 미분불연속의 크기를 결정하는 계수이다. 식 (6)은 쐐기함수를 포함하고 있어서 한번 미분을 하면 점프를 갖지만 근사식 자체는 C0 연속이기 때문에 식 (3)을 만족시킨다. 즉, u S ( x , x ¯ ) [ u S ] Γ = 0 이고 [ u S n ] Γ 0 이다. bs (x) 를 계산 하기 위해 식 (4)를 다음과 같이 쓸 수 있다.

    κ Γ [ u n ] Γ + [ κ ] Γ u n Γ = h ¯
    (7)

    여기서, <⋅>ΓΓ에서의 평균값을 의미한다. [κ]Γ는 계면 경계에서 열전달계수의 차이, <κ>Γ는 계면경계에서 열전달계 수의 평균값이다. 식 (7)에 식 (6)을 적용하면, p m T ( x ¯ Γ , x ¯ ) 는 불연속성이 없으므로 u n Γ = p m T ( x ¯ Γ , x ¯ ) n a ( x ¯ ) 이다. 이때, 다항식 벡터를 계면경계 위로 투영한 후 수직방향에 대해 미분한 값은 다음과 같다.

    p m T ( x ¯ Γ , x ¯ ) n = ( 0 , n x , n y , ( x Γ ¯ x ¯ ρ ) , ( y Γ ¯ y ¯ ρ ) , ( x Γ ¯ x ¯ ρ ) n x , ( y Γ ¯ y ¯ ρ ) n x + ( x Γ ¯ x ¯ ρ ) n y , ( y Γ ¯ y ¯ ) ρ ) n y )
    (8)

    여기서, xΓΓ를 이산화한 접선평면 위로의 x 에 대한 투영 점이고 이산화된 계면경계 상에서 수직방향으로 x 와 가장 가까운 점이 된다. 식 (8)은 x 의 투영점 xΓ가 적용된 후로는 더 이상 변수가 아닌 상수값이 되기 때문에 미분의 영향을 받지 않는다. 이제 식 (7)로부터 계산된 bs (x) 를 식 (6)에 다시 대 입하면 다음과 같이 식 (4)의 계면경계조건을 항상 만족시키 는 확장된 근사식을 얻는다.

    u S ( x , x ¯ ) = p Γ T ( x , x ¯ ) a ( x ¯ ) + h ¯ 2 κ Γ b Γ ( x , x ¯ )
    (9)

    식 (4)를 만족시키는 과정에서 변형된 다항식 기저함수 p Γ T ( x , x ¯ ) 는 다음과 같다.

    p Γ T ( x , x ¯ ) = p m T ( x , x ¯ ) [ κ ] Γ 2 κ Γ p m T ( x ¯ Γ , x ¯ ) b Γ ( x , x ¯ )
    (10)

    이동최소제곱 차분법에서는 절점이 적절하게 분포된 해석영 역에서 이동최소제곱법을 이용하여 미지의 벡터 a(x)를 구하면 실제 미분계산 없이도 m차까지의 미분근사식을 얻을 수 있다. 이제, a(x)를 구하기 위해 아래와 같은 가중잔차식을 정의한다.

    J = I w ( x I x ¯ ρ ) [ u s ( x I , x ¯ ) u I ] 2
    (11)

    여기서, w ( x I x ¯ ρ ) 는 이동최소제곱 차분법에서 일반적으로 사용하는 가중함수이다. 잔차식에 대한 정지조건(stationary condition) J a = 0 로부터 미지의 벡터 a(x) 구한 후 xx로 치환하면, α 차 미분근사식을 다음과 같이 얻는다.

    D x α ( x , x ) = I Φ I α ( x ) u I + h ¯ 2 κ Γ ( D x α b Γ ( x , x ) I Φ I α ( x ) b Γ ( x I , x ) )
    (12)

    여기서, Φ I α ( x ) α 차 형상함수를 가리키며, α = ( α 1 , , α n ) 이고, n은 공간차수이다. 예를 들어, α = (1,1)은 2차원 문제 에서 2 x y 를 의미한다. 보다 자세한 유도과정은 Yoon 등 (2007) 또는 Yoon과 Song(2014a)을 참고할 수 있다. 식 (12)에서 유일한 미지계수는 uI 인데, 이것은 이산화된 지배 미분방정식에 대한 계 방정식을 풀어서 얻는다. 기존의 수치기 법들처럼 근사식을 유도한 후 그 식을 직접 미분하면 계산이 매우 복잡한 반면, 식 (12)는 실제로 D x α b Γ ( x , x ) 만 계산하면 되기 때문에 계산이 매우 간단하다. 한편, 식 (12)의 형상함수 Φ I α ( x ) 는 식 (5)로부터 계산된 일반적인 형상함수와 달리 계면 경계조건 식 (4)를 만족시키기 위한 항들이 내재적으로 포함 되어 있는 새로운 형태를 갖는다. 이와 같이 수정된 형상함수는 확장된 기저함수의 미분을 포함하며 다음과 같이 표현된다.

    Φ I α ( x ) = ( D x α p Γ T ( x , x ) ) M 1 ( x ) B ( x )
    (13)

    여기서, M(x) 과 B (x) 는 ‘모멘트 행렬’과 ‘B 행렬’로 불리는 무 요소법에서 나타나는 전형적인 행렬식이며, 상세한 형태는 Yoon 등(2007) 또는 Yoon과 Song(2014a)을 참고할 수 있다. 예를 들어, 형상함수에 대한 x방향 1차 미분근사와 x방향 2차 미분근사는 각각 Φ I ( 1 , 0 ) ( x ) = ( D x ( 1 , 0 ) p Γ T ( x , x ) ) M 1 ( x ) B ( x ) , Φ I ( 2 , 0 ) ( x ) = ( D x ( 2 , 0 ) p Γ T ( x , x ) ) M 1 ( x ) B ( x ) 로 쓸 수 있다. 결과 적으로 식 (12)가 계면경계조건을 이미 만족시키고 있기 때문에 계면경계조건에 대한 추가적인 고려 없이도 각 절점에서 보존 방정식에 대한 차분식을 구성하는 것만으로 계 방정식을 구성 할 수 있으며, 이 계 방정식을 풀어서 절점해인 미지계수 uI 를 얻는다.

    3. 지배 미분방정식에 대한 차분식의 구성

    열전달 방정식의 이산화를 위해 보존방정식 식 (1)과 경계 조건식 식 (2)에 대한 차분식을 구성해야 한다. 이산화 과정에서 고려하는 절점의 영향영역이 계면경계를 일부라도 포함하면 식 (12)를 사용하고, 그렇지 않은 경우에는 p m T ( x , x ¯ ) 로 구성한 기존의 미분근사식을 사용한다. 본 연구에서 제안한 근사식이 계면경계조건을 자동으로 만족시키기 때문에 계면경계조건에 대한 별도의 차분식은 필요없으며, 다음과 같은 차분식의 구성만으로 문제의 해석이 가능하다.

    κ ( x J ) Δ u ( x J ) + κ ( x J ) · u ( x J ) = f ( x J ) , x J Λ Ω
    (14)

    u ( x K ) = u ¯ ( x K ) , x K Λ Ω
    (15)

    여기서, ΛΩΛ6Ω는 각각 내부영역의 절점과 경계절점의 집합을 나타낸다. 계면경계 근방에서 라플라시안(laplacian)과 그라 디언트(gradient)는 식 (12)를 이용하여 계산한다. 차분식 구성시 미지계수인 uI 와 관련된 항들만 좌변에 남기고 나머지 기지값들을 모두 우변으로 넘기면 다음과 같은 차분방정식을 얻는다.

    κ ( I ( Φ I ( 2.0 ) ( x K ) + Φ I ( 2.0 ) ( x K ) ) u I ) + D x ( 1 , 0 ) κ I Φ I ( 1 , 0 ) ( x K ) u I + D x ( 0 , 1 ) κ I Φ I ( 0 , 1 ) ( x K ) u I = f ( x K ) κ h ¯ 2 κ Γ ( D x ( 2 , 0 ) b Γ ( x K , x K ) I Φ I ( 2 , 0 ) ( x K ) b Γ ( x I , x K ) + D x ( 0 , 2 ) b Γ ( x K , x K ) I Φ I ( 0 , 2 ) ( x K ) b Γ ( x I , x K ) ) h ¯ 2 κ Γ ( D x ( 1 , 0 ) κ { D x ( 1 , 0 ) b Γ ( x K , x K ) I Φ I ( 1 , 0 ) ( x K ) b Γ ( x I , x K ) } + D x ( 0 , 1 ) κ { D x ( 0 , 1 ) b Γ ( x K , x K ) I Φ I ( 0 , 1 ) ( x K ) b Γ ( x I , x K ) } )
    (16)

    계산의 편의상 전체 계 방정식은 보존방정식과 경계조건식에 대한 차분식을 식의 종류에 관계없이 절점 순서대로 조립 (assemble)한다. 본 연구에서는 계 방정식의 계수행렬의 조건수(condition number)가 크지 않은 것을 확인하였으며, 이것은 라플라시안 항이 계수행렬의 타원성(ellipticity)을 잘 유지시켜 주기 때문인 것으로 판단된다. 한편, 계면경계조건을 만족시키는 근사함수를 사용하더라도 계 방정식의 구성 시 고려 하는 절점에 대한 계면경계로의 투영점을 찾아야 하므로 계면 경계의 이산화는 필요하다. 그러나 계면경계 상에 부여된 미지 계수가 없기 때문에 계 방정식의 크기는 증가하지 않는다. 또한, 계면경계를 모델링하는 정밀도가 보존방정식에 대한 차분식의 정확도에 영향을 미치므로 절점분포의 밀도를 반영하고 계면 경계의 기하학적 형상을 잘 묘사할 수 있도록 계면경계를 적절한 간격으로 분할해 주어야 한다.

    4. 수치예제

    4.1 원형 계면경계를 갖는 열전달 문제

    본 절에서는 원형 매입물을 갖는 이종재료의 열전달 문제에 대한 해석을 수행하여 제안된 수치기법의 정확성을 검증했다. 열전달 계수가 식 (18)과 같이 주어질 때, 이론해는 식 (17)과 같이 정해진다.

    u ( x , y ) = { r 2 , r 2 < 1 / 4 1 4 ( 1 1 8 c 1 c ) + 1 c ( r 4 2 + r 2 + 1 10 l o g ( 2 r 2 ) ) , r 2 1 / 4
    (17)

    k ( x , y ) = { k = r 2 + 1 , r 2 < 1 / 4 k + = c , r 2 1 / 4
    (18)

    여기서, r = ( x x c ) 2 + ( y y c ) 2 이고, (xc, yc) 는 매입물 중심의 좌표이다. 해석 시 c = 100을 적용하여 열전도계수의 차가 최대 100배 가까이 나도록 설정했다. 이 경우 많은 수치 기법에서 미분방정식의 가해성(solvability)이 상당히 나빠 지기 때문에 계면경계조건의 정확한 처리없이 정확성을 확보 할 수 없다.

    이동최소제곱 차분법은 적분식을 사용하지 않기 때문에 L2 오차 대신 L 오차를 사용하여 수렴률을 조사한다. L2 오차가 해석영역 전체에 대한 오차를 대표하는데 반해, L 오차는 국부 적인 영역의 오차를 대표하기 때문에 특이성을 갖는 문제에서는 수렴이 상대적으로 더 어렵다. 수렴률 조사를 위해 400개, 1,600개, 6,400개의 절점을 갖는 모델에 대해 해석을 수행했다. Fig. 1은 온도와 온도 그라디언트에 대한 L 상대오차의 수렴률을 나타내는데, 온도와 그라디언트에 대해 각각 1과 2 이상의 기울기를 나타내고 있는 것을 통해 본 연구에서 제시한 수치기법이 2차 정확도를 확보하는 것을 알 수 있다. 계면경계 문제에서는 해의 미분에 존재하는 불연속성으로 인해 2차의 정확도를 얻기 어렵기 때문에 개발된 수치기법이 불연속성을 매우 정확하게 처리한다는 것을 알 수 있다. Fig. 2는 6,400개 절점 모델로 계산한 온도분포이고, 원형 매입물 경계를 따라 쐐기형상의 해가 잘 묘사되는 것을 볼 수 있다. Fig. 3에는 온도 그라디언트의 x방향 성분을 도시하였으며, 원형 매입물의 경 계에서 점프가 잘 표현되고 있다. 불연속의 크기가 상당히 큼 에도 불구하고 점프가 뭉개지는 smearing 현상이나 떨림 (oscillation) 현상도 나타나지 않는데, 이것은 개발된 미분근 사식이 계면경계조건을 정확히 만족시킨다는 것을 간접적으로 보여준다.

    4.2 사각형 계면경계를 갖는 열전달 문제

    본 절에서는 사각형 모서리를 갖는 계면경계의 기하학적인 효과가 불연속 재료상수를 갖는 열전달 방정식의 해석에 미치는 영향을 조사한다. 실제로 계면경계조건 만족을 위해 접선평면 위로의 투영 개념을 사용하는 미분근사식의 성능을 함께 검증 하게 되며, 앞 절의 부드러운 계면경계를 갖는 경우와 비교를 통해 부드럽지 않은 계면경계 형상을 갖는 문제에 대한 수치기 법의 강건성을 확인할 수 있다. 사각형 모양 계면경계로 둘러 싸인 영역을 Ω+, 나머지 사각형 영역에 둘러싸인 영역을 Ω- 로 가정할 때, 이 문제에 대한 이론해는 다음과 같다.

    u ( x ) = { r 2 + ( x 1 2 ) ( x 3 2 ) ( y 1 2 ) ( y 3 2 ) , i n Ω r 2 , i n Ω +
    (19)

    여기서, 불연속성을 갖는 재료상수는 다음과 같이 가정한다.

    κ = { 1 , i n Ω 10 , i n Ω +
    (20)

    Fig. 4는 400, 1,600, 6,400개의 일정하게 분포된 절점을 갖는 모델들을 사용하여 온도와 온도 그라디언트에 대해 계산 한 L 상대오차의 수렴률을 보여준다. 온도와 온도 그라디언 트에 대한 수렴률 모두 이론적 최적 수렴률인 1과 2 보다 작은 데, 이것은 해의 미분이 갖고 있는 불연속성 외에 사각형 계면 경계의 모서리에 의한 기하학적 형상이 추가적인 특이성을 유 발하기 때문이다. Figs. 5(a)~(c)는 6,400개 절점 모델로 계산한 온도와 x, y방향 온도 그라디언트를 3차원으로 표현한 것으로 계면경계를 따라 쐐기형상이 잘 묘사되는 것과 계면경 계 위치에서 smearing 현상없이 미분의 불연속성이 잘 표현되 는 것을 확인할 수 있다.

    4.3 쐐기형상의 계면경계를 갖는 열전달 문제

    본 절에서는 사각형 모서리보다 더 날카로운 형상을 갖는 계면경계를 갖는 문제를 통해 계면경계 형상에 대한 기하학적 영향을 좀 더 살펴본다. 쐐기모양의 계면경계는 해석영역의 중 심에 그 꼭지점을 두고 있으며 쐐기의 내부각도는 45°이다. 쐐기모양 계면경계로 둘러싸인 영역을 Ω+, 그리고 나머지 해석 영역에 둘러싸인 영역을 Ω-로 가정할 때, 이 문제의 이론해는 다음과 같다.

    u ( x , y ) = { r + c ( 1 s ) ( x x c ) ( y y c ) ) ( 1 1 s ( x x c ) ( y y c ) ) , i n Ω r , i n Ω +
    (21)

    여기서, cs는 임의의 값을 갖는 상수이고, 불연속 재료상수는 앞 절의 예제와 동일하게 가정했다. 본 연구에서는 계면경계가 절점들과 상관없이 임의의 위치에 위치할 수 있다.

    Fig. 6은 수치모델의 내부와 경계에 6,400개의 절점들을 일정하게 분포시킨 해석모델을 보여준다. Fig. 7은 400, 1,600, 6,400개의 절점을 사용한 모델로부터 계산한 L 상대 오차의 수렴률이다. 이 경우에도 계면경계의 기하학적 형상이 갖는 특이성으로 인해 온도와 온도 그라디언트에 대한 수렴률 모두 이론적 최적 수렴률 보다 작게 나타나고 있다. 그러나 이 동최소제곱 차분법은 특이성을 갖는 문제에 대해서도 여전히 좋은 수렴률을 유지하고 있다는 점을 주목할 필요가 있다. Figs. 8(a)~(c)는 6,400개 절점 모델로 계산한 온도와 x, y방향 온도 그라디언트를 3차원으로 표현한 것이다. 수치해가 이론해와 잘 일치하기 때문에 계면경계를 따라 쐐기형상이 잘 묘사되고 있으며, 계면경계 위치에서 미분의 점프도 잘 표현되고 있다. 불연속의 크기가 상당히 큼에도 불구하고 smearing 현상이나 떨림 현상은 나타나지 않았다. 결과적으로 새롭게 개발된 이동최소제곱 차분법이 기하학적 형상의 특이 성을 포함한 계면경계의 특이성을 효과적이고 정확하게 처리 하고 있음을 알 수 있다.

    4. 결 론

    본 논문에서는 계면경계조건이 근사함수 안으로 매입되어 계면경계조건을 자동으로 만족시키는 근사미분식을 제시하고 그것을 이용하여 계면경계 문제를 효율적으로 강정식화하고 해석할 수 있는 이동최소제곱 차분법을 제시하였다. 개발된 수 치기법은 접선평면으로의 투영을 이용한 쐐기함수를 사용하는 확장된 근사미분식을 이용해 절점에 의존하지 않고 계면경계를 자유롭게 모델링할 수 있기 때문에 계면경계가 갖는 미분불연속 특성을 효과적으로 묘사할 수 있다. 이종재료의 열전도 방정식은 규칙적으로 분포된 절점들 위에서 구성된 차분식을 계 방정식을 통해 강형식으로 이산화되었다. 계면경계에 대한 추가적인 방 정식을 만들 필요가 없기 때문에 계 방정식의 크기는 절점의 개수와 일치하게 되어 효과적인 강정식화가 가능하다. 이와 같은 방식으로 구성된 계 방정식은 약 정식화를 사용하는 수치기법 에서 필수적인 수치적분이나 계면경계조건과 관련된 미지계수의 추가적인 도입이 필요없어 계산속도가 빠르다. 기존 수치방법 들이 계면경계조건의 만족을 위해 미지계수를 추가로 도입하고 그 결과로 전체 계방정식의 크기가 증가되었던 것과 비교하면 수치 알고리즘 개발 측면에서 상당히 개선된 부분이다.

    다양한 이종재료의 열전달 문제에 대한 수치해석을 통해 L 상대오차의 수렴률을 조사하여 개발된 수치기법이 2차 정 확도를 갖는 것을 확인하였다. 개발된 수치기법이 부드러운 형상을 갖는 계면경계 문제에서 이론적 최적값에 상응하는 수렴률을 보여줄 뿐만 아니라, 날카로운 형상의 계면경계를 갖는 문제에 서도 우수한 성능을 발휘한다는 것을 보였다. 날카로운 계면경 계는 해의 미분이 갖는 불연속성 외에 계면경계의 기하학적 형상으로 인한 특이성이 추가되기 때문에 수치적으로 정확한 계면경계의 처리없이는 높은 정확성 확보가 불가능한데, 본 연구에서 제시한 수치기법은 이러한 계면경계의 복합적인 특이 성을 잘 묘사하고, 지배 미분방정식도 정확하고 효율적으로 풀 수 있다는 것을 보였다. 본 연구에서 제안된 수치기법은 향후 추가적인 연구를 통해 계면경계의 형상이 시간에 따라 변화하는 이동경계 문제와 균열이 전파하는 불연속 고체역학 문제 그리고 재료 비선형을 고려한 대변형 문제에도 적용할 수 있을 것으로 기대된다.

    감사의 글

    이 논문은 정부(미래창조과학부)의 재원으로 한국연구재단의 지원을 받아 수행된 기초연구사업임(NRF-2017R1A1A1A0 5001196).

    Figure

    COSEIK-32-4-215_F1.gif

    Convergence rates of the relative L errors for the heat conduction problem with a circular interface

    COSEIK-32-4-215_F2.gif

    Surface plot for the computed temperature with a contour(circular interface case)

    COSEIK-32-4-215_F3.gif

    Surface plot for the x-directional temperature gradient(∂u∂x)

    COSEIK-32-4-215_F4.gif

    Convergence rates of the relative L errors for the heat conduction problem with a square interface

    COSEIK-32-4-215_F5.gif

    Surface plots for the computed temperature and temperature gradients obtained for 6400 node model (square interface case)

    COSEIK-32-4-215_F6.gif

    6,400 node model for the heat conduction problem with a non-smooth interface

    COSEIK-32-4-215_F7.gif

    Convergence rates of the relative L errors for the heat conduction problem with a wedge-shaped interface

    COSEIK-32-4-215_F8.gif

    Surface plots for the computed temperature and temperature gradients obtained for 6,400 node model(wedge-shaped interface case)

    Table

    Reference

    1. Cordes, L.W. , Moran, B. (1996) Treatment of Material Discontinuity in the Element-Free Galerkin Method, Comput. Methods Appl. Mech. & Eng., 139(1-4), pp.75~89.
    2. Hou, S. , Liu, X.-D. (2005) A Numerical Method for Solving Variable Coefficient Elliptic Equation with Interfaces, J. Comput. Phys., 202(2), pp.411~445.
    3. LeVeque, R.J. , Li, Z. (1994) The Immersed Interface Method for Elliptic Equations with Discontinuous Coefficients and Singular Sources, SIAM J. Numer. Anal., 31(4), pp.1019~1044.
    4. Moes, N. , Dolbow, J. , Belytschko, T. (1999) A Finite Element Method for Crack Growth without Remeshing, Int. J. Numer. Methods Eng., 46(1), pp.131~150.
    5. Osher, S. , Sethian, J.A. (1988) Fronts Propagating with Curvature-dependent Speed: Algorithms based on Hamilton-Jacobi Formulations, J. Comput. Phys., 79(1), pp.12~49.
    6. Tu, C. , Peskin, C.S. (1992) Stability and Instability in the Computation of Flows with Moving Immersed Boundaries: A Comparison of Three Methods, SIAM J. SCi. Statist. Comput., 13(6), pp. 1361~1376.
    7. Yoon, Y.-C. , Kim, D.W. (2007) Heat Transfer Analysis of Bi-Material Problem with Interfacial Boundary Using Moving Least Squares Finite Difference Method, J. Comput. Struct. Eng. Inst. Korea, 20(6), pp.779~787.
    8. Yoon, Y.C. , Song, J.-H. (2014a) Extended Particle Difference Method for Weak and Strong Discontinuity Problems: Part I. Derivation of the Extended Particle Derivative Approximation for the Representation of Weak and Strong Discontinuities, Comput. Mech., 53(6), pp.1087~1103.
    9. Yoon, Y.C. , Song, J.-H. (2014b) Extended Particle Difference Method for Weak and Strong Discontinuity Problems: Part II. Formulations and Applications for Various Interfacial Singularity Problems, Comput. Mech., 53(6), pp.1105~1128.