KR20130138702A - 고점성 유체의 시뮬레이션 방법 - Google Patents

Abstract

본원에는 챔버 내의 고점성 유체를 시뮬레이션하는 전산화된 시뮬레이션 방법이 개시되어 있으며, 이 시뮬레이션 방법에서는 유체의 모델을 챔버의 모델 내에 세팅하고, 유동 계산을 수행한다. 상기 유동 계산에서는, 유체 모델이 접촉하는 챔버 모델의 벽의 접촉면에 대하여, 유체 모델의 슬립 속도를 특정 식에 의해 정한다.

Description

고점성 유체의 시뮬레이션 방법{METHOD FOR SIMULATING HIGH-VISCOSITY FLUID}

본 발명은 고점성 유체를 시뮬레이션하는 전산화된 방법에 관한 것이다.

최근에는, 밴버리 혼합기 등의 혼합기의 챔버 내에서 혼련되는 비가교 고무 또는 수지 복합체 등의 가소성 유체의 유동 상태를, 네비어-스톡스 방정식을 이용하여 계산 또는 시뮬레이션하는 방법을 비롯한 다양한 전산화된 시뮬레이션 방법이 제안되어 있다.

특허문헌 1 : JP-A-2011-27593

비특허문헌 1 : "Numerical and Experimental Study of Dispersive Mixing of Agglomerates" V. Collin, E. Peuvrel-Disdier et al.

공기 등의 저점성 유체의 흐름을 계산하는 경우, 유체가 유동하는 공간의 벽의 표면에서 유체의 유속은 제로로 설정될 수 있다. 그러나, 아직 가황처리되지 않은 비가교 고무 복합체 등의 고점성 유체의 흐름을 계산하는 경우, 유체가 유동하는 공간의 벽의 표면에서, 유체는 유동 방향으로 소정값의 유속을 가질 수 있다. 따라서, 유체는 벽에서 미끄러진다. 이러한 고점성 유체에 대해 전산화된 시뮬레이션을 작성하는 경우, 벽에서의 슬립 현상 등을 고려해야만 한다.

지금까지는, 벽면에서의 유속(이하, 슬립 속도)은 벽면에서의 전단 응력의 함수로서 정의된다. 슬립 속도가 선형인 경우, 네비어 법칙에 따른 전단 응력 τw은 하기의 식(1)에 의해 주어진다:

τw = F_slip[v_slip-v_wall] (1)

슬립 속도가 비선형인 경우, 전단 응력 τw은 하기의 식(2)에 의해 주어진다:

τw = F_slip[v_slip-v_wall]|v_slip-v_wall|^{eslip -1} (2)

여기서

"v_slip"은 벽면에 있어서 유체의 벽면과 평행한 방향의 속도이고,

"v_wall"은 벽면의 이동 속도에 있어서 벽면과 평행한 방향의 성분이며,

"F_slip"은 유저에 의해 정해지는 상수이고, 그리고

"e_slip"은 유저에 의해 정해지는 상수이다.

상수 "F_slip"으로 설정된 값은 해당 유체의 고유의 값이며, 슬립 유발의 용이성과 관련이 있다. 대개, 상기 값은 예를 들어 특허문헌 1에 개시된 장치를 채용한 실험을 통하여 결정된다.

이러한 실험을 통해, 유체가 유동하는 공간의 벽면에서의 전단 응력 τw와, 벽면에서의 슬립 속도 "v_slip", 그리고 벽면의 이동 속도 "v_wall"을 계측한다.

그 후에, 이상으로부터 "F_slip"의 값을 결정한다. 보다 구체적으로, X축에 슬립 속도 "v_slip"을 플롯하고 전단 응력 τw를 Y축에 플롯한 이중 로그 도표를 준비한 후, 플롯한 점에 대한 누승(累乘) 근사 곡선을

y = a·x^b로 구한다.

"F_slip" 및 "e_slip"은 각각 계수 "a" 및 "b"에 의해 결정된다.

그 후에, 하기의 식 (3)에 의해 벽면에서의 슬립 속도(v_slip)를 구한다:

v_slip = v_wall + τw/F_slip (3)

이렇게 구한 슬립 속도(v_slip)를 솔버측에 제공하고, 이에 따라 솔버는 수렴 계산을 행하여 해를 구한다.

상수 "F_slip"의 값이 작거나 또는 전단 응력 τw의 값이 비정상적이면, 식 (3)에서 구하여 솔버에 제공한 슬립 속도 "v_slip"의 값은 비정상적이게 된다.

따라서, 이러한 슬립 속도가 종래의 유체 시뮬레이션 방법에 도입된다면, 계산이 불확실해질 가능성이 있어, 반복 계산이 수렴하지 않거나 발산한다.

따라서, 본 발명의 목적은, 고점성 유체의 유동 계산을 안정시켜, 고점성 유체의 유동 상태를 정확하게 시뮬레이션하는 전산화된 방법을 제공하는 것이다.

본 발명에 따르면, 벽을 구비한 챔버 내의 고점성 유체를 시뮬레이션하는 전산화된 시뮬레이션 방법은,

상기 챔버의 챔버 모델을 정하는 단계와,

상기 유체의 유체 재료 모델을 정하는 단계, 그리고

상기 유체 재료 모델을 상기 챔버 모델 내에 세팅하고 미리 정해놓은 조건 하에서 유동 계산을 행하는 단계

를 포함하고, 상기 유동 계산에서는, 상기 유체 재료 모델이 접촉하는 상기 챔버 모델의 벽의 접촉면에 대하여, 이 접촉면과 평행한 방향의 상기 유체 재료 모델의 속도인 슬립 속도 "v_slip"을 하기 식 (a)로 정하며,

v_slip = α·v_t+(1-α)·v_wall (a)

여기서

"v_t"는 상기 접촉면으로부터 그 법선 방향으로 거리 "d_wall"를 두고 이격된 위치에 있어서 상기 접촉면과 평행한 방향의 상기 유체 재료 모델의 속도 성분이고,

"v_wall"은 상기 접촉면과 평행한 방향의 상기 접촉면의 속도 성분이며,

"α"는 0 내지 1의 범위의 변수인 슬립율이고,

여기서

"α"는 하기의 식 (b) 또는 (c)를 충족시키며,

α/(1-α) = μ/(d_wall·F_slip) (b)

α/(1-α) = μ/{(d_wall·F_slip) |v_slip-v_wall|^{eslip -1}} (c)

여기서

"μ"는 상기 유체 재료 모델의 점도이고,

"Fslip"은 상수이며, 그리고

"eslip"은 상수인 것이다.

예를 들어, 상기 챔버는 밴버리 혼합기의 케이싱과 이 케이싱 내에 회전 가능하게 배치된 적어도 하나의 로터 사이에 형성된 밴버리 혼합기의 혼련 공간이고, 상기 유체는 로터에 의해 혼련되는 비가교 고무 또는 수지 재료이다.

따라서, 상기 유동 계산에서는 벽에서의 유체의 미끄러짐이 계산된다. 따라서, 비가교 고무 또는 수지 등의 고점성 유체의 유동을 정확하게 시뮬레이션할 수 있다.

식 (a)에서 슬립율 "α"가 0인 경우, 슬립 속도 "v_slip"은 접촉면의 속도 "v_wall"과 동일해진다. 속도 "v_wall"이 0인 경우, 유체 재료 모델의 속도는 0이 된다. 이는, 유체 재료 모델의 미끄러짐이 접촉면에서 일어나지 않는 상태에 해당된다.

식 (a)에서 슬립율 "α"가 1인 경우, 슬립 속도 "v_slip"은, 상기 접촉면으로부터 그 법선 방향으로 거리 "d_wall"를 두고 이격된 위치에 있어서 상기 접촉면과 평행한 방향의 상기 유체 재료 모델의 속도 "v_t"와 동일해진다. 이는, 유체 재료 모델이 접촉면에서 완전히 미끄러지는 상태, 즉 유체 재료 모델과 접촉면 사이에 마찰이 없는 상태에 해당된다.

본 발명에 따르면, 슬립 속도 "v_slip"은 안정적으로 계산될 수 있어, 시뮬레이션에서 행해지는 유동 계산을 확실하게 수렴시킬 수 있다.

도 1은 가소성 재료를 혼련하는 밴버리 혼합기의 주요 구성요소의 단면도이다.
도 2는 본 발명의 실시형태인 유체 시뮬레이션 방법의 흐름도이다.
도 3은 렌더링된 챔버 모델을 사시도로서 보여주는 도면이다.
도 4는 챔버 모델의 단면도이다.
도 5는 기능 부품으로 분할된 챔버 모델의 단면도이다.
도 6은 유체 재료 모델 및 기상 모델이 배치되어 있는 챔버 모델의 상태를 보여주는 개략적인 단면도이다.
도 7은 유동 계산의 흐름도이다.
도 8은 벽면 근방에서 유동하는 유체 재료 모델의 속도를 설명하는 도면이다.
도 9는 슬립율의 값과 식 (b) 또는 (c)의 우변의 값 사이의 관계를 보여주는 그래프이다.
도 10은 본 발명에 따른 시뮬레이션 방법에 의해 구한 가소성 재료의 혼련 상태를 보여주는 도면이다.

이제 본 발명의 실시형태를 첨부 도면과 함께 상세히 설명한다.

본 발명은 챔버 내에서 유동하는 유체의 유동 상태를 컴퓨터(도시 생략)를 사용하여 시뮬레이션하는 방법에 관한 것이다.

여기서, 대상 유체는 챔버의 벽의 접촉면에서 비교적 큰 슬립이 일어나는 고점성 유체이다. 가소성 재료가 충분히 혼련되고 완전히 혼합되어 안정적인 가소성 재료의 유동이 일어나는 상태로 된다는 전제 하에서, 예컨대 비가교 고무, 수지 및 엘라스토머 등의 가소성 재료가 상기 고점성 유체에 해당한다. 특히, 온도가 약 80℃까지 상승한 고무는 이러한 상태인 것으로 고려될 수 있다.

챔버는 벽면에 의해 둘러싸여 있고 고점성 유체가 이동하는 공간이다. 챔버는 다양한 형상 또는 구성을 가질 수 있다. 대표적인 챔버는, 재료가 혼련 및 혼합되는 밴버리 혼합기의 혼련 공간이다.

도 1은 밴버리 혼합기(1)의 개략적인 단면도이다.

본 예에서의 밴버리 혼합기(1)는, 케이싱(2), 케이싱(2) 내에 회전 가능하게 배치된 한 쌍의 로터(3), 및 케이싱(2)과 두 로터(3)의 사이에 형성되며 8자 형상을 지닌 혼련 공간으로서의 챔버(4)를 포함한다. 물론 챔버(4)는 8자 형상 이외의 다양한 형상을 가질 수 있다.

도 2는 본 발명의 실시형태인 시뮬레이션 방법의 흐름도를 보여준다.

* 챔버 모델을 정하는 단계 S1

본 실시형태의 시뮬레이션 방법에 따라, 우선 챔버(4)의 유한 요소 모델인 챔버 모델을 정하고 컴퓨터에 저장한다.

도 3은 렌더링된 챔버 모델(5)을 사시도로서 보여주는 도면이다. 도 4는 이 챔버 모델의 단면도이다.

챔버 모델(5)은 유한개의 3차원 요소로 분할된 3차원 공간으로서, 케이싱(2)의 내주면에 의해 형성된 외주면(5o)과, 회전 가능한 각 로터(3)의 외주면에 의해 형성된 내주면(5i), 그리고 로터(3)의 축방향 양단측에 있는 단부면(5s)으로 둘러싸인 공간이다.

이들 면(5o, 5i 및 5s)은 전술한 챔버 모델(5)의 벽면이다.

전술한 챔버 모델(5)의 3차원 요소는 오일러(Euler) 요소이다. 예를 들어, 4면체 요소, 6면체 요소 등의 한 종류, 또는 복수 종류의 다면체 요소가 사용된다.

이후에 설명하는 바와 같이, 각 요소에 대하여, 예컨대 압력, 온도, 속도 등의 가소성 재료(유체 재료 모델)의 물리량을 계산한다.

챔버 모델(5)의 외주면(5o)과 단부면(5s)은 변형될 수 없는 것으로 정해진다.

챔버 모델(5)의 내주면(5i)은 변형될 수 없고 로터(3)의 회전에 따라 회전 가능한 것으로 정해져, 챔버 모델(5)의 형상의 변화가 허용된다.

본 실시형태에서는, 챔버 모델(5)의 형상 변화를 효과적으로 시뮬레이션하기 위해, 챔버 모델(5)을 도 5에 도시된 바와 같이 한 쌍의 회전부(5A, 5B), 이들 회전부 사이에 있는 연결부(5C) 및 전술한 부분(5A, 5B, 5C)을 에워싸는 외측 프레임부(5D)의 4개의 부분으로 나눈다.

회전부(5A/5B)는 각각 원형 외주면(5Ao/5Bo)과 하나의 로터(3)의 외주면에 대응하는 내주면(5i)을 갖는 관 모양의 것이다.

회전부(5A, 5B)는 외측 프레임부(5D) 내에 놓이고, 로터(3)의 회전에 의해 야기되는 챔버(4)의 용적 형상의 변화를 표현하기 위해 각 중심축(Oa, Ob)을 둘레로 회전 가능한 것으로 정해진다.

이에 반하여, 중간 연결부(5C)는 회전부(5A, 5B)의 사이에 정지된 상태로 있고, 각 회전부(5A, 5B)에 접하는 2개의 오목면(j)을 갖는다.

오목면(j)과 각 원통형 외주면(5Ao, 5Bo)에서는, 경계 조건이 미끄럼면으로서 정해진다. 이로써, 회전부(5A, 5B)에서 일어나는 물리적 작용(힘, 열 등)이, 상기 중간 연결부(5C) 내에 있는 유체 모델에 오목면(j)을 통하여 전달될 수 있게 된다.

외측 프레임부(5D)는 관 모양의 것이고 상기 회전부(5A, 5B)와 상기 중간 연결부(5C)를 둘러싼다. 외측 프레임부의 축방향 양단부는 2개의 단부면(5s)에 의해 폐쇄되어 있다.

외측 프레임부(5D)와 회전부(5A, 5B) 사이의 계면 및 외측 프레임부(5D)와 중간 연결부(5C) 사이의 계면에서는, 경계 조건이 미끄럼면으로서 정해진다. 이로써, 회전부(5A, 5B)에서 일어나는 물리적 작용(힘, 열 등)이, 상기 두 계면을 통하여 외측 프레임부(5D)에 전달될 수 있게 된다.

로터의 작용에 의해, 외측 프레임부(5D)는 비교적 큰 전단력을 받게 된다. 따라서, 재료의 속도 등을 보다 상세히 계산하기 위해서는, 외측 프레임부(5D)를 구성하는 요소를 회전부(5A, 5B) 및 중간 연결부(5C)보다 작은 크기로 만드는 것이 바람직하다. 이로써, 챔버 모델(5)의 내면 근방에서의 유체 재료 모델의 속도 프로파일 등을 보다 상세히 계산할 수 있다.

* 유체 모델을 정하는 단계 S2

이어서, 컴퓨터를 이용하여 유체 재료 모델을 정하고 저장한다.

유체 재료 모델은 챔버(4) 내에서 유동 또는 이동하는 유체(가소성 재료)의 모델이다. 유체 재료 모델에 있어서, 전단 점도, 비열, 열전도율 및 비중 등의 가소성 재료의 물리적 물성치를 정하고 컴퓨터에 저장한다.

전단 점도에 관해서는, 분석 대상(가소성 재료)을 복수의 온도 조건 하에서 점탄성 특성(G', G")에 대해 측정하고, 이 점탄성 특성을 Cox-Merz 법칙에 따라 변환함으로써, 전단 점도를 구한다.

전술한 바와 같이 구한 전단 점도 μ를, 누승 법칙에 따른 하기의 식 (f)으로 근사한다.

μ = mγ'^n-1 (4)

여기서

m: 절대 온도 T의 함수인 계수이고,

γ': 전단 속도이며,

n: 계수이다.

유체 재료 모델의 비열에 관하여는, 분석 대상인 가소성 재료를 단열형 연속 가열법(@ 25℃)으로 측정하고, 측정한 비열값을 미리 컴퓨터에 입력 및 저장한다.

유체 재료 모델의 열전도율에 관하여는, 분석 대상인 가소성 재료를 열선법(@ 25℃)으로 측정하고, 측정한 열전도율값을 미리 컴퓨터에 입력 및 저장한다.

* 기상 모델을 정하는 단계 S3

이어서, 컴퓨터를 이용하여 기상 모델을 정하고 저장한다.

기상 모델은 챔버 내에 제2 유체로서 존재하는 기체의 모델이다.

본 실시형태에서, 가소성 재료의 챔버에의 충전율은 100% 미만이다. 따라서, 유동 계산을 가능하게 하기 위해, 챔버 모델(5)의 일부에는 유체 재료 모델을 채우지 않고, 기상 모델을 채운다.

기상 모델에 있어서, 기체의 점도 및 비중을 정하고 컴퓨터에 저장한다.

대개, 특정 온도에서의 공기의 점도 및 비중의 실제값이 기상 모델에 세팅된다. 그러나, 실제값과는 다른 값이 세팅될 수도 있다.

점도가 서로 크게 다른 기상 모델(공기)과 유체 재료 모델(가소성 재료)의 혼상 유동 분석을 실시하는 경우, 기상 모델과 유체 재료 모델의 계면에서 전단열 발생이 증가되어, 유동 계산이 불안정해질 가능성이 있다.

본 실시형태에서는, 안정적인 유동 계산을 위하여, 기상 모델에 세팅된 점도의 값을, 계산 결과에 악영향을 미치지 않는 한도 내에서, 가능한 실제값에서 증가시킨다.

본원의 발명자들은, 챔버 모델(5) 내에서의 압력장을 비교하기 위해, 기상 모델에 세팅된 점도의 값만을 변화시킨 다양한 조건 하에서 유동 계산(유동 시뮬레이션)을 실시하였다.

그 결과, 기상 모델의 점도의 값이 공기의 점도의 실제값의 10배를 초과한다면, 챔버 모델(5) 내에서의 압력의 값이 증가되고 압력장에 악영향을 미치는 것으로 확인되었다.

한편, 기상 모델의 점도의 값이 공기의 점도의 실제값의 5배 미만이면, 유동 계산이 불안정해지는 것을 완전히 방지하기 어렵다.

따라서, 공기의 점도의 실제값에 대해 5배 내지 10배의 범위 내의 값이 기상 모델에 세팅되는 것이 바람직하다.

* 경계 조건을 정하는 단계 S4

이어서, 유동 계산을 실시하는 데 필요한 경계 조건 등의 여러 조건을 정한다.

경계 조건으로는, 챔버 모델(5)의 벽면에서의 유속 경계 조건 및 온도 경계 조건 등이 있다.

유속 경계 조건에 관해 살펴보면, 유체 재료 모델이 접촉하는 챔버 모델(5)의 벽면의 접촉부에는, 유동 계산 과정에서 유체 재료 모델이 상기 접촉부에서 소정의 속도 또는 슬립 속도를 가질 수 있도록, 벽면 슬립 조건이 제공된다.

온도 경계 조건에 관해 살펴보면, (a) 열이 챔버 모델(5)의 표면을 통하여 외부로 빠져나가지 않는 단열 조건, 또는 (b) 챔버 모델(5)의 전체 표면이 일정한 온도(예컨대, 50℃)를 갖는 조건이 시뮬레이션의 목적, 요구되는 정확도 등에 따라 정해질 수 있다.

또한, 상기 조건에는 유체 재료 모델의 초기 온도가 포함될 수 있다. 본 실시형태에서, 초기 온도는 20℃로 세팅되어 있다.

또한, 상기 조건에는 밴버리 혼합기의 로터의 회전수에 대응하는 챔버 모델의 회전부(5A, 5B)의 회전수가 포함될 수 있다.

더 나아가, 상기 조건에는 챔버 모델(5)의 전체 용적에 대한 유체 재료 모델의 충전율이 포함될 수 있다.

또한, 유동 계산의 초기 상태, 계산의 시간 간격, 내부 처리에서의 반복 횟수, 계산(반복)의 최대 기간 등이 조건으로서 포함될 수 있다.

초기 상태에 대해서는, 예를 들어 도 6에 도시된 바와 같이, 챔버 모델(5)을 가로질러 연장되는 것으로 정해진 수평 경계면(S)의 상측의 영역 A를 기상 모델이라 하고, 상기 수평 경계면(S)의 하측의 영역 M을 유체 재료 모델이라 하는 것으로 정할 수 있다.

따라서, 상기 경계면(S)의 레벨을 변경함으로써, 유체 재료 모델(가소성 재료)의 충전율이 조정될 수 있다.

*유동 계산을 실행하는 단계 S5

본 실시형태에서는, 도 6에 도시된 바와 같이, 유체 재료 모델(영역 M)과 기상 모델(영역 A)이 챔버 모델(5) 내에 세팅되어 있고, 전술한 바와 같이 정해진 조건에 따라 유동 계산이 행해진다.

유동 계산에서는, 적어도 다음 5가지 미지량, 즉 유체 재료 모델의 3개의 좌표축 방향(x, y, z 방향)의 3가지 속도 성분과 유체 재료 모델의 압력(p) 및 온도(T)가 계산된다.

본 실시형태에서는, 비압축성 유동의 경우에 네비어-스톡스 방정식에 기초하여 유동 계산을 실시한다.

따라서, 유동 계산의 과정에서, 기상 모델 및 유체 재료 모델의 밀도는 일정한 것으로 취급된다.

본 실시형태에서는, 유동 계산의 과정에서, 유체 재료 모델은 전체 온도 범위에 걸쳐 유체인 것으로 취급된다. 따라서, 풀어야 하는 유체 방정식은 연립 방정식(네비어-스톡스 방정식, 질량 보존 방정식 및 에너지 방정식)이다.

본 실시형태의 유동 계산에서는, 챔버 모델(5) 내에 존재하는 2종류의 유체, 즉 기상 모델 및 유체 재료 모델을 한꺼번에 취급할 필요가 있다.

이러한 목적으로, 본 실시형태에서는, 자유 계면을 갖는 유동을 계산하는 데 사용되는 VOF(Volume of Fluid)법이 채용된다.

VOF법은 2종류의 유체간의 계면의 움직임을 직접 계산하지 않는다. VOF법에서는, 챔버 모델(5)의 각 요소의 체적 내에 있어서 유체 재료 모델의 충전율에 상당하는 체적분율을 정함으로써 자유 계면을 표현한다.

지배 방정식은 다음과 같다.

* [운동 방정식]

본 실시형태에서는, 기상 모델과 유체 재료 모델이 챔버 모델(5) 내에서 함께 유동하는 2상 유동을 단상 유동으로 취급한다. 이러한 경우에, 3개의 좌표축 방향(x, y, z)에서 풀어야 하는 운동 방정식은 하기의 식 (5)이다. 이는, 2개의 상을 VOF법에 의해 평균화하고 단일 상으로 취급한 결과로서 가능해진다.

식 (5):

여기서

u: 혼상 유체 모델의 속도,

p: 혼상 유체 모델의 압력,

ρ: 혼상 유체 모델의 밀도,

g: 중력 가속도,

T: 혼상 유체 모델의 절대 온도,

F: 외력.

유체 재료 모델 및 기상 모델을 포함한 요소의 밀도 ρ 및 점성 계수 μ는, 하기의 식 (6)에 나타내어진 바와 같이, 각 상(즉, 유체 재료 모델 및 기상 모델)이 차지하는 체적에 의해 가중치가 부여된 후, 평균화된다.

식 (6):

여기서

δq: 체적분율

ρq: 각 요소에 있어서 각 상의 밀도

μq: 각 요소에 있어서 각 상의 점도

*[질량 보존 방정식]

질량 보존 방정식(연속 방정식) 및 압력 방정식에 관해서는, 단 하나의 세트의 3개의 좌표축 방향의 방정식만을 풀면 충분하다. 따라서, 본 실시형태의 시뮬레이션 방법에 따르면, 컴퓨터는 유동장을 혼상임에도 불구하고 단상으로서 계산할 수 있다. 다시 말하자면, 위치(체적분율)에 따라 재료 물성치가 달라지는 유동을 푼다.

상기 계산의 결과로서 구해지는 체적분율의 분포에 따라 각 상의 위치를 판단할 수 있다.

* [에너지 방정식]

유체 재료 모델의 온도는 하기의 식 (6)에 의해 구할 수 있다.

식 (7):

여기서

E: 엔탈피,

k: 열전도율,

S: 소스 항.

* [체적분율의 수송 방정식]

체적분율의 분포는 두 상(또는 기상 모델과 유체 재료 모델) 사이의 계면의 위치를 결정한다.

이 체적분율 δq는 하기의 식 (8)을 정확히 푸는 것에 의해 구할 수 있다.

식 (8):

챔버 모델(5)의 임의의 요소(e)에 있어서 체적분율 δq가 0이면, 이는 q상 (q상=유체 재료 모델)이 이 요소(e)에 존재하지 않는 것을 의미한다.

체적분율 δq가 1이면, 이는 요소(e)의 전체 용적에 q상이 채워져 있는 것을 의미한다.

0<δq< 1이면, 이는 요소(e)에 q상(유체 재료 모델)과 다른 상(기상 모델)이 채워져 있고, 즉 요소는 계면을 갖는 혼상임을 의미한다.

이 방정식은 "ANSYS Fluent User's Manual, 26.2.9 Modified HRIC Scheme"에 상세하게 설명되어 있는 Modified-HRIC (Implicit)에 의해 풀릴 수 있다.

본 실시형태에서, 전술한 방정식 각각은 압력 베이스의 분리법에 의해 풀린다.

압력 방정식과 운동 방정식의 결합에는, SIMPLE(Semi-Implicit Method for Pressure-Linked Equations) 알고리즘이 사용되는 것이 바람직하다.

도 7은 컴퓨터에 의해 행해진 유동 계산의 일례의 흐름도이다.

** 스텝 S51

본 예에서는, 우선 압력 구배 및 속도 구배에 대한 상한 및 하한을 세팅하고, 속도 구배 및 압력 구배를 계산한다.

** 스텝 S52

이어서, 현재 압력장의 이산화를 통하여 운동 방정식을 정하고, 반복법에 의해 상기 운동 방정식의 해를 구한다. 즉, 재료 모델(또는 혼상)의 3개의 좌표축 방향의 속도를 계산한다. 반복법으로서는, 가우스-자이델법을 사용할 수 있다.

** 스텝 S53

이어서, 전술한 속도가 질량 보존 방정식을 충족시키는가의 여부를 검사하기 위해, 우선 챔버 모델(5)의 요소의 표면에 있어서 유체 재료 모델의 수정전 질량 유량을 계산한다.

여기서, "수정전 질량 유량"은 SIMPLE 알고리즘의 루프의 초반에 일시적으로 사용되는 질량 유량이다.

이러한 질량 유량은 큰 오차를 가질 수 있으므로, "수정전 질량 유량"이라 한다.

요소의 질량 유량은, 해당 요소의 전체 표면을 통과하는 해당 물질의 질량 유량이다.

** 스텝 S54

이어서, SIMPLE 알고리즘을 이용하여, 속도장과 압력장을 결합하고, 압력장을 수정하기 위한 하기의 압력 보정 방정식 (9)를 생성한다.

식 (9):

** 스텝 S55

이어서, 예컨대 AMG 솔버, CG, 또는 Bi-CG 등의 반복법에 의해, 압력 보정 방정식의 해를 구하고, 압력 보정량 p'를 계산한다.

** 스텝 S56

이어서, 이렇게 구한 해에 기초하여, 압력장을 하기 식 (10)에 의해 보정한다:

p^{n +1} = pⁿ + ωp' (10)

여기서

p: 압력,

n: 현재의 시간 스텝의 수,

ω: 완화 계수.

본 예에서, 완화 계수는 0.3으로 세팅되지만, 완화 계수에 다른 값이 세팅될 수도 있다.

** 스텝 S57

이어서, 계면의 경계 조건을 수정(업데이트)한다. 여기서, 계면은, 유체 재료 모델 및 기상 모델을 포함하는 유체 모델과, 그 외의 고체 모델(로터 및 챔버) 사이의 계면이다. 구체적으로, 보정된 압력장으로부터 압력 구배를 구한다. 이렇게 구한 압력 구배는 경계 조건으로서 제공된다.

** 스텝 S58

이어서, 하기의 식 (11)에 의해 요소의 표면에서의 질량 유량을 구한다:

m_f ^{n +1} = m_{f *} + m'_f (11)

여기서

m_f ^{n +1}: 수정전의 요소의 표면에서의 질량 유량,

m_{f *}: 수정후의 요소의 표면에서의 질량 유량,

m'_f: 질량 유량의 수정값.

** 스텝 S59

이어서, 하기의 식 (12)에 의해 속도장을 수정한다:

v^{n +1} = v*-(V▽p'/∂_pV) ---(12)

여기서

V: 요소의 체적,

v*: 운동 방정식에서 구한, 수정 이전의 중간 속도장,

∂_pV: 운동 방정식의 매트릭스의 대각 성분,

▽p': 압력 보정량의 구배.

** 스텝 S60

이어서, 전술한 에너지 방정식 (7)의 해를 구함으로써, 재료 모델의 온도 및 점도를 계산한다.

** 스텝 S61

또한, 챔버 모델(5)의 벽면의 슬립 속도 "v_slip"을 식 (a)에 의해 계산한다:

v_slip = α·v_t+(1-α)·v_wall

도 8에 도시된 바와 같이, 식 (a)에 있어서,

"v_t"는 챔버 모델의 접촉면(CP)으로부터 그 법선 방향으로 거리 "d_wall"를 두고 이격된 위치에 있어서 상기 접촉면(CP)과 평행한 방향의 상기 유체 재료 모델의 속도 성분이고,

"v_wall"은 접촉면(CP)의 속도에 있어서 접촉면(CP)과 평행한 방향의 성분이며,

"α"는 0 내지 1의 범위의 변수인 슬립율이다.

슬립 속도가 도 8에 도시된 바와 같이 선형인 경우, 슬립율 "α"는 하기의 식 (b)를 충족시킨다.

슬립 속도가 비선형인 경우(도시 생략), 슬립율 "α"는 하기의 식 (c)를 충족시킨다.

식 (b):

α/(1-α) = μ/(d_wall·F_slip)

식 (c):

α/(1-α) = μ/{(d_wall·F_slip) |v_slip-v_wall|^{eslip -1}} (c)

식 (b) 및 (c)에서,

"μ"는 유체 재료 모델의 점도이고,

"F_slip"은 상수이며,

"e_slip"은 상수이다.

"Computational Methods For Fluid Dynamics", P256-P257, 식 (8.74) 및 (8.75) (저자: Joel H. Ferziger 및 Milovan Peric, 출판사: Springer)에 따르면, 슬립 속도가 선형인 경우, 접촉면(CP)에서의 전단 응력 τw은 일반적으로 하기의 식(13)에 의해 주어진다:

식 (13):

식 (1)과 식 (13)으로부터, 하기의 식 (14)가 얻어진다.

식 (14):

식 (14)에 식 (a)를 대입함으로써, 하기의 식 (15)가 얻어진다.

식 (15):

식 (15)는 식 (b)로 표현될 수 있다.

식 (15) 또는 (b)에서, "μ"는 상기 유체 재료 모델의 점도로서, 고유의 값이고, "dwall"은 알려진 값이다. 예를 들어, 접촉면(CP)에 접하는 요소의 중심(대표점)에 대한 접촉면(CP)으로부터의 거리가 "dwall"로 세팅될 수 있다. 따라서, 이러한 경우에 "dwall"의 값은 챔버 모델(5)의 메싱(meshing)에 의해 결정된다.

식 (b)에서, "Fslip"은 유저에 의해 정해지는 기지의 상수이다.

따라서, 슬립율 "α"는 식 (b)로부터 쉽게 결정될 수 있다.

슬립 속도 "v_slip"은, 슬립율 "α"를 식 (a)에 대입함으로써 쉽게 구해질 수 있다.

이러한 식들을 이용함으로써, 유체 재료 모델의 슬립 속도 "v_slip"이 비정상적인 값을 나타내는 것을 피할 수 있다. 따라서, 유체 시뮬레이션의 로버스트성이 향상되고, 급증한 속도에 기인하는 정확도의 악화가 방지될 수 있다. 또한, 계산의 수렴성이 향상되어, 계산 시간을 단축할 수 있다.

슬립 속도가 비선형인 경우, 식 (2)와 식 (13)으로부터, 하기의 식 (16)이 얻어진다.

식 (16):

식 (16)과 식 (a)으로부터, 전술한 식 (c)가 얻어진다.

슬립 속도 "v_slip"과 접촉면에서의 속도 "v_wall" 사이의 차가 매우 작은 경우, 식 (c)의 분모의 값이 매우 작아진다. 이러한 경우, 식 (c)의 계산 과정에서, 소위 "0으로 나눔"이 야기될 가능성이 있다.

따라서, 비선형 슬립 속도가 정해진 경우에는, 비선형 슬립 속도에 기초하여 안정적인 계산을 실시하기 위해, 식 (c)에서의 속도차(v_slip-v_wall)가 미리 정해놓은 한계값(예컨대, 약 1.0×10^-8)보다 작아질 때, 식 (b)를 사용하는 것이 바람직하다.

도 9는 식 (b) 또는 (c)의 우변의 값과 슬립율 "α"의 값 사이의 관계를 보여주는 그래프이다.

항상 양의 값인 우변의 값이 크게 변하더라도, 슬립율 "α"의 값의 변경 범위는 0과 1의 사이이다.

따라서, 본 발명에 따르면, 식 (b)를 사용하는 한, 슬립율 "α"은 결코 비정상적인 값을 갖지 않는다.

** 스텝 S62

이어서, 상기 계산의 해가 수렴하는지를 판단한다. 여기서, 상기 해는, 본 실시형태의 유체 시뮬레이션의 지배 방정식인 네비어-스톡스 방정식[식 (5)+연속 방정식]의 해이다.

수정된 질량 유량의 총계가 미리 정해놓은 오차 범위 내에 있는가의 여부에 기초하여, 수렴을 판정할 수 있다.

수렴하지 않는 것으로 판단된다면(스텝 S62에서 N), 컴퓨터는 스텝 S51과 후속 스텝을 다시 수행한다.

** 스텝 S63

수렴한 것으로 판단된다면(스텝 S62에서 Y), 시간 스텝을 1 증분하여, 시각을 하나의 시간 간격만큼 앞당긴다.

** 스텝 S64

또한, 현재의 시간 스텝이 미리 정해놓은 시간 스텝에 도달하지 않았는가의 여부를 판단한다.

현재의 시간 스텝이 미리 정해놓은 시간 스텝에 도달하지 않았다면(스텝 S64에서 N), 컴퓨터는 스텝 S51과 후속 스텝을 다시 수행한다.

본 실시형태의 유동 계산에서, 유체 재료 모델의 전단 가열은 계산되지만, 기상 모델의 전단 가열은 계산되지 않는다.

가소성 재료이기 때문에, 유체 재료 모델의 전단 발열은 고려되어야만 한다.

유체 재료 모델과 기상 모델 사이의 계면에서, 전단 속도는 대개 증대된다.

전단 발열은 전단 점도에 전단 속도의 제곱을 곱한 것에 상당하므로, 전단 속도의 값이 매우 큰 경우, 온도의 계산(에너지 방정식)에 부정적인 영향을 미칠 가능성이 있다.

그래서, 본 실시형태에서는, 전단 발열의 계산을, 유체 재료 모델에만, 보다 구체적으로 체적분율 δq가 일정값(본 실시형태에서는 0.90) 이상인 요소에만 한정하여, 계산의 안정화를 도모한다.

비교 시험

본 발명의 유익한 효과를 확인하기 위해, 도 3에 도시된 챔버 모델과 미가황 고무의 유체 재료 모델을 사용하여, 전술한 방법에 따라, 실시간으로 20초간 지속된 혼련의 시뮬레이션을 행하였는데, 이 시뮬레이션에서 미가황 고무의 챔버 모델에의 충전율은 70%로 하였고, 로터의 회전수는 30 rpm으로 하였으며, 계산 시간 간격은 1.973×10^-3로 하였다.

식 (a) 및 (b)를 사용하여 슬립 속도를 계산하는 과정에서, 비정상적인 슬립 속도의 값은 발생하지 않았으며, 유동 계산은 안정적으로 수렴하였다.

이 시뮬레이션의 결과가 도 10에 도시되어 있는데, 도면에서 혼련된 고무의 상태는 시간순으로 도시되어 있고, 흑색 부분은 유체 재료 모델이며, 회색 부분은 기상 모델이다.

한편, 본원의 종래 기술 섹션에서 설명한 방법에 따르면, 슬립 속도의 계산 과정에서 비정상적인 값이 빈번히 발생되었으며, 유동 계산은 수렴하지 않았다.

Claims (2)

1. 벽을 구비한 챔버 내의 고점성 유체를 시뮬레이션하는 전산화된 시뮬레이션 방법으로서,
   상기 챔버의 챔버 모델을 정하는 단계와,
   상기 유체의 유체 재료 모델을 정하는 단계, 그리고
   상기 유체 재료 모델을 상기 챔버 모델 내에 세팅하고 미리 정해놓은 조건 하에서 유동 계산을 행하는 단계
   를 포함하고,
   상기 유동 계산에서는, 상기 유체 재료 모델이 접촉하는 상기 챔버 모델의 벽의 접촉면에 대하여, 이 접촉면과 평행한 방향의 상기 유체 재료 모델의 속도인 슬립 속도 "v_slip"을 하기 식 (a)로 정하며,
   v_slip = α·v_t+(1-α)·v_wall (a)
   여기서
   "v_t"는 상기 접촉면으로부터 그 법선 방향으로 거리 "d_wall"를 두고 이격된 위치에 있어서 상기 접촉면과 평행한 방향의 상기 유체 재료 모델의 속도 성분이고,
   "v_wall"은 상기 접촉면과 평행한 방향의 상기 접촉면의 속도 성분이며,
   "α"는 0 내지 1의 범위의 변수인 슬립율이고,
   여기서
   "α"는 하기의 식 (b) 또는 (c)를 충족시키며,
   α/(1-α) = μ/(d_wall·F_slip) (b)
   α/(1-α) = μ/{(d_wall·F_slip) |v_slip-v_wall|^{eslip -1}} (c)
   여기서
   "μ"는 상기 유체 재료 모델의 점도이고,
   "F_slip"은 상수이며,
   "e_slip"은 상수인 것인 시뮬레이션 방법.
2. 제1항에 있어서, 상기 챔버는 밴버리 혼합기의 케이싱과 이 케이싱 내에 회전 가능하게 배치된 적어도 하나의 로터 사이에 형성된 밴버리 혼합기의 혼련 공간이고,
   상기 유체는 상기 적어도 하나의 로터에 의해 혼련되는 비가교 고무 또는 수지 재료인 것인 시뮬레이션 방법.

Images