[FLOW-3D 이론] Auxiliary Model/Elastic Membrane and Elastic Wall Model탄성막과 탄성벽 모델

제한된 유체구조 상호작용(FSI) 기능이 탄성막과 벽모델을 이용하여 FLOW-3D 에서 가능하다. 이 모델에서 탄성막이나 탄성벽의 변형은 인접 유체유동에 영향을 미치고 교대로 유체 압력은 변형에 작용한다. 이런 상호작용이 완전히 결합된 방식으로 FLOW-3D 내에 기술되어 있다.

이 모델의 주요한 한계는 변형이 작다, 즉 각 막이나 탄 성벽은 그 휨이 크기에 비해 아주 작다고 가정된 것이다. 이는 대신에 모델을 유용하게 단순화를 가능케 한다. 박막과 탄성벽의 형상은 계산시 시간에 따라 변하지 않는다고 가정되고 한편 유체유동에 대한 이 변형의 효과는 유체구조 경계면 상에 분포된 체적 소스나 싱크로 기술된다. 압력이 막표면에 균일하게 작용한다는 추가 가정 하에 더 나은 계산 효율을 위해 구조해석 알고리즘보다 오히려 수치해가 박막의 변형을 결정하는데 이용된다.

Elastic Membrane 타성막

FLOW-3D 에서 탄성 막은 외력의 작용하에 작은 탄성 변형을 겪는 직각 또는 원형의 얇은 판이다. 이의 두께와 재질은 균일하다고 가정된다. 이의 모서리들은 단순히 지지되어 있거나 고정되 있을 수 있다. 단순 지지 모서리는 휨과 순수모멘트가 0인 모서리이다. 그러나 고정 모서리는 휨과 1차 미분값은 0이나 힘 모멘트는 일반적으로 0이 아니다. 어느 경우던지 이 모델은 막이 모서리를 따라서 다 같은 조건을 가지는것을 필요로 한다. 망 격자 내에 막의 위치는 제한이 없으나 막의 표면은 x, y, 또는 z 축에 수직이어야 한다.

이 모델은 막에 작용하는 두 외력을 고려한다: 수리력과 작동기 힘. 수리력은 막의 양편에 작용하는 압력을 적분하여 얻어진다. 그런 후에 전체 막 위에 균일하게 분포된 힘으로 전환된다.

작동기 힘은 마이크로-펌프 유동, 잉크젯 방울 형성같은 많은 응용에 존재한다. 한 예는 막에 부착된 압전 작동기이다. 전압이 가해지면 압전작동기는 막표면에 수직한 방향으로 힘을 가하는데 이것이 소위 작동기 힘이라고 불린다. 사용자는 작동기 힘을 사인파나 구간선형 인 시간의 함수로 지정할 수 있다. 이 모델에서 작동기는 항상 막의 한편 중심에 위치한다고 가정되고 작동기의 막과 접촉하는 면적의 형상(반드시 크기일 필요는 없지만)은 막의 형상과 같다. 다른 말로 막과 작동기는 직각 또는 원형이어야하고 같은 대칭축을 가져야한다. 작동기 자체가 보통 막보다 낮은 강직성을 가지므로 작동기의 힘은 항상 접촉 면적상에 균일하게 작용한다고 가정된다. 접촉 면적이 없으면 작동기 힘은 막의 중심에 작용하는 집중된 힘으로 처리된다.

더 나은 계산효과를 위해 더 나은 계산 효율을 위해 구조해석 알고리즘보다 오히려 해석해가 휨의 계산에 이용된다. 임의의 시간과 점에서 막은 수리력, 작동기 힘 그리고 막의 강성의 균형에 의해 정의되는 평형을 이루고 있다고 가정된다. 해석해는 작은 변형을 가지는 박판의 평형방정식을 해석함으로써 얻어진다.

(120){\nabla ^2}{\nabla ^2}w = \frac{f}{D}

여기서

  • w 는 휨
  • f 는 막의 단위 면적당 순수외력
  • D 는 굴곡 강성

(121)D = \frac{E{h^3}}{12\,(1 - {\nu ^2})}

여기서

  • E 는 영의 탄성계수
  • ν 는 포아송 비이며
  • h 는 막두께이다.
  • z 축에 수직인 표면을 가지는 직사각형의 막을 고려한다. 데카르트 좌표에서 시스템 식 (10.120) 은 다음과 같이 표현된다.

(122)\frac{{{\partial ^4}w}}{{\partial {x^4}}} + 2\frac{{{\partial ^4}w}}{{\partial {x^2}\partial {y^2}}} + \frac{{{\partial ^4}w}}{{\partial {y^4}}} = \frac{f}{D}

편리하도록 좌표 원점을 막중심에 잡는다. a b 를 x 와 y 의 각기 방향에서의 막의 길이라고하자. 단순지지 모서리에 대한 막의 경계조건은

(123)\begin{gathered}
  w = 0 \quad \text{    and    } \quad \frac{{\partial ^2}w}{\partial {x^2}} = 0 \quad \text{    for    } \quad x =  \pm \frac{a}{2} \\
  \\
  w = 0 \quad \text{    and    } \quad  \frac{{\partial ^2}w}{\partial {y^2}} = 0 \quad \text{    for    } \quad y =  \pm \frac{b}{2} \\
  \end{gathered}

막이 고정된 모서리를 가진다면 막에대한 경계조건은

(124)\begin{gathered}
  w = 0 \quad \text{    and    } \quad \frac{\partial w}{\partial x} = 0 \quad \text{    for    } \quad x =  \pm \frac{a}{2} \\
  \\
  w = 0 \quad \text{    and    } \quad  \frac{\partial w}{\partial y} = 0 \quad \text{    for    } \quad y =  \pm \frac{b}{2} \\
  \end{gathered}

이다.

막에 대해 식 (10.120) 을 원점이 막의 중심에 위치한 원통좌표계로 쓰는 것이 편리하다,

(125)\frac{1}{r}\frac{d}{{dr}}\left\{ {r\frac{d}{{dr}}\left[ {\frac{1}{r}\frac{d}{{dr}}\left( {r\frac{{dw}}{{dr}}} \right)} \right]} \right\} = \frac{f}{D}

막의 반경이 a 일 때 단순 지지된 모서리에 대한 경계조건은

(126)w = 0 \quad \text{    and    } \quad \frac{{{\partial ^2}w}}{{\partial {r^2}}} = 0 \quad \text{    for    } \quad r = a

이다.

고정 모서리에 대한 경계조건은

(127)w = 0 \quad \text{    and    } \quad \frac{\partial w}{\partial r} = 0 \quad \text{    for    } \quad r = a

이다.

이 모델에서 사용되는 막의 휨에 대한 모든 해석해는 상기 식들을 만족한다. 일부 해석해는 Timoshenko (1959) 에서 찾을 수 있으나 다른 해들은 중첩 방식을 이용하여 Timoshenko 해로부터 얻어진다. 이런 해의 상세한 내용은 Flow Science Technical Note #79 를 참조하라.

유동에 대한 막의 운동효과를 고려하기위해 연속방정식이 우측에추가되는 체적 소스(또는싱크)항 S 로써 수정된다,

(128)\frac{{{V_f}}}{\rho }\frac{{\partial \rho }}{{\partial \,t}} + \frac{1}{\rho }\nabla  \cdot (\rho \,\vec uA) = S

계산 유효체적, 또는 망 셀에서,

(129)S = \frac{S_{\rm{mb}}}{V_{\rm{cell}}}\, {\vec V_{\rm{mb}}} \cdot \vec n

여기서 Vcell 는 셀 체적, Smb, ~n V~mb 는 망 셀에서의 각기 표면적, 단위 외부 법선 벡터와 막 표면의 속도이다. V~mb 는 처짐의 변화율로부터 얻어진다. VOF 함수의 이송방정식 또한 소스항 FS 를 이용하여 수정된다,

(130){V_f}\frac{\partial F}{\partial \,t} + \nabla  \cdot (F\vec u{A_f}) = FS

모멘텀, 에너지, 난류 그리고 스칼라들의 이송방정식은 변하지 않는데 그 이유는 FLOW-3D 에서 이 방정식들은 비보존 형태로 사용되기 때문이다. 이들이 연속방정식을 고려하여 보존형태로 유도될 때 막운동에 의한 소스항들이 상쇄된다.

Elastic Wall 탄성벽

FLOW-3D 에서 탄성벽은 임의 형상의 탄성물체이며, 그 표면 변형은 작고 수압에 비례한다. 즉,

(131)w =  - \frac{{(p - {p_{\rm{ref}})}}}{K}

여기서

  • w 는 표면 법선방향에서의 지역 휨
  • p 는 지역압력
  • pref 는기준압력
  • K 는 단위면적당 강성계수

이러한 종류의 탄성변형은 벽재질의 포아송비가 0이면 발생하며 이는 법선응력이 횡 변형을 일으키지 않는 것을 뜻한다. 포아송비가 0이 아니지만 작으면 탄성 벽 변형에 대한 좋은 근사법이다. 횡 변형 항이 무시되면 Hooke법칙은 다음과 같이 환원된다.

(132)\varepsilon  = \frac{\sigma }{E}

여기서

  • ε 는 변형
  • σ 는 수직응력
  • E 는 Young 탄성계수
  • 식(10.131) 에서 K 는 사용자-기술 변수이다. K 값을 정하기위해 한쪽은 고정되고 다른쪽은 수직력은 받는 평형에 있는 판을 고려한다. 포아송비는 무시할 만하고 수직응력은 표면상의 한 위치에서 σ 이다. 힘의 균형으로부터 판내의 수직응력은 그 위치를 통과하는 수직선상을 따라 같은 σ 이다. 판의 두께를 h 라하고 휨을 w 라고 표시한다. 변형은 이때 w/h 이고 식(10.132) 에서Hooke 의 법칙은 식 (10.131) 를 준다.

(133)w = \frac{\sigma }{{{E \mathord{\left/
{\vphantom {E h}} \right.
\kern-\nulldelimiterspace} h}}}

이 경우 이라는 것을 가리킨다. 일반적으로 K 는 다음으로 추정된다

(134)K = \frac{E}{L}

여기서L은 탄성벽의 두께에 비교될 만한 길이 규모이다. 또한 실험 측정이나 실제규모의 구조해석 같은 다른 방법으로 부터도 구해질 수 있다.

유체 유동에 대한 이 탄성벽 운동 효과는 벽표면에 체적 소스(또는 싱크)를 추가함으로써 기술된다. VOF 에대한 연속 및 이송방정식들은 탄성막에 대한 절에서 기술된 것과 같은 방식으로 수정된다.

Model Limitations

이 모델은 다른 형상, 크기, 방향 및 물리적 양을 가지는 다수의 탄성막과 탄성벽 물체를 허용한다.이는 대부분의 다른 FLOW-3D 모델들과 같이 사용될 수있다. 예를들면 막과 탄성벽을 통과하는 열전달이 포함될 수있다. 그러나 제약과 제한도 있다.

이 모델에서 탄성막과 탄성벽이 작은 휨을 가져야한다. 막에 대한 작은 휨은 휨이 막의 크기보다 작다는 것을 뜻한다. 한 탄성벽에서 그 변형은 계산 셀 크기와 비교하여 작아야한다. 이 가정은 변형되는 표면의 위치에 실제 차이를 무시하게끔 한다. 다른 말로 막과 변형하는 벽의 형태는 그들의 초기에 지정된 정의대로 계산하는 동안 고정된다. 위에 보여준 바와 같이 유체 운동에 대한 표면 변형의 효과는 변형하는 표면상에 유체 소스 분포를 통해 모델링된다. 계산 결과가 보여질 때, 변형은 의 등고선을 그림으로써 가시화될 수있다. 탄성막/벽은 동시에 다공체이거나 이동체일 수가 없음을 또한 주목한다. 이동체가 막이나 탄성벽과 충돌하면 후자는 충돌모사에서 고정된 강성체로 간주되어서 충격은 단지 이동체의 운동에만 영향을 줄 것이다.

[FLOW-3D 이론] Auxiliary Model/Elasto-visco-plastic & Viscoelastic Model 탄성-점성-소성 및 점탄성모델

어떤 물질은 전체응력 상태에서 압력과 점도외에 추가의 응력이 존재한다: 탄성응력. 점성응력은 물질의 변형율의 함수인 반면에 탄성응력은 물질의 지역 전체 변형의 함수이다. 이런 물질은 탄성응력이 항복 한계까지 증가하며 그 때에 점성유체 같이 거동하는 Bingham 유체, 그리고 지역 변형율에 따라 물질안에서 점성 및 탄성응력이 계속적으로 변하는 점탄성 물질을 포함한다. Bingham 유체의 예는 치약 같은 고상입자들이 작은 변형하에서 서로 “잠겨져” 있고 이 결합이 부서지면 액체같이 흐르는 고상으로 채워진 부유물 형태이다. 중합체 융해와 용매는 점탄성 유체의 예인데 중합체의 긴분자들이 서로 얽혀 급격한 전단력을 받으면 변형에 저항하지만 전단이 천천히 일어날 경우 쉽게 다른 분자들을 지나 미끄러질 수 있다.

점성응력은 순전히 어느 한 시점에 영역 내 모든 점에서의 순간 변형율에 달려있어서 모사 기간 중에 FLOW-3D 에서 저장될 필요가 없지만 탄성응력은 물질 내 한 점에서의 시간에 대한 이력이 현재 상태의 탄성응력을 계산하기 위해 알려져야 할 필요가 있으므로 저장되어야 한다. 더구나 대부분 물질은 자연적으로 등방성이고 Cauchy 응력은 대칭이며 6개의 고유한 탄성응력텐서 성분이 있다. 디폴트로 이들은 이류 양들로 저장되어 탄성응력은 영역전반에서 물질을 따라 이류된다. 최소한의 유동이 예상되는 경우에 탄성응력의 이류는 중지될 수 있다.

단지 탄성응력의 편향 부분이 저장된다; 등방성 부분은 단지 압력이고 이는 FLOW-3D 에서 점성 액체에서와 마찬가지로 연속방정식의 시행을 통해 해석된다. 편향 부분은 물질의 전단 및 인장과 관련된 탄성응력이다. 응력의 등방성 부분은 (거의) 비압축성 물질에 대한 압력과 많이 유사하게 해석된다; 모든 경우에 대해 대부분 물질의 높은 체적 탄성율로 인한 작은 시간간격을 피하기 위해 내재적으로 해석된다. 응력의 편향 부분은 탄성응력(Elastic stress) 솔버의 활성화시킴으로써 내재적 또는 외재적으로 별도로 해석된다.

 Model Formulation 모델형성

FLOW-3D 에 포함된 증분적 탄성응력 모델은 선형Hookean 이론을 사용하는 탄성응력을 계산한다. 비록 이 간단한 구성 방정식은 응력에 선형 반응을 예측하지만 증분 모델로서의 시행은 상당히 비선형적 반응의 예측을 허용하는데 이는 각시간 단계 내의 증분적 변형은 작으므로 각각의 작은 시간 단계내의 반응이 선형으로 근사될 수 있기 때문이다.

탄성응력, 점성응력 그리고 항복의 관계는 선택된 모델에 의존한다.

탄-점소성모델에대해 밑의 그림은 어떻게 전체 응력이 항복응력 한계로 인한 미끄러짐에따라 점성 응력과 탄성 응력의 합이 되는지를 그림을 통해 보여준다. 이와같이 이 모델은 전체응력 상태를 각기 감쇄기와 스프링으로 나타내진 점성응력과 탄성응력의 합으로 예측한다. 이 모델에 의해 기술된 바와같이 이런 재질의 봉이 아주 단기간에 부과되는 일정량의 변형을 겪는다고 가정하자. 이 모델은 변형에대해 선형적으로 비례하는 탄성응력의 상응 증가를 예측한다. 또한, 변형이 부과되는 단기간에 점성응력도 상당히 클 것이며 단지 부과되는 변형이 멈출때 0으로 떨어질 것이다. 더 큰 변형이 탄성응력이 항복응력을 능가하는 점까지 부과되면 재질은 항복 변형되어 액체처럼 흐르기 시작할 것이다.

Pictorial view of Elasto-viscoplastic model, showing the relationship between the elastic and viscous stresses.그림 10.1 탄성과 점성 응력간의 관계를 보여주는 탄성-점소성 모델의 그림도해

점탄성 모델(Oldroyd-B 모델 또는Giesekus 모델) 이 선택되면 점성과 탄성 응력 사이의 도해적 관계가 하기 그림에서와 같이 보여진다. 이 모델에서 스프링과 감쇄기는 직열로 연결되어 물질 내의 임의의 점에서 힘의 전체상태는 점성응력과 탄성응력이 둘 다 같다는 뜻에 주목한다. 그러나 변형과 변형율은 매우 다를 수 있다. 급격히 인장되는 이 재질로 된 봉은 급격한 이동에 의해 변형률이 크므로 큰 점성응력을 가질 것이다. 탄성응력도 마찬가지로 커서 재질은 고체같이 거동할 것이다. 같은 전체 변형이 장기간에 적용된다면 결과적으로 변형율이 아주 작아져서 점성력도 매우 작을 것이다, 마찬가지로 전체 변형이 매우 크더라도 탄성력이 매우 작음에 틀림없다. 이런 상황에서 재질은 아주 점성 액체같이 거동할 것이다.

Pictorial view of viscoelastic model, showing the relationship between the elastic and viscous stresses.그림 10.2 탄성과 점성응력의 관계를 보여주는 점탄성 모델의 도해

변형과 변형율을 전체 응력상태에 연결해 주는 모델은 구성 방정식들이며 FLOW-3D 에서 3가지 선택이 있다(탄성력이 전혀 없다는 표준 가정과 함께).

(103){\boldsymbol \sigma } = -p {\mathbf{I}}+ {{\boldsymbol \tau }_V} + {{\boldsymbol \tau }_E}

Navier-Stokes 시스템 방정식의 해에 필요한 전체 응력상태는 σ = −pI + τV + τE 이며 여기서 p 는 지역압력, τV 는 점성응력의 편향성분, τE 는 탄성응력의 편향성분이다. p τV 는 이미 현재의 모델에서 고려되고 있다; 이 모델의 목적은 τE를 계산하는 것이다.

Elasto-viscoplastic model

점소성 모델에서 각 시간 단계 사이에 작은 변형의 Hookean 가정을 사용하므로 증분적 응력(한 계산 사이클에서 다음까지의 변화)은 증분적 변형에 관련되어 있다.

(104){{\boldsymbol \tau }_E} = \left( {K - \tfrac{2}{3}G} \right)e{\mathbf{I}} + 2G{\mathbf{E}} - 3\alpha K\left( {\Delta T } \right){\mathbf{I}}

여기서

  • τE 는 탄성응력 텐서
  • E 는변형텐서
  • I 는 단위등방 텐서
  • G 는전단탄성울
  • K 는 체적탄성율
  • α 는 선형열팽창계수
  • T 는 지역 열증분이며
  • e 는 체적변형;
  • 이는 재질의 등방성 팽창이나 수축을 특성화하며 변형 e 의 대각성분 합

(105){\textup{tr}} \left( {\mathbf{E}} \right) .

과 같다.

식(10.104) 은 한 선형 모델이지만 각 증분 변형이 선형가정을 따르기에 충분히 작도록 작은 시간 동안에 증분적으로 계산된다면 비선형 과정도 예측될 수 있다. 식(10.104) 은 τE 의 (정의에 의해 τE 는 대칭) 편향 과 등방 성분 둘 다에 대해 쓰여질 수 있다. 편향 성분은

(106){{\boldsymbol \tau '}_E} = 2 G {\mathbf{E '}}

여기서 τE 는 탄성응력의 편향 성분이고, G 는 전단 탄성율, 그리고 E는 증분 변형의 편향 성분이다:

(107){\mathbf{E'}} = \Delta t \left\{ {\raise0.5ex\hbox{$\scriptstyle 1$} \kern-0.1em/\kern-0.15em \lower0.25ex\hbox{$\scriptstyle 2$}} \left[ {\nabla {\mathbf{u}} + {{\left( {\nabla {\mathbf{u}}} \right)}^T}} \right] - {\raise0.5ex\hbox{$\scriptstyle 1$} \kern-0.1em/\kern-0.15em \lower0.25ex\hbox{$\scriptstyle 3$}} \left( \nabla \cdot {\mathbf{u}} \right) {\mathbf{I}} \right\}

여기서 u 는 지역의 물질속도( FLOW-3D에서 계산되는 것 같은)이고 ∆t 는 시간단계이다

식 (10.104)의 등방성 부분은 p = −Ke + 3αK, 이며 여기서 p 는 평균 등방응력 또는 “압력”의 음이고 −tr (τE)/3와 같다. 식(10.106) 을 미분방정식으로 표현할 수 있다:

(108)p = - Ke + 3\alpha K\left( {\Delta T} \right)

(109)\underbrace{\frac{\partial {\boldsymbol \tau '}_E}{\partial t}}_{\substack{\textrm{Change in stress} \\ \textrm{at fixed point in space}}} + \underbrace{\nabla \cdot \left( {\mathbf u} {\boldsymbol \tau '}_E \right)}_{\substack{\textrm{Change in stress}\\ \textrm{due to movement}\\ \textrm{of material}}} = \underbrace{2G {\mathbf{D'}} \left( {x, t} \right) }_{ \substack{\textrm{Change in stress}\\ \textrm{due to stretching}\\ \textrm{or shearing}}} + \underbrace{{\boldsymbol \tau '}_E \cdot {\mathbf{W}} + {\mathbf{W}}^{ \mathbf{T}} \cdot {\boldsymbol \tau '}_E}_{ \substack{\textrm{change in stress due to}\\ \textrm{solid-body rotation}}}

여기서 D는 변형 율 텐서의 편향성분이다. 3차원에서 이것은

(110){\mathbf{ D'}} = {\raise0.5ex\hbox{$\scriptstyle 1$} \kern-0.1em/\kern-0.15em \lower0.25ex\hbox{$\scriptstyle 2$}}\left[ {\nabla {\mathbf{u}} + {{\left( {\nabla {\mathbf{u}}} \right)}^T}} \right] - {\raise0.5ex\hbox{$\scriptstyle 1$} \kern-0.1em/\kern-0.15em \lower0.25ex\hbox{$\scriptstyle 3$}} \left( \nabla \cdot {\mathbf{u}} \right){\mathbf{I}}

이며 W와동 텐서이다

(111){\mathbf W} = \frac{1}{2} \left[ {\nabla {\mathbf{u}} - {{\left( {\nabla {\mathbf{u}}} \right)}^T}} \right] .

식 (10.109) 은 FLOW-3D 에서 탄성응력 성분을 저장하기 위해 스칼라 배열을 이용하여 해석될 수있다; 이 성분들은 이미 물질의 운동과 함께 이류되고 증분적 응력이 각 시간단계에서 응력 텐서 성분에 추가된다. 유사하게, 식 (10.108) 은 다음과같이 다시 쓰여질 수 있다.

(112)\frac{{\partial p}}{{\partial t}} + \nabla \cdot \left( {{\mathbf{u}}p} \right) = - K\dot e + 3\alpha K\left[ {\frac{{\partial T}}{{\partial t}} + \nabla \cdot \left( {{\mathbf{u}}T} \right)} \right]

식 (10.112) 의 이류항(좌측2번째항)은 압력파가, 특히 고체 내에서, 물질이 이류하는 속도보다 훨씬 빠르므로(모든 아음속 모사에서) p 의 시간 미분량에 비해 상대적으로 작다고 가정된다.

응력의 현재값은 물질요소의 과거 이력의 함수이다. 영역으로 들어오는 물질은 영역경계에서 별도로 지정되지 않는한 의미없는 응력을 갖는다. 또한 모사 시작시의 초기조건은 초기조건에 의해 지정되지 않는한 의미없는 응력을 갖는다. 항복효과를 예측하기 위해 Mises항복 조건이 사용된다. 이 조건은

In order to predict yielding effects, the Mises yield condition is used. This condition is

(113)I{I_{{\boldsymbol \tau '}_{\mathbf E}}} = \frac{Y^2}{3}

이며,  여기서 IIτE 는 탄성 응력텐서의 2차 불변항이고 Y 는 사용자 지정변수인 항복응력 한계이다. 탄성응력이(IIτE 에 의해 측정되는) 항복 범주를 벗어나는 물질의 지역에서 탄성응력은 식 (10.113) 의 조건이 만족되도록 완화된다:

(114){\boldsymbol \tau '^{^*}}_{\mathbf E} = \sqrt{\frac{2Y^2}{3I{I_{{\boldsymbol \tau '}_{\mathbf E}}} } {\boldsymbol \tau '}_{\mathbf E} }

여기서 τ*E 는 항복-한계의 탄성 응력텐서이다; 유체의 모멘텀 균형을 위한 Navier-Stokes 방정식에 적용되는 것이 이 텐서이다.

Viscoelastic models 점탄성모델

점탄성 물질에서는 유체내 점성력과 탄성력의 연속적인 균형이 이루어지므로 다른 구성 방정식이 필요하다. 식(10.104) 대신에 항복 조건인 식 (10.113)과 함께 Oldroyd-B 모델에 대한 τE 를 정의하기 위한 구성방정식은

(115)\lambda \frac{D {\boldsymbol \tau '}_{\mathbf E}}{D t} + {\boldsymbol \tau '}_{\mathbf E} = 2 \eta {\mathbf D'}

이다.

여기서 λrelaxation time 이고 ηelastic viscosity 이다

\frac{D {\boldsymbol \tau '}_{\mathbf E}}{D t}  {\boldsymbol \tau '}_{\mathbf E}의  upper-convected derivative 이다. 또는 η 로 정의될 수 있는데 여기서 G 는 고체같이 거동하는 물질의 탄성계수이다.  그러므로

(116)\frac{D {\boldsymbol \tau '}_{\mathbf E}}{D t} = 2 G {\mathbf D'} - \frac{1}{\lambda} {\boldsymbol \tau '}_{\mathbf E}

식(10.116) 은 우측의 추가 선형항을 제외하고는 식 (10.109) 와 같은 형태로 확장된다. 완화시간 λ 의 큰 값은 이 항이 작고 탄성응력이 천천히 완화된다는 것을 뜻한다. 작은 λ는 이항이 크고 탄성응력이 아주 급격히 완화된다는 것을 뜻한다

유사하게 Giesekus 모델에서도 식 (10.115)에 추가 비선형 항이 더해진다:

(117)\lambda \frac{D {\boldsymbol \tau '}_{\mathbf E}}{D t} + {\boldsymbol \tau '}_{\mathbf E} - \alpha \frac{\lambda}{\eta} {\boldsymbol \tau '}_{\mathbf E} \cdot {\boldsymbol \tau '}_{\mathbf E} = 2 \eta {\mathbf D'}

여기서 αdimensionless mobility factor이다. 이 인자는 물질의 점탄성 거동에 대한 추가 맞춤 변수를 제공한다. 다른 변수들은 Oldroyd-B 에 대해서와 같은 의미를 갖는다.

식 (10.117), 를 재정리하여 다음을 얻는다:

(118)\frac{D {\boldsymbol \tau '}_{\mathbf E}}{D t} = 2 G {\mathbf D'} - \frac{1}{\lambda} {\boldsymbol \tau '}_{\mathbf E} + \frac{\alpha}{\lambda G} {\boldsymbol \tau '}_{\mathbf E} \cdot {\boldsymbol \tau '}_{\mathbf E}

Oldroyd-B Giesekus 모델에서 구성방정식 (10.116)과 (10.118)들은 FLOW-3D 스칼라 이류 방정식을 이용하여 해석되고 식 (10.109), (10.116), 또는 (10.118) 들의 추가 항(선택된 모델에 따라)들은 내내적 또는 외재적으로 풀어진다. 다음절은 이를 더 자세히 기술한다.

Solution Method

변형율텐서 E˙ 의 해법은 점성 응력 알고리즘에 사용된 것과 유사하다. E 의 법선성분은 셀 중심에서 계산돠고 전단 성분은 셀 면에서 계산된다. 하기 그림은 전형적인 계산셀과 응력성분 계산 위치를 보여준다. 전단성분은 계산의 편리성과 수치안정성 때문에 셀 면에서 계산된다.

Computational cell showing locations of components of strain rate tensor

Figure 3변형율텐서의 인자의 위치를 보여지는 수치해석 셀

 

자유표면에서 외부유체는 유체면에 무시할 만한 응력을 미치는 기체기포로 가정된다. 그러므로 이런 경계면에서 n·τE = 0이다. FLOW-3D 에서 자유표면과 인접한 셀의 주방향은 알려져 있다. 이 정보로부터, τE 의 관련된 성분은 위에서 언급된 조건을 만족시키기 위해 자유표면에서 0으로 지정된다. 고체벽과 물체와의 경계면에서 벽(또는 물체)과 인접하는 셀의 유체속도와 벽(또는물체)의 속도가 식 (10.110)를 이용하여 변형율을 계산하는데 이용된다.

외재적 방법이 선택되면 탄성응력은 전 시간단계의 정보를 이용하여 식 (10.109), (10.116) 또는 (10.118) 들에서 갱신되고 τE 의 결과값이 다음 단계의 새속도를 계산하기 위해 식 (10.9) (점성 가속도에의 추가항으로) 에서 사용된다.

이 방법은 각 시간 단계에서 단순하고 빠른 장점을 갖는다. 그러나 탄성효과가 다른 효과들을 능가하는 물질들에서는 수치안정성 제약이 아주 작을 수 있다. 그러므로 이런 문제들에서 시간단계의 제약을 없애기 위해 내재적 방법이 선택될 수 있다. 이 방법의 단점은 각 시간단계에서 훨씬 더 큰 계산적 노력이 필요하며 물리적으로 존재하는 진동을 감쇠시킬 수 있다.

내재적 알고리즘이 사용되면 속도 예측은 외재적방법과 같이 계산된다. 그러나 속도가 등방성 응력식(10.112) 과 함께 모멘텀 방정식의 잉여속도(외재적 근사와 가장 최근의 속도의 근사와의 차이)에 의거하여 갱신된다: 

(119)\delta {\mathbf{u}} = \frac{{\Delta t}}{\rho }\left[ { - \nabla \left( {\delta P} \right) + \nabla \cdot \left( {\delta {\boldsymbol \tau '}_{\mathbf E} } \right)} \right]

여기서 τE 의 값은 가장 최근의 값이다.

식 (10.119) 을 해석하는데 이용되는 방법은 Jacobi 방법이다; 이 방법에서 영역 전체에 대해 갱신된 모든 탄성응력텐서는 전반복에서 계산된 속도에 근거하여 해석된다; 인접한 최근에 계산된 셀들의 성분들은 다음 반복때까지 반영되지 않는다. 이는 한 반복에서 다음 반복때까지 속도의 차이가 미리 지정된 수렴공차 보다 작아질 때까지 반복된다. 원하면 이 방법은 점성효과가 아주 큰 물질들에 대한 Jacobi 내재적 점성 응력 방법과 함께 사용될 수 있다. 현재, 내재적 탄성응력모델은 ADI (Alternating Direction Implicit) 나GMRES (Generalized Minimum Residual Solver) 내재적 점성 응력 알고리즘과 함께 사용될 수 없다.

압력 P 는 식 (10.112)을 이용하여 별도로 계산되며 해법 알고리즘은 액체압력에 대해서와 똑 같다. 압력과속도 성분은 식 (10.112) 의 잔여가 미리 지정된 수렴 공차 밑으로 떨어질 때까지 각 반복에서 갱신된다. 수치적 방법은 SOR (Successive Over Relaxation), ADI, 또는GMRES일 수 있다. SOR 알고리즘은 현 반복단계에서 갱신된 인접 값들을 포함하여 최근에 갱신된 값이 사용된다는 점을 제외하고는 이미 언급된 Jacobi 방법과 아주 유사하다. ADI 알고리즘은 동시에 선택된 방향의 줄 전체의 셀을 해석하며 GMRES 알고리즘은 완전히 결합된 방정식 시스템을 해석한다.

Computation of parameters 계산변수

탄성계수 G , 항복응력한계 Y , 완화시간 λ 그리고 유동인자 α 는 유체 1과 2모두에 정의 될 수 있다. 추가로 표로 시간에 따라 변하는 데이터를 제공할 수도 있다.

두 점탄성 유체가 한 계산 셀내에 있을 때 그 셀내의 변수는 두 유체 물성치의 가중 선형보간 값을 이용한다. 표상의 온도의존 데이터에 대해서는 FLOW-3D 내의 다른 온도 의존 양에서와 마찬가지로 구간 선형 보간이 두 데이터 점 사이에서 가정된다.

[FLOW-3D 이론] Auxiliary Model/High Concentration Granular Media Model 고농축 입상매질모델

여기서 고농축 입상유동이라는 명칭은 입상율의 체적율이 50% 또는 그 이상일 때를 뜻한다. 고농축시 강한 결합이 고체입자와 주변유체간에 작용하여 이 혼합물은 단일 혼합유체로 잘 근사될 수 있다. 이 혼합물은 이를 순수기체 지역과 분리시키는 자유 표면을 가지는 비압축성 유체로 간주된다. 두 물질의 속도 차이로 인한 혼합물 내의 2상 효과는 표류-플럭스 근사를 이용하여 고려된다.

두 종류의 입상 유동 모델이 있다. 하나는 고체 입자를 둘러싼 연속유체가 가스이며, 다른 하나는 연속유체가 액체이다(혼합물은 슬러리로도 불린다).

일반적으로, 혼합유체는 불균일 밀도를 가진다. 유체 혼합물의 밀도는 비압축성 가정 때문에 체적 혼합물 유동시에 변할 수 없다. 그러나 밀도변화는 표류-플럭스 모델에서 기술된 바와같이 혼합물 내에서 고체와 가스간의 상대유동 때문에 발달할 수 있다.

탄탄한 패킹은 입자들이 함께 압착될 때 밀려난 유체가 제거되어야한다. 기체/입상 혼합 모델에서 혼합체적의 손실은 혼합물 자유표면에서의 기체의 손실로 간주된다. 기체는 기체와입자의 상대운동으로 인해 자유표면상에서 주변공간으로 빠져 나갈 수있다. 모든 제거된 기체는 주위공간 지역으로 전달된다. 어느 경우던지 모든 제거된 기체는 주위 기포로 전달된다. 액체/입상 혼합 모델에서 고체물질을 다지면 자유표면을 가지는 순수액체 공간이 만들어진다.

이산 고체는 단지 밀도가 순수 고체보다 작은 지정된 패킹 한계까지만 쌓여질 수있다. 균일 크기의 좁은 패킹을 가지는 구형입자의 일반적 고상 체적율 값은 is fspk = 0.63이다. 표류-플럭스 모델에서 발생할 수 있는 최대 고상율을 선택할 수 가있다. 이 제한은 임계 체적율에 도달할 때 표류 속도를0으로 지정함으로써 적용된다. 이 처리는 고 체적율을 가지는 분산 물질 내 발생하는 입자간의 상호작용의 경험적 설명인 Richardson-Zaki상관식에 의한 표류 속도 제약과 유사하다.

고농도의 입상 물질내 전단응력은 분산된 고체를 운반하는 유체의 점성 전단응력 보다 훨씬 크다고 알려져있다. [Bag05] 의 1941년도 연구에서 시작하여 광범위한 양의 연구가 요약되어 있으며 [Mih99]에 의해 더 큰 실험범위로 확대되었다. 고농도상태에서 전단응력의 발생의 주요한 원인은 입상간의 충격력(즉, 충돌)이다. 두 번째, 그리고 일반적으로, 작은 부분이 유체에 영향을 주는 분산된 고체에 기인한다. Mih의 유효 동적점성 표현식([Mih99] 의 식 39로부터 변경된)은,

(91){\mu _{\rm{eff}}} = 7.8\mu \left( {\frac{{{\lambda ^2}}}{{1 + \lambda }}} \right) + {\rho _s}\left( {\frac{{0.015}}{{1 + 0.5\rho /{\rho _s}}}} \right)\left( {\frac{{1 + e}}{{{{\left( {1 - e} \right)}^{0.5}}}}} \right){\left( {\lambda d} \right)^2}\left| {\frac{{du}}{{dy}}} \right|

여기에서:

µ ρ 는 연속적 유체 점도 및 농도는,

ρs 는 모래농도

e 는 입상 충격과 관련된 반발계수,

d 는 입자 직경이며

λ 는 고상체적율 fs 로나누어진 최대 고상 체적율 fsmx = 0.63 의 함수이다.

(92)\lambda  = \frac{1}{{{{\left( {1.032{{f_s^{{\text{mx}}}} \mathord{\left/
{\vphantom {{{\text{f}}_{\text{s}}^{{\text{mx}}}} {{f_s}}}} \right.
\kern-\nulldelimiterspace} {{f_s}}}} \right)}^{1/3}} - 1}}

물리적으로, λ = d/S 이며 여기서 S 는 입상 중심들의 평균거리에서 입상 직경 d 를 뺀 값으로 정의된다. 입상이 접촉하면(즉, 패킹되면) S=0 이고 λ 는 무한대가 된다. 현재의 목적을 위해 단순 전단율 du/dy 는 변형율 eij 의 크기로 대체되고 모래에 대한 일반 반발계수 0.7은 타당한 일반값으로 가정된다. 이렇게 변경하면 식(10.91) 은 다음으로 축소된다.

(93){\mu _{\rm{eff}}} = 7.8\mu \left( {\frac{{{\lambda ^2}}}{{1 + \lambda }}} \right) + 0.066{\rho _s}{\left( {\lambda d} \right)^2}\left| {{e_{ij}}} \right|.

 

이 식은 ”shear thickening” 점도이다. 농밀화 특성은 saltation 이라고 불리는 과정인 입상 사이의 충격 영향으로부터 발생한다. 즉, 가라앉거나 느리게 움직이는 입상들의 층 위에서 빠르게 움직이는 입상들은 서로 부딪혀서 저속의 입상들을 위의 흐름으로 이동하게 하며 여기서 저속 입상들은 더 많은 충돌을 통해 흐름으로부터 모멘텀을 흡수한다.

충격 점도항(식(10.93)의 2번째항)은 일반적 대와류 모사법(LES) 점도 형태를 가진다. 식(10.93) 의 경험적으로 유도된 표현은 난류유동 조건을 포함하므로 입상유동 모델에 난류이송 모델을 사용하는 것은 불필요하다(아마 일치하지 않을 것임).

유체/입상 혼합물의 2상 유동을 모델링하는 중요한 요소는 모래가 적용된 압력이나 체적력을 받을 때 패킹되고 운동에 저항하기위한 한 방편을 갖는 것이다. 이 목적을 위해 발포되는 모래의 관찰에 기반한 유동저항 모델이 이용된다. 입상물질이 각 개의 입상물질이 서로 접촉하기 시작하는 밀도로 패킹되면 혼합물이 유동하기에 더 어려워진다. 이런 상태를 일종의 기계적 방해라고 하며 전형적인 체적율 fsjam = 0.61를 가진다. 고체의 패킹에 해당하는 아직 더 높은 고상밀도에서 입상은 이웃 입상들과 접촉하며 유동이 불가능할 것이다.

Bagnold 에 의하면 바닥에 가라앉고 패킹된 입상들은 단지 모래유동이 한계값 uthrs 보다 큰 속도를 가질 때 바닥에서 벗어나거나 모래 유동 내로 포함된다. 중력 g 와 입상 직경 d 를 포함하는 경험적 연구와 한 근사 이론에 의해 한계속도는

(94){u_{\rm{thrs}}} = 1.41{C_{\rm{drg}}}\sqrt {d\left| \textbf {g} \right|{{\left( {{\rho _s} - {\rho _a}} \right)}/{{\rho _a}}}}

여기서Cdrg는 1이다. 이식에서 Cdrg 의 포함은 아마 추후에 입상간의 간섭효과를 고려하는 한계 속도 조정을위해서 일 것이다.

유체를 포함하는 망 요소 내의 패킹 저항은 요소내 유체속도가 이 한계속도 보다 크거나 고상물체의 체적율이 기계적인 방해에 대한 필요보다 작으면 0이라고 가정된다. 이 저항값보다 큰 체적율 및 한계값보다 작은 속도에 대해 유동저항은 다음과 같이 주어지는 음의 가속도로 유체 모멘텀방정식의 우편에 추가되며,

(95)- {S_{\rm{drg}}}\sqrt {\frac{{\left| g \right|}}{{\,4d}}} \left( {\frac{{{u_{\rm{thrs}}} - \left| {\vec u} \right|}}{{{u_{\rm{thrs}}}}}} \right)\left( {\frac{{{f_s} - f_s^{\rm{jam}}}}{{f_s^{\rm{mx}} - f_s^{\rm{jam}}}}} \right)\vec u

 

여기서 Sdrg 는 차수가 1인 무차원 계수이다. 이 저항은 한계속도를 초과하는 유동속도 및 또한 모래의 체적율에 비례한다.

고상물질의 체적율이 약0.99fspk 의 값에 도달하거나 초과하면 유동속도는 0으로 되고 물질은 완전히 패킹됬다고 간주된다. 이의 예외가 내부로 향하는 표면법선과 체적력(즉, 중력 같은)의 방향간의 각도로 정의된 안식각보다 큰 구배를 가지는 표면요소에 적용된다. 이런요소는 유동저항이 없다.

입상 유동에 대해 현재의 모델에 더 현실성을 주는 안식각 및 근사 패킹한계의 예외가 있다 한 예외는 안식각보다 약 2도정도 큰 값의 이동각을 포함하는 것이다. 표면이 패킹되고 안정화 되있을 때 표면구배가 이동각보다 클 때까지 유동하지 못하나 구배가 안식각보다 큰 동안에 한해서 계속 유동할 것이다. 이동각은 유동이 발생하기 전에 극복되어야 할 일종의 정지 마찰력 같이 거동한다. 유사한 “정지마찰” 개념이 이전에 패킹된 지역을 굳지않게 하기 위한 혼합유체의 내부에 이용된다.

고체 벽경계에서 입상물질이 경계로부터 분리(표류)될 수 있어 고상 체적율이 50% 보다 작아질 수 있다. 이런 가스/입상 혼합물의 경우 특수한 취급이 필요하다. 이 경우 초기에 자유표면이 없어서 일종의 “공동” 과정을 통해 생성되어야 한다. 50% 보다 작은 고상율이 자유표면으로 인식되지 않은 지역에서 발생하면 그 격자요소 위치의 압력이 외부의 가스압으로 완화되고 이 격자요소는 압력이 정의되고 비압축성이 적용되지 않는 특수요소로 표식된다. 이로써 요소는 열려서 새로운 순수가스 지역이 된다.

최대 패킹 체적율 fspk = 0.63 은 변할 수 있는 입력변수이다. 변경이 되면 fspk = 0.63에서도 그랫듯이 마찬가지로 상응하는 고상율 fsjam fsmx 가 자동적으로 변경된 fspk 와 같은 비례값을 가지도록 조절된다.

여기서 기술된 입상모델은 패킹된 고체 내를 통하는 유체의 자세한 유동을 다루지 않는 것 같은 약간의 제약이 따른다. 그럼에도 불구하고 많은 실제 양상에 유용한 결과를 보여준다. 예를들면, 금속 주조시 모래코어 제조에 적용할 때 정성적으로 실험 데이터와 잘 일치하고 있다.

Sand Core Blowing Model 모래코어 발포모델

금속주조 응용을 위한 모래 코어의 생성은 일반적으로 원하는 코어의 형태의 공동을 가지는 상자로 모래/공기의 혼합물을 불어넣어 만들어진다. 코어 전체 체적으로 균일하게 모래충진을 하기 위해 보통 다수의 공기 배출구가 상자 내에 주어진다. 이 배출구들은 많은 경우에 모래입자들 위에 입혀진 점결제의 열 또는 화학적 경화를 효과적으로 하기 위해 공기나 촉매가스가 모래 내로 유입되도록 해준다.

모래코어 모델은 특정하게 정의된 공기 배출구가 보충된 입상유동 모델의 응용이다. 이 모델의 출력은 공기 배출구의 위치 및 크기를 고려한 코어 상자 내의 충진 형상의 이력을 포함할 뿐만 아니라 완성된 코어 내에 존재할지도 모르는 밀도변화의 정보도 제공한다.

기본유동 모델의 상세내용은 High Concentration Granular Media Model 의 기술에서 찾을 수 있다. 모래코어 발포에 필요한 배출구에는 두 가지 형태가 있다; 다이 내부의 공기 포켓들과 다이 외부의 공기를 직접 교환하는 정규 배출구와 또 하나는 모래 내의 포켓으로부터 모래에 의해 덮혀진 출구로 가는 간접적 공기 유동의 결과인 출구가 그것이다. 후자의 배출구는 “일반적 배출”이라고 불린다.

정규 배출구는 둘 다 계산 격자내 순수 가스지역(기포)과 격자외부의 압력 지역간의 가스를 교환하는 기능을 가진면에서 밸브와 유사하다. 배출구는 위치, 유동면적 A, 외부압력 Pext 그리고 손실계수 Vc. 에 의해 특성화 된다. 순수가스 지역이 배출구를 둘러싸면 압력차이에 따라 가스를 얻거나 잃을 것이다. 배출구를 통한 체적유량 Qv, 은 다음과 같이 정의된다.

(96)Q_v = V_c \cdot A \cdot \sqrt {\frac{{2\left( p - P_{\rm{ext}} \right)}}{\rho }}

제곱근인자는 압력강하로 인해 배출구를 통과하는 유동속도에 대한 Bernoulli 근사이다. 변수 p ρ 는 출구에 인접한 가스의 압력과 밀도이다. 이는 간단히

(97)Q_v = V_{c2} \cdot \sqrt {p - P_{\rm{ext}}}

로 표현되고

여기서 배출구 계수, 배출구 면적, 제곱근 표시 안에 있는 인수 2 그리고 가스밀도 모두 배출구 계수 Vc2안에 결합되어 있다.

전형적인 모래코어 발포 응용에서, 모래 막은 뚜렷한 모서리를 갖는 틈이나 구멍 같은 출구를 가지는 배출구 채널상에 놓여진다. 이런 출구를 통한 유동손실에 대한 합리적 추정은 이다. 이와 같이 배출구 계수는

(98)V_{c2} = A \sqrt {\frac{1}{{2\rho }}}

로 주어지며

여기서A는 배출구(즉, 막)의 실제 개방된 유동 면적이다.

임의의 수의 배출구가 정의될 수 있다. 배출구 위치를 포함하는 계산 격자 셀 안에서의 혼합 유체분율이0.5보다 크면 배출구는 유체에 의해 막혀있고 더 이상 배출구를 통한 가스유동이 없다고 가정된다.

모래코어 발포에서는 심지어 패킹된 모래에서도 가스의 유동이 있다. 코어 상자 안에 배출구가 있거나 상관 없이 상자 내의 모든 공기기포가 배출되도록 해주는 이러한 가스유동의 단순한 표현은 일반 배출구라고 불리어진다. 일반 배출구는 이 배출구의 외부압력 Pg 가 모든 정규 외부압력의 평균이고 일반배출구 계수 Vg 는 모든 정규 Vc2 값의 이라는 것을 제외하고 각기 배출구로서의 같은 공식을 갖는다. 추가로 일반배출구의 유동계수는 모래 통과시의 평균 유동 손실에 더해진 모래에 의한 막힘으로 인한 배출구 면적의 감소를 나타내는 인수 Vg 로 곱해진다.

단순한 추정을 일반배출 계수 승수 Vg 에 대해서 할 수 있다. 표준 배출구에서 배출구를 통과하는 유동속도는 유동손실 인자와 함께 Bernoulli 근사를 사용하여 추정된다. 일반 배출구에 대해 공기는 실제 배출구를 통해 빠져나가기 전에 다공성 모래를 통과하여야 한다. 모래의 투과성은 수직 배출구에대한 Bernoulli 표현과 직접 비교를 가능케한다는 가정인 형상 손실이 주를 이룬다. 이 가정하에 배출구로 모래를 통과하는공기속도는 다음식에 의해 표현된다,

(99)\left| u \right| = {\left( {\frac{{d{\left( {1 - f_s^{\rm{pk}}} \right)}^2}}{6f_s^{\rm{pk}}L}} \right)^{1/2}}\sqrt {\frac{2\left( {p - P_g} \right)}{\rho}}

이 표현에서 p 는 다이 내의 공기압력이고, Pg 는 외부압력, ρ 는 공기밀도, d 는 모래입자의 직경은 공기 주머니와 배출구간의 평균거리 그리고 fspk 는 완전히 패킹된 모래의 체적율이다. 패킹된 고상율을 포함하는 인자들은 모래 내의 다공도를 고려한다. 단위 시간당 이 일반 배출구를 통과하는 공기의 체적Q는 수직 배출구의 수 Nv 와 모래에 의한 막힘으로 줄어든 배출구 면적 Av 와 상기 속도를 곱한 값이다.

(100)Q = N_v{A_v\left( {1 - f_s^{\rm{pk}}} \right)\left| u \right|}

이 비율을 수직 배출율과 비교하면 식. (10.97) 와 (10.98)는 수직배출 손실계수들의 합계를 곱해야 하는 모래를 통과하는 유동과 관련된 유효 추가 손실계수를 나타낸다.

(101)V_g = {\left( {\frac{{2{{\left( {1 - f_s^{\rm{pk}}} \right)}^4}}}{{3f_s^{\rm{pk}}}}} \right)^{1/2}}{\left( {\frac{d}{L}} \right)^{1/2}}

최대 모래패킹율인 fspk = 0.63, 통상 입상 직경인 0.02 cm 그리고 L = 1.0 cm 의 추정치에 대해 이 표현은 Vg = 0.02값을 준다. 이 값은 계산적으로 합리적인 결과를 주는것으로 알려져 있다.

Slurry Model 슬러리모델

 

고상입자와 액체의 혼합물은 가끔 슬러리라고 불려진다. 슬러리들은 제조과정에서 많이 볼 수 있고 자연적으로는 폭우에 따른 산의 협곡에서의 석편류에서 나타날 수있다. 슬러리에서의 중요한 양상은 (운반체인 유체와 고체간의 밀도차이가 상대적으로 작기 때문에) 서로 부딪혀 분리되는 전단 유동내의 입자들의 많은 충돌에 의해 발생하는 분산 압력이다. 충분히 큰 전단응력 하에서 분산압력은 입자들이 침전하고 패킹되는 것을 방지한다. 마찰 각 Θfriction 은 입자충돌시 발생하는 접선과 법선응력을 연결시켜 준다. 일반적으로 이 각도는 안식각 보다 2~8도 큰 정도이다. 더 큰 마찰각은 분산압력을 감소시킨다.

[Bag05] 는 모래와 공기의 유동과 관련하여 분산압력의 개념을 소개했다. 높은 고상밀도에서 고체/액체 혼합물의 주 점도 메커니즘은 입자간의 충돌에 의한 것이며 전단유동 내 접선 모멘텀의 전달에 매우 효과적이다. 그러나 입자가 충돌할 때 서로 반발하고 접선 방향 모멘텀 교환에 추가로 전단에 수직한 방향으로 서로로부터 멀어져간다. 또한 입자가 그들 사이에서 멀리 떨어져 나가고 큰 공간을 만드는 경향이 전단유동의 초기에 발생할 때 Reynolds팽창이라고 불려진다. Bagnold 는 일반적으로 충돌에 의해 생성된 전단응력은 항상 상응하는 수직응력을가지는 것을 인식했고 이는 고상입자들을 흩어지게 하는 경향이 있으므로 분산 압력이라고 불렀다. 이 둘 사이의 관계는 간단히

(102){\tau} = {tan(\Theta_{friction})P_{dispersive}}

로 나타나는데 여기서 τ 는 접선 충돌응력, Pdispersive 는 상응하는 분산압력 또는 수직응력이다. 마찰각 Θfriction 는 충돌하는 입자간에 작용하는 마찰응력과 관련된 입상물질의 물성치이다. 이 각도는 입상이 가스보다 액체에 의해 둘러싸여 있을 때 추가된 점도 유체의 영향때문에 더 크다.

[FLOW-3D 이론] Numerical Approximations 수치근사 – Notation 표기

Numerical Approximations 수치근사

Notation 표기

지배방정식을 수치적으로 해석하는데 이용되는 유한 차분망은 폭dxi, 깊이dyj그리고 높이dzk의 직교 셀들로 구성되어 있다. 활성화된 망 지역은 색인i로 명명된 x방향에서 IBAR셀, 색인j로 명명된 y방향 에서 JBAR셀, 그리고 색인k로명명된z방향 에서 KBAR셀을 갖는다. 이 지역은 망 경계조건을 지정하기 위해 사용되는 경계셀이나 가상 셀의 층에 의해 둘러싸여 있다. 이와 같이 전체망에서 통상 전체 셀의 수는 (IBAR + 2)(JBAR + 2)(KBAR + 2)이다. 그러나 주기적 또는 지정된 경계조건이 한 주어진 좌표 방향에서 주어지면 경계셀에 하나의 추가 층이 그 방향에서 사용된다. 이 사실은 차원 설정을 지정할 때만 명심되어야 한다. 전처리는 자동적으로 모든 경계조건을 만족시키는 데 필요한 경계 셀의 수를 초기화할 것이다. 하기 그림은 셀 표시 명명법의 도해이다.

그림 10.11 망 배열및 표식 관례

유체속도와 압력은 밑의 그림에서 보이는 전형적인 셀에서 보여준 바와같이 엇갈린 망위치에 위치한다: x방향에 수직한 셀 표면의 중심에서의 속도u와 면적 분율Ax , y방향에 수직한 셀 표면의 중심에서의 속도v와  면적 분율Ay  그리고z방향에 수직한 셀 표면의 중심에서의 속도w와  면적 분율Az 압력(p), 유체분율(F), 체적분율(VF ), 밀도(ρ), 내부에너지(I), 에너지(q),소산(D),  그리고 점도(µ)의 난류양은 셀 중심에 있다.

여기서 사용된 유한 차분표기는 분수 색인 값이 사용될 수없는 코드에서 사용되는 표기에 상응한다. 관례는 모든 분수 색인은 가장 가까운 정수로 귀결된다. 예를들면, 셀(i, j, k) 과 (i + 1, j, k) 사이의 셀면상에 위치한i + 1/2에서의u속도는 로 표기된다. 윗 첨자nn번째 시간단계 값을 뜻한다. 유사하게,

  (10.311)

ρ, I, q, D, µ 와 유사하게

  (10.312)

면적및 체적분율도 다음과같은 표기로 나타난다.

  (10.313)

자유표면이나 유체경계면이 존재할 때 비어 있거나, 표면이 있거나 또는 한 유체가 가득찬 이러한 셀들을 구별하는 것이 중요하다. 정의에 의해 표면 셀은 유체1을 가지고 적어도 한개의 비어 있거나 유체 2로 가득한 인접(i ± 1,j ± 1,k ± 1에서)셀을 가지는 셀이다. 1보다 작은 F값을 가지며 비어있는 인접 셀들이 없는 셀은 1유체 문제에서 가득찬 셀이라고 간주된다. 표식NFi,j,k는 셀을 분류하기 위해 그리고 표면 셀의 경우에 어느 인접셀이 표면에 내부로 향한 수직인 방향으로 놓여 있는지를 지정하기 위해 사용된다.

   (10.314)

NF는 유사하게 두 유체 사이의 표면의 방향을 표시하는데 사용된다.

Outline of Finite Difference Solution Method 유한차분 해석 해법의 개요

시간의 한 증분δt을 통해 해가 전진하는 기본 과정은 세단계로 구성되어 있다. :

1. 모멘텀 방정식(10.9)의 외재적 근사가 모든 이류, 압력 그리고 다른 가속도들에 대한 초기 조건이나 전시간 단계의 값을 이용하여 새로운 시간 단계에서의 속도에 대한 첫 근사를 계산하기 위해 사용된다.

2. 연속방정식(10.1), (10.5), 또는 (10.8)을 만족시키기 위해 내재적 선택이 사용되면 압력은 반복적으로 각 셀에서 조절되고 각 압력변화에 의해 유발된 속도의 변화는 단계(1)에서 계산된 속도에 추가된다. 한 셀내에 야기된 압력의 변화는 6개의 인접셀 내의 균형을 변동시키므로 반복이 필요하다.압축성 문제에 대한 상태방정식을 만족시키기 위해 외재적 계산에서 각 셀내에서 한 반복이 그래도 필요할 수도 있다.

3. 마지막으로, 자유표면이나 유체 경계면이 존재할 때 새 유체의 형상을 주기 위해 식(10.19)을 이용하여 갱신되어야 한다. 압축성 문제에서 밀도, 식(10.1)와 에너지, 식(10.21)은 이류, 확산 및 소스과정을 반영하기 위해 갱신되어야 한다.  난류량과 벽온도 또한 이 단계에서 갱신되어야 한다.

이와 같은 단계의 반복으로 요구되는 시간단계를 통해 해가 전진할 것이다. 물론 각각의 단계에 모든 망, 구조물, 자유수면 경계에서 적절한 경계조건이 부여되어야만 한다. 이러 단계와 경계조건의 상세 내용은 다음의 세부항목에서 주어진다.

Momentum Equation Approximations 모멘텀방정식 근사

식(10.9)의 유한 차분근사 일반적 형태는(분수 색인 관례를 기억하며) 다음과 같다.

  (10.315)

여기서 예를 들면, 다음과 같이

  (10.316)

   (10.317)

이류, 점성 그리고 가속도항들은 분명한 의미를 갖는다, 즉, FUX는x방향에서의u의 이류 유속량; VISX는x성분 점성가속도; BX는x방향에 수직한 배플에 대한 유동손실; WSHX는x방향에서의 점성벽 가속; 그리고Gx는x방향에서의 중력, 회전 그리고 일반 비관성가속을 포함한다.

First-Order Method 1차적 방법

가장 간단한FLOW-3D유한 차분근사는 시간및 공간 증분에 대해 1차적으로 정확하다. 이 경우 이류및 점성항은 모두 속도에 대해 전시간 단계(n) 값을 이용하여 평가된다. 벽전단응력은 하기에 기술된 바 와 같이 내재적으로 평가된다. 시간단계n+1 에서 압력은 일반적으로 사이클 초기에 알 수 없으므로 이 방정식은 직접n+1시간단계의 속도를 직접 평가할 수 없고 연속방정식과 결합되어야 한다. 한 해의 첫단계에서 이 방정식들에서의pn+1의 값은 새 속도를 위한 첫 추측을 얻기 위해pn로 대체된다. 외재적 근사에서는 식(10.315)에서 압력구배가 시간n에서 평가되므로p에 대한 추가 조정이 un+1의 평가에 영향을 미치지 않는다.

그림 10.13 U모멘텀에 대한 유한 차분근사에 사용된 x-z평면에서의 유한체적(점선)

식(10.315)에서의 다양한 가속도항에 대해 어느 특정한 근사도 수치적으로 안정된 알고리즘에 이르는 한 상대적으로 중요하지 않다. 그러나 불균일한 셀의 크기를 가지는 망에서의 근사에는 특별히 주의를 기울여야 한다. 이 문제는 다른 곳(즉 참고[HN81])에서 논의되고 있으나 가끔 간과되는 문제 때문에 여기서 반복된다. 데카르트좌표를 사용하는 원래의MAC방식 에서 사용된 근사과정을 고려하자 [HW65], [WHS+66].  또한 간단하게 하기 위해 모든 체적및 면적분율이 1이라고 가정한다. MAC방식에서 연속및 모멘텀방정식, 식(10.8)와 (10.9)은 우선 결합되어 이류 유속항들은 발산형식, 즉 uu대신에 ∇uu로)으로 쓰여질 수 있다. 이와 같이, 예를들면, FUX는uδu/δx라기보다δu2/δx일 것이다. MAC방식 에서 발산형태가 선호되는 이유는 차분근사에서의 모멘텀의 보존을 간단히 확실하게 해주기 때문이다. 이는 상기 그림에서 점선으로 표시된 것과 같이ui에대해 사용된 유한 체적의x-z단면을 고려함으로써 보여질 수있다. 발산 형태를 가지고Gauss정리를 이용하여 유한체적 내의FUX적분값을 체적 경계면에서의 유속량으로 전환될 수 있다. 이와 같이 한 유한 체적을 떠나는 유속량을 인접 유한체적이 받아들임으로써 이류 시의 보존이 확실하게된다.

이 개념은 균일한 망의 사용을 위해 개발된 원래의MAC방식에서는 잘 작동했다. 그러나 불행히도 불균일 망에서는 보존이 자동적으로 정확성을 의미하지는 않는다. 이를 보기 위해 상류 또는 공여-셀 차분근사가 FUX =  에 대해 사용된다고 가정하는데 이는 조건적 안정 알고리즘을 제공한다. 이 공여 셀 근사는 다음과 같다.

  (10.318)

 (10.319)

    (10.320)

식(10.318)의 정확도를 체크하기 위해FUX가 평가되는x-위치에 대한Taylor급수 전개를 하면 아래로 평가된다. (u-속도가 양이라고 가정하여)

  (10.321)

0차항에서의 계수는 셀 폭이 같지않은 한, δxi = δxi+1, 부정확하다. 다른 말로 변동 망은 보존차분 근사차이를 한 차수 감소시키며 이 경우에 부정확한 0차수 결과가된다. 차라리 공여셀보다 중앙 차분 근사가 사용되었더라면 결과는 균일 망에서와 같이2차가아니라 1차수에 정확하였을 것이다.

실제로 변동 망들이 항상 덜 정확하다는 것이 반드시 위의 해석에 기인하는 것은 아니다. 예를 들면 이들은 유동변수들이 급격히 변하는 지역에서 개선된 해상도를 가능하게 한다. 그럼에도 불구하고 변동 망 사용시에는 상당히 주의를 기울여야 한다. 예를 들면 근사 차수의 감소를 최소화하기 위해서는 셀 크기의 점진적 변화를 사용하는 것이 최선이다. 또한 변동 망에 적용될 때 공식적인 정확도를 잃지않기 위한 다른 근사를 구할만한 가치도 있다. 이점에서 이류항의 보존형태가 정확도를 잃는 이유는 유한 체적이ui의 위치에 대해 중심에 있지 않기 때문이라는 것에 주목한다. 유한 체적의 이동은식(10.321)에서 계산된 정확도의 직접적인 감소로 이어진다.

FLOW-3D에서 수정된 공여 셀 근사는 변동망에서 정확도를 유지하도록 개발되어 왔고 망이 균일할 때 보존차분 형태로 축소된다[HS85]. 이 방법은 비보존 형태의 이류 유속항 uu,을 근사하는데 이는 위에 언급된 보존적 근사의 원천적 어려움 때문에 필요하다. 이 근사에서 공여셀과 중앙 차분근사를 각 근사의 상대적 양을 조절하는한 변수를 갖는 하나의 식으로 결합하는 것이 또한 가능하다. FUX = 에 대한 이 근사의 일반적 형태는 다음과 같다.

  (10.322)

  (10.322)

(10.322)

망이 균일하면 이 근사는α = 0일때 공간적으로 2차 정확도를 가지는 중앙 차분으로 축소된다. α = 1이면 1차 정확도를가지는 공여셀 근사로 돌아온다. 어느 경우에나 이 방식은 변동망에서 정확한 0차 식으로 축소된다. 균일망에서는 이류 유속항 근사가 보존 근사형태∇uu로 환원되도록 보여질 수 있다.

식(10.322)에 있는 기본 개념은 하류값 보다 유속양이 들어오는 상류값을 가중하는 것이다. 가중 인자는 각기 상류방향과 하류방향에 대해 (1 + α) 와 (1 − α)이다. 1차 근사방식의 이 형태는 식 (10.315)에서 나타나는 모든 이류 유속항에 대해 사용된다. 모멘텀 방정식의 모든 가속도항은 표준 중앙차분에 의해 근사된다.

Second-Order Method 2차적 방법

위에서 기술된 근사는 시간증분δt의 1차멱수에 그리고 α ̸= 0이거나 망이 불균일하면 공간증분dx, dy dz의 1차 멱수에 비례하는 단절 오차를 가진다. 이러한 1차 근사의 장점은 단순하고 계산적 안정성을 유지하는게 쉽다는 것이다. 그러나 어떤 경우에는 너무 고비용이 소요되어 정확한 1차 정확 해석을 위한 필요한 망 해상도를 사용할 수 없다는 것이다. 이런 경우가 발생하면 이류나 점성 가속도들에 대해 2차 정확 근사를 사용하는 것이 유용할 수 있다.  모멘텀 방정식에 대해 두 가지의 선택적 2차 근사 방식이 입력데이터로 요청될 수 있다. 첫 방식의 본질은 이류와 점성 서브루틴을 두번 통과하는 것이다. 첫번째 통과시 공여셀 변수α = 1.0을 가지는1차 방식이 사용된다.  이들 새 속도들은 전시간의 속도 열들에 저장된다.  두번째에서는 1차계산이 변수α = −1.0로 지정된 후 반복된다. 마지막으로 두 계산의 결과가 새로운 시간에서의 2차근사를 얻기 위해 평균된다.

이 근사는 첫번째는 시간n의 속도를 이용하고 두번째는 시간n+1에서의 속도를 위한 (1차)근사를 하므로 시간에 대해 2차이다.  이때 평균은 시간 n+1/2를 가지는데 이는δt에 대해 2차이다. 마찬가지로 첫번째에α = +1.0를 사용하고 두번째에서는α = −1.0를 사용하면 평균α-값 0이되며 이는 망이 균일하면 dx, dy, dz에대해 2차이다(2차근사식에 대한 더 완전한 해석에 대해서는[Hir78]를 참조하라).

이 알고리즘은FLOW-3D에서 이용 가능한 세가지 이류 방식 중에서 가장 덜 수치적으로 확산적이다. 그러나 이는 표준 상류차분 방식에서 유동 교란이 원래 위치에서 하류로 전달되는 것을 확실케 해주는 전달 특성을 가지지 않는다. 추가로, 이 방식은 가장 CPU 집중적이다. 마지막으로, 이 방식은 천이적 자유표면 유동에서 수치적 불안정성을 가끔 발생시킨다.

Second-Order Monotonicity Preserving Method 2차 단조 보존방식

FLOW-3D의 다른 고차 이류 방정식은 2차 단조 보존 풍상(upwind) 차분기법[VanLeer77]에 기반을 둔다. 이는 원래의1차 이류기법 만큼 강력하다. 대부분의 경우에 차이는 무의미 하지만 1차방식 보다 약간 더 많은 CPU시간이 소요된다.

단조보존방식은FLOW-3D에서 모멘텀 이류뿐만 아니라 밀도, 에너지 그리고 유체분율 이류를 근사하는데도 적용될 수 있다. 고차 이산 방식은 각 좌표 방향에서의 이류량들에 대한 2차 다항근사를 이용하여 유도된다[HB91]. x방향에서 이류된 한 변수 Q에대해 셀 면을 통과한 유속값Q*은 다음과 같으며,

 (10.325)

여기서

    • Qi는 셀 중심에서의 값
    • C는Courant수,
    • δxi는 셀 크기이며,
    • A는 셀내의 위치에서 의 1차 미분에 대한 2차근사이다.

(10.326)
계수A는 이런 미분들이 2차로 정확하다면 선형 보간에 의해 두 인접1차미분들로부터 쉽게 계산할 수 있다. 후자는Qi위치들 사이의 중간에서 미분을 구함으로써 얻어질 수있다; 예를 들면,
  (10.327)
위 식에서 Qi Qi + 1사이의 점Q의 2차 정확1차 미분이다. 이런 접근으로 고차 단조 보존 방식을 불균일 망으로 확장하는 것은 쉽다.
식(10.325)의 장점은 우편의 첫항이 통상적인 1차 공여셀 근사를 준다는 것이다. 이와 같이 식(10.325)의 우편의 둘째항을 추가하여 생성되는 2차 근사가 쉽게 컴퓨터 프로그램에서 선택으로 될 수 있다.
단조를 확실케하기 위해 미분 A의 값을 제한하는것이 필요하다. 참고[VanLeer77]에 의하면 A의 값이 이 계산에서 사용되는 중앙Q미분의 최소 크기의 두배가 넘지 않도록 되어야 한다.
  (10.328)
더구나, Qi가 지역 최대 또는 최소값-즉, 식(10.328)에서 나타나는 두 개의 중앙 미분이 서로 다른 부호라면-이때 A는 0이 되고 공여 셀 근사가 사용된다.

Locally Implicit Approximation of Advection 이류의 지역적 내재적 근사

시간에 따라 진화하는 방정식 해를 구하는데 이용되는 외재적(즉, 시간에 따라 전진하는) 유한 차분근사는 시간증분의 크기에 제약을 받는다. 이런 형태의 근사에서는 모든 종속변수들의 값이 작은 시간 증분의 연속을 통하여 시간에 따라 전진한다. 이 때 시간 단계n까지의 모든 단계의 변수값이 알려져있으므로 이들이 변수들의 변화율을 추정하는데 이용될 수 있다. 이는 작은 시간δt 후인 tn+1 = tn + δt에 해당하는 단계n + 1에서의 종속 변수값들이 무엇인지를 추정하는 것을 가능하게 한다.

이러한 근사 형태에서의 작은δt의 필요성은 종속변수의 변화율이 보통 이 변수와 공간 바로 이웃의 값과 차이의 견지에서 평가되기 때문에 발생한다. 증가하는 시간 단계 크기에서 한 변수는 바로 인근들뿐만 아니라 더욱 멀리에 있는 값들에 의해서도 영향이 미쳐질 것으로 예상된다. 그러므로 안정과 정확성을 확실케하기 위해 시간 전진 크기δt에 제한이 있어야 한다.

>외재적 유한 차분 방정식은 사용하기에 편하지만 시간단계 크기를 제한해야 하는 필요성 때문에 오랜 계산시간이 소요될 수 있다. 이 제약은 극복하기 위해 문제를 일으키는 시간 단계 제약을 제거하기 위해 내재적 유한 차분법에 의존하는 것이 바람직할 수 있다. 기본개념은 진전된 시간단계n + 1 에서 변수의 값및 변화율의 근사를 포함하는 것이다. 지금 계산하려고 하는n+1수준 값이 같은 값에 의존하므로 차분방정식은 내재적이라고 불린다.

유체 이송방정식의 이류항을 근사하기 위한국소적 내재적 방법이FLOW-3D에 존재한다. 내재적 처리는 시간과 공간에 대해 선택적으로 적용될 수 있으므로 전반적으로 내재적방법 사용시의 간접 비용및 부정확성을 감소시킬 수 있다.단순한 1차원 스칼라 이류식을 고려해보면,

(10.329)

공여 이류(u > 0)에 기초한 상응하는 풍상(upwind) 유한 차분식은 다음과 같으며,

  (10.330)

여기서   (10.331)

시간 수준은 변화율Aj에 대한 근사에서 지정되지 않고 있다. S에 대한n단계값을 사용하는 외재적 근사는δtδx/u보다 작은 것을 필요로 하는 안정성 한계를 가진다.

내재성을 더하기 위해 Sj에 관한Taylor급수로 Aj를 확장한다. 1차 보다 큰 차수의 항들을 무시하면, 다음과 같은 식을 얻을 수 있다.

  (10.332)

식(10.331)으로부터
  (10.333)

지금 식 (10.330) 은 다음과 같이 표현될 수 있다.

  (10.334)

지금 이류 유속은 이에 내재적 정도를 더하는 를 포함한다. 이 근사의 안정성을 체크해 보자. 무조건적 안정성을 체크하는 점근적δt테스트는δt로 나누고 무한대로 보내는 것이다.

이 극한에서 첫 항은 사라지며 결과는S의 변화에 대한 새 방정식이며,

  (10.335)

이는 수정된 식이 무조건 안정적이라는 것을 뜻한다(∂A/∂S가 0이 아닌 한). 사실, 이 결과 식은 시간 단계만 ∂A/∂S의 역수로 치환된 것을 제외하고는 원래의 외재적 표현과 같아 보인다. 이  새 ”시간단계”는 정확히 외재적 안정성, δx/u을 위한 극한 시간단계이므로 식(10.335)의 근사는 필연적으로 안정적이다.

무한 시간단계의 제한을 취하지 않으면S의 전진을 위해 새로운 국소적으로 내재적이고 무조건적 안정적인 유한 차분방정식을 얻으며, 다음과 같이 쓰여질 수 있다.

  (10.336)

이 근사는 효과적으로 시간단계를 외재적 안정성을 위해 필요한 값보다 항상 작은 새로운 값으로 치환하기 때문에 모든 시간 단계크기에 대해 안정적이다. 식(10.336)은Sn+1 = Sn일 때 아직도 같은 정상상태, 즉A = 0, 해를 가진다. 이 방정식에서의 유속항은 이웃한 셀들의 시간n + 1의 값에 결합되지 않았으므로 이를 구하기 위한 반복법이 필요하지 않다.

방법은 무조건적 안정성을 가지는 이류 유속량들의 내재적 수치근사를 생성하기 위한 간단한 기법이다. 이류가 예기되는 상류의존을 유지하며 유속을 본질적으로 시간 단계당 하나의 격자 셀로 제한한다. 활성화되면 이 방식은 모든 유체 이송 방정식에 적용된다.차분근사는 단지 이웃 양들의 값을 수반하고 반복과정을 포함하지 않으므로 국소적으로 내재적 이류는 시간 단계 크기가 너무 커지면, 정확하지 않을 수 있다. 즉 제약이 없는 내재성은 부정확성에 이를 수 있다. 자유표면상의 유체 분율같은 이류 양의 큰 구배는 내재적방법 사용 시 커다란 수치 에러를 유발할 수 있다.
이러한 부정확성을 감소시키기 위해under-relaxation은 내재적으로 처리되지 않으면 계산상의 불안정성이 발생할 수 있는 위치에서 만으로 제한된다. 다른 말로 국소적으로 내재적인 근사는 단지 시간단계 크기가 셀 내의 속도와 셀의 크기에 의해 정의되는 지역 한정성 한계를 넘을 때와 넘는 장소에 있는 셀들에서만 적용된다. 모든 다른 시간 및 장소에서는 정상적인 보존 외재근사가 이용된다.자유표면 처리시의 부정확성을 방지하기 위해 내재적 이류 방식은 움직이는 자유표면을 포함하는, 즉 자유표면에 수직인 커다란 속도 성분을 가지는 셀들에는 적용되지 않는다.. 이 제한은 시간단계 크기가 이런 셀들에서는 통상적인 방식으로 제약되어야 한다는 것을 뜻한다.

Density Evaluation 밀도 평가

식(10.9)에서 필요한 밀도 는 다루려는 문제의 유형에 따라 평가된다. 한 물질의 비압축성 유동에서 모든 밀도값은 유체 1의 밀도로 지정된다. 2유체 비압축성 유동에서는 밀도는 체적분율F에 의거하여 두 유체 밀도(ρ1 와 ρ2)의 가중평균이다.

  (10.337)

유체밀도ρ1 와 ρ2는 기화 잔류물 절에서 논의된 바와 같이 온도에 의존할 수 있다.

1유체문제의 자유표면에서 유체는 상변화 모델이 사용될 때 기화하거나 응축될 수있다. 이는 경계면의 기체 쪽이 일정 압력을 가지거나 균일한 기포지역일 때의 두 경우에 해당한다. 유체가 하나 또는 더 많은 용질을 포함하면 용질들의 농도는 액체의 증감에 따라 바뀌어야한다. 일반적으로 용질은 액체의 기화시 더 농축된다.

응축이나 기화의 경우에 유체와 관련된 스칼라 농도는 상변화에 따라 변화된 농도를 가질 것이다. 표면요소에 유체가 반보다 적게 차있으면, 이 때 농도변화는 농축지역이 표면요소 두께의 반에 해당하는 정도로 표면요소의 주 이웃으로 퍼져갈 것이다.

기화가 충분히 발생하여 용질의 농도가 충분히 높아지면 표면에 막이 생기거나 액체가 완전히 증발하면 고체표면에 잔류물이 발생할 수 있다. 이를FLOW-3D에서 모사하기 위해 잔류 모델이 선택되어야 한다. 이 모델의 활성화는 용질이 일단 농도가 사용자가 지정한 최대 패킹밀도에 도달하면 움직이지 못하는 잔류물이 형성되게끔 한다. 한가지 이상의 용질이 존재하면 잔류모델은 잔류물의 원인이된 모든 용질 전체를 기록한다.

압축문제에서 밀도는 연속방정식(10.1)으로부터 평가되며, 아래와 같이 체적 가중 상태밀도 방정식에 해당한다.

  (10.338)

여기서ρ1 과 ρ2는 압력과 내부에너지의 함수이다.

밀도에 대한 압축유동 상태방정식 함수는 간단히 수정을 허용하게끔 프로그램되어 있다; 그러나 주어진 함수는 유연성을 갖는다. 두 유체 모두 일정한 비열(사용자가 지정할 수있는)을 가져서 비 내부에너지는 다음과 같다.

   (10.339)

유체 1은 비압축성으로 가정되며 이의 밀도는 RHO1으로 지정된다. 유체 2는 다음식과 같이 이상가스이며,

  (10.340)

여기서RF2는 사용자 지정 가스상수이다.

Wall Shear Stresses 전단응력

식(10.11)에서의 벽응력은 커다란 벽면적과 작은 유동공간을 가지는 셀들에서 발생할 수 있는 수치 불안정성을 피하기 위해 내재적으로 포함되어 있다. 예를들면w-속도 방정식에 대한 기본 접근은 다음과 같다. w에 영향을 미치는 벽전단이w를 둘러싼x y셀면에 위치한 벽면적으로부터 발생할 수 있다. 이런 면들 중 하나에서 유동 면적율A가 1보다 작으면 나머지 면적율(1 − A)이 응력이 생성되는 벽이라고 간주된다. 예를들면w우편의x면상에서 층류시의 벽전단에의한 가속wsz∂/∂x(µ(∂w/∂x))의 근사치이며,

   (10.341)

여기서Ax, δx µw가 위치한 셀내에서 평가된다. 면적율Azw가 위치한 같은 유한 체적면에서 평가된다. 속도w는 0이거나 또는 이동 요소나 망 경계에서의z방향 접선속도와 동등하다. w는 두 셀사이 경계에 있으므로Ax는 이 셀들에 대한 평균이다.

유사한 응력들이 4개의 주변 셀 벽의 각각에서 평가되고 이들의 합이 전체응력이 된다. 식(10.341)에서 나타나는w가 시간 단계n에서의 값이면 근사는 외재적이나 수치 불안정성이 있을 수 있다. 이 가능성을 피하기 위해w는 전체 모멘텀 방정식으로부터 계산된 새 값을 취한다. 이는w에 대한 모멘텀 방정식을 내재적으로 만들지만w에 대해 선형이므로 해를 구하는 것은 문제가 되지 않는다.  난류 전단 응력에서는 식(10.341)이 다음과 같이 바뀌며,

  (10.342)

여기서u+ = u*/u는 식(10.299) 으로부터 계산된다.

접선속도에 대한 2차식은 시간 단계 n + 1에서 취해진 우측의 괄호 안에서의w와 함께 선형화된다.물체표면으로부터의 질량주입에 따라 발생하는 유효 벽전단응력은 대체로 같은 방식으로 모델링된다. 유효응력은, 예를 들면,

   (10.342)

여기서QSRw가 정의된 경계의 한편에 있는 망 셀들에서의 평균질량 소스율이다. 소스가 유체가 지역속도를 갖기 때문에 음의값(싱크)이면 응력이 사용되지 않는다. 여기서 다시 새w속도가 이용되기 때문에 모멘텀방정식의 기여는 증강된 수치 안정성을 위해 내재적이다.

벽에서의 미끄러짐 형상은 단지 층류에서만 허용되고 마찰인자κ를 이용하여 모델된다. 미끄러짐을 위해 수정된 식(10.341)에의해 주어지는 벽전단응력은

 (10.344)

여기서ws는 벽에서의 유체 미끄럼 속도이다.

  (10.345)

벽표면에서의 난류및 층류 두 벽전단응력은 벽조도를 정의하여 수정될 수 있다. 조도는 길이의 차원을 갖는다. 어떤 의미에서 길이는 조도 요소들의 크기에 비례한다. 이는 분자점도에ρ, ξ, 그리고 w의 곱을 추가함으로써 통상적 전단응력 계산에 포함되어 있으며 여기서ξ는 표면조도이며w는 지역 유체속도와 벽 속도간의 차이이다.

이의 실행에서 층류에서의 벽전단응력은ρ(ν + ξw)w/δx와 같다. 난류에서는 ξ가 두 특정길이 규모 중에서 클 때 점도의 변화(즉, ν에서ν +ξw로)가ν/w에 의해 정의된 특정 길이 규모에 대한 로그 대수의존도를ξ로 자동적으로 전환하는 경우를 제외하고는 벽 법칙관계는 부드러운 벽에서와 같은 형태를 유지한다.

조도는 더 큰 값이 사용될 수 있지만 물체 경계에서 격자 셀크기 보다 작아야 한다.  마찰 조도변수ROUGH의 값은 물체/유체 열전달에 영향을 주지 않는다(또한Surface Area Evaluation를 참조).

Porous Baffle Flow Losses 다공 배플 유동 손실

배플에서의 유동손실에 대해 벽전단응력에 대해서와 유사한 내재적 해가 사용된다. 이 경우의 단지 차이는w (또는 u v)의 내재적 방정식이 2차이고 그 해의 제곱근을 취할 필요가 있을 수 있다는 것이다.

Generalized Drag Forces 일반 항력

다공질내 유동및 응고/용융 모델은 속도의 1차멱수, −Fdu, 에 비례하는 항력을 이용한다. 이런 모델들은 인위적으로 아주 큰 항력계수Fd,의 가능성을 가정하므로 이 경우에 모멘텀방정식은 항력과 압력의 균형을 이루게 될 수있다. 비압축성 유동에대한 이렇게 높은 항력의 한계에서 계산하기 위해 모멘텀방정식에서 뿐만 아니라 연속방정식에서도 항력항을 내재적으로 처리할 필요가 있다. 이는 항력에서시간 전진에서의 속도를 사용하고 대수적으로 새 속도에 대한 미분방정식을 해석함으로써 이루어진다. 이 결과는 새 속도의 모든 기여를 (1+Fddt)항으로 나눈 값이다. 모든 속도/압력 보정 동안에 이 추가 항의 효과를 유지하는 것은 압력구배와 항력사이의 균형이 달성되는 것을 확실케 해주며 연속방정식을 또한 만족시킨다.

Non-Inertial Reference Frame비관성 기준계

FLOW-3D는 움직이는 탱크내의 유동 해석을 위한 두 가지 계획이 있다: General Moving Objects Model 방법과 비관성계 방법이다. 전자는 관성계에 고정인 계산 망내에서 탱크를 이동시킨다(이 장과Model Reference에있는General Moving Objects Model 참조) 후자는 이 절에서 보여주는 탱크와 함께 움직이는 비관성계의 계산 망를 가진다.

비관성계에서는 기준계의 변환에 따라 유체가 관성계에 대한 비관성계의 움직임을 기술하는 일련의 독립변수에 연관될 수 있는 “관성” 가속도를 겪는다. FLOW-3D에서 이 변수들은 테카르트 좌표로 쓰인 3개의 선형 가속도 성분, 3개의 각속도 (시간당Radian) 성분 그리고3개의 각가속도(시간 제곱당Radian) 성분이다. 이 모든 9개의 변수들은 시간의 함수이다. 또한 각속도와 각가속도 성분들 간에 만족되어야 하는 세가지 관계식이 있다. 그러므로 실제로는 단지 6개의 독립적인 시간의 함수가 있을 수 있다.

추가로 어떤 형태의 운동 지정을 간단히 하기 위한 세개의 (시간종속)변수가 있다. 이들은 세 개의 “offset”벡터(RCX, RCY 그리고 RCZ)의 데카르트 성분이다. 이 성분들은 좌표 원점에서 회전이 기술되어야 하는 점까지의 벡터를 지정한다. 이는, 예를들면, “외력”을 겪지 않고 움직이는 물체의 회전을 기술하는데 유용한데 이 경우 질량 중심이 편리하게 “회전”점으로 간주될 수 있다.

유체가 겪는 국소적 ”관성”력의 평가는 아주 일반적인 세개의 서브루틴accxcl, accycl, 그리고 acczcl에 의해 되어진다. 일반적으로 이 루틴들은 사용자에 의한 수정이 필요하지 않다. 그러나 9개의 시간 종속변수들의 평가는 완전히 일반화될 수가 없다. 그러므로 이 평가는 하나의 서브루틴motion으로 분리 되어지며, 이는 필요시 사용자에 의해 수정될 수 있다. 또 다른 특별한 경우는 Coupled Rigid Body Dynamics 기술에서 논의되며 이 절에서는 더 이상 논의되지 않을 것이다.

See also: 또한 참조하라;

  • 비관성 기준계 표기
  • 비관성 기준계 모델을위한 강체역학 알고리즘
  • 비관성 기준계 운동 방정식
  • 비관성 기준계 강체역학
  • 비관성 기준계 응용예제: 원심주조
  • 중력
  • 비관성 기준계 의 충격운동
  • 비관성 기준계 운동
  • 표로 주어지는 부드러운 운동

Gravitational Force in Non-inertial Frames 비관성계내의 중력

많은 경우에 균일한 중력장내 기준계의 운동에 관심이 있을 수있다. 한 예는 해상의 유조선에서의 화물의 출렁임이다. 이런 해석을 용이하게 하기 위해 특정 “중력 벡터” 모델이 사용될 수 있다. 이 모델에서 사용자는 namelist MOTN에서 입력변수GRAVX, GRAVY, 그리고 GRAVZ를 통해 중력의 초기 성분들(데카르트 계산기준계에 상대적인)을 지정해야 한다. 이 때 루틴motion은 기준계가 회전할 때 유체 해석 알고리즘에서 요구되는 순간 성분들을 평가할 것이다. 기준계의 병진은 이 중력성분들에 영향이 없음을 주목한다.

Specifying Non-Inertial Reference Frame Motion비관성계 운동의 지정

어느계산에서 비관성 기준계모델을 이용하기 위해 우선 이 선택이 IACCF = 1을 지정하므로써 켜져야 한다. 이는 FLOW-3D가 MOTN namelist를 읽고 처리하여 계산시 적절한 시점에 motion, accxcl, accycl, 그리고 acczcl 서브루틴을 불러오게 한다.

기준계의 이동을 motion에 지정하기 위한 4가지의 선택이 있으므로 이들 중에서의 선택을 위해 새 입력 변수가MOTN에 추가되어 있다. 이 값들은 다음을 뜻한다:

IATYPE = 0 : x축에대한 단순 회전 또는 가속도 또는z축에대한 회전

IATYPE = 1 : 선형가속도와 각 속도/가속도의 조화진동

IATYPE = 2 :기준계 변수의 표를 통한 지정

다음 절은 이들 각각에 대해 상세히 기술한다.

Simple Rotation or Acceleration along the X-Axis  X축에 대한 단순회전및 가속도

이 선택은FLOW-3D의 이전 버젼에 호환성을 위해 주어진다. 같은 기준계운동이 원하면 다른 선택을 이용하여 이루어질 수 있다. 이 선택에 관련된 입력 변수들이다: RCX, RCY, RCZ, A0, OMG0, SPIN 그리고 RPS. 세가지 운동이 이 선택으로가능하고 원하면 결합될 수 있다:

데카르트 x축에 대한 가속도는 다음과 같다.

  (10.346)

상쇄좌표RCX, RCY 및 RCZ를통한 회전은 rad/초로 주어진다. x에평행한 축에대한 회전은 아래와 같이 주어지며,

  (10.347)

z에 평행한 축에대한 회전은 다음과 같다.

  (10.348)

Harmonic Oscillation 조화진동

이 선택은 선가속도와 각속도를 갖는 각 성분에 대한 독립적인 단순 조화함수을 이용을 가능케 한다. 이 선택으로 중력성분의 자동평가가 이루어질 수 있다. 이 선택에 관련된 입력변수들이다: RCX, RCY, RCZ, GRAVX, GRAVY, GRAVZ, TFREQX, TFREQY, TFREQZ, TMAGX, TMAGY, TMAGZ, TPHIX, TPHIY, TPHIZ, RFREQX, RFREQY, RFREQZ, RMAGX, RMAGY, RMAGZ, RPHIX, RPHIY 그리고 RPHIZ.

각 선형가속도의 성분은 다음과 같이 주어지며,

 (10.349)

여기서

  • f는 주기입력(TFREQX, 등등),
  • d는 변위입력(TMAGX, 등등),
  • φ는 위상각 입력(TPHIX 등)이다.

변위함수를 요구되는 선형 가속도로 변환하는 주파수의 제곱과 부호 변화에 주목하라. 각속도의 각성분은 다음으로 주어지며:

  (10.350)

  (10.351)

여기서 각속도에서와 같은 입력변수가 호환성을 확실히 위해 사용된다.

Tabular Specification 표를 통한 지정

계산시 서브루틴motion에 의해 읽혀질 데이터를 제공함으로써 6개의 기준계 변수에 대한 데이터를 표를 통해   지정할 수 있다. 기준계변수의 3개는 항상 선가속도의 데카르트 성분을 갖는다. 다른 3개의 변수는 각가속도(IATYPE = 2 or 4)의 성분이거나 각속도(IATYPE = 3 or 5)의 성분일 수 있다. 나머지 3개의 변수는 호환성을 맞추기 위해 적절히 수치 미분이나 적분을 이용하여 계산된다. 이 선택으로 중력성분은 자동적으로 평가된다.

관련 입력변수들이다: RCX, RCY, RCZ, GRAVX, GRAVY, GRAVZ, ACONV, 그리고 OMCONV.

IATYPE = 4 또는 5값은 충격 간격뿐만 아니라 순조롭게 변하는 가속도의 지정을 허용한다. 또한 같은 값은 원하면 특정한 시간에 사용될 시간 단계 크기를 지정하게끔 한다. IATYPE = 2 또는 3는 단지 순조롭게 변하는 운동만을 허용한다.

실행중 서브루틴motion 은 데이터 파일movin(이 파일의 이름은 사용자가 정의한다; 또한 이 위치로의 완전한 경로가 사용자에 의해 정의될 수 있다)으로부터 필요한 데이터를 읽으려고 할 것이다.

이 파일은 자유 포맷 구조를 갖는다. 처음 카드는motion에 의해 건너 뛰어야 할 줄의 수를 나타내는 정수를 가져야 한다. 이는 고유한 서식및 검사의 주석을 달게끔 해준다. 이 “두서(header)” 줄 뒤에 일련의 줄들이 있어야 한다. 각 줄은 특정 시간에서의 기준계 변수를 기술한다. 이 시간들은 엄격하게 증가하는 순서로 되어야 한다!

시간은 이 줄에서 첫 입력이어야 한다. 시간은 또한FLOW-3D내의 시간과 일치하여야 한다. 재시작 계산에서 시간은 0이아닌 선택된 값, TREST으로부터 시작한다. 다음 세 입력은 선형가속도의x-, y-, 그리고 z-성분이다.

이 줄에서의 다음 세 입력은 단위 시간당 radian인 각속도의 세x-, y-, 그리고 z성분 이다(IATYPE = 3 또는 5이면). 다음 세 입력은 단위시간 제곱당 radian인 각 가속도의 세x-, y-, 그리고 z성분이다(IATYPE = 2 또는 4이면).

 

IATYPE > 3의 경우에만 읽혀지는 다음 두 입력은 원하는 시간 단계 크기(이 값은 양수일 경우에만 사용된다)와 충격간격의 기간이다. 충격 간격 기간 변수는 순조로운 보간점들(양이아닌)과 충격구간 데이터(양수인) 사이를 구별한다.

서브루틴motion은 선형 가속도 값을 입력변수ACONV로 그리고 각의 변수들을OMCONV로 곱함으로써 입력파일 내의 데이터 단위를 전환할 것이다. 전환된 값은 솔버 요약 파일hd3out file에 인쇄된다. 값들은 이용 가능한 데이터 사이에서 보간된다. 맨 처음시간 값 전에 그리고 맨 마지막 시간값 후에는 일정한 값들이 사용되는데 이 경우 메시지는 데이터 파일값들이 외삽되고 있는 것을 보여주는hd3out에서 인쇄된다.

위에서 암시된대로, IATYPE > 3이면 충격및 순조롭게 변하는 가속도를 결합할 수 있다. 이는 어떤 카드 집단에서는 충격기간에 대해 양의 값 그리고 다른 카드 집단에서는 0 또는 음의 값을 지정함으로써 이루어진다. 서브루틴motion은 이런 종류의 데이터 파일을 우선 부드럽게 변하는 데이터 점의 변수들을 현재의 시간으로 보간하고 이 때 현재의 시간이 충격 안에 있다면 충격값을 추가한다. 이는 급격한 변화를 가지는 가속도를 가능케한다. 이는, 예를들면, 우주선에 탑재된 제어 시스템 작용에 대응하여 발생할 수 있다. 가속 충격에 마주치기 전에FLOW-3D가 시간단계 크기를 감소시킬 필요를 예상할 수 있도록 시간 단계크기 조절이 주어진다. 이는 또한 FLOW-3D가 선택된 시간 단계 크기가 충격 간격보다 크기 때문에 충격기간 동안에 뜻하지 않게 건너뛰는 것을 방지한다. 자동 시간단계 선택을 할 수 있게 되면 시간 단계크기는 충격 중에 큰값으로 회복되어 전체 계산 비용을 크게 감소시킨다.

Coupled Rigid Body Dynamics 결합 강체역학

비관성좌표계 운동에서 기술된 바와 같이FLOW-3D는 관성계에 대해 상대적으로 움직이는 좌표계에 탱크를 내장함으로써 탱크 내의 유체 유동을 탱크를 연구하는데 사용되어질 수 있다. 어떤 경우에는FLOW-3D가 탱크 내부 유체에 의해 작용하는 힘과 모멘트의 영향을 포함하는 결합된 방식으로 이 기준계의 운동을 또한 해석할 수 있다. 이는 “Coupled RigidBody Dynamics” 모델로 알려져 있고 이절에서 기술된다.

이모델의 기본 접근은 유체지역을 관성공간에서 움직일 수 있는 더 큰 강체의 부분으로 간주하는 것이다(Rigid Body Dynamics참조). 이 물체는 다양한 힘과 토크를 받으며 잘 알려진 방식에 따라 이들에 반응할 것이다. 계산망은 강체내에 내장되어 있으므로 이의 움직임은 강체의 움직임과 같다.

물체에 작용할 수있는 한 세트의 중요한 힘과 토크는 이에 인접한 유체의 움직임에 의해 형성된 것이다. 대여섯 가지의 다른 결합방식이 입력에 의해 선택될 수 있고 하기에 간단히 기술되어 있다.  많은 다양한 힘들이 물체에 또한 작용할 수 있다. FLOW-3D는 이 힘들을 세 가지 부류로 나눈다:

  • 환경적 힘과 토크
  • 힘과 토크 조절
  • 중력가속도에 의한 힘

처음 두 종류는 문제의 특정해석을 위해 프로그램될 수 있는 두 개의 서브루틴과 연결해줌으로써 처리된다. 환경과 조절 부류의 차이는 단지 편리상이며 임의로 해석자에 의해 결정될 수 있다. 점 중력체로부터의 중력의 효과는 좌표계와 물체의 상대운동의 효과를 포함하여 자동적으로 계산된다. 이 모델은 다른 것들 중에서도, 항공업계의 사용이 의도 되었으므로  중력의 방향과 크기는 계산을 통하여 조절된다.

어떤 경우들에서는 하나 이상의 유체 탱크가 강체의 운동에 영향을 줄 지도 모른다. 이런 탱크들은 다른 물성(밀도, 점도 등)을 지닌 유체를 포함할지도 모른다. 다중 탱크 모델이 이런 경우의 연구를 위해FLOW-3D에 포함되어 있다.

Input Data 입력데이터

결합된 강체역학 모델은XPUT namelist에서IACCF = 2로 함으로써 활성화 된다. 또한 강체의 건조질량 및 관성 모멘트를 지정하는 것이 필요하다; 이들은namelist RBDATA에서 주어진다. 다른 입력데이터는 물체의 초기위치, 방향 그리고 운동과 중력체의 변수들을 포함한다.

입력변수IRBACC는 사용될 결합 알고리즘과 선형및 각 가속도의 초기값이 계산되는 방식을 조절한다. 유체 반응력이 초기에는 알려져 있지 않으므로 가속도의 초기 평가가 필요하다. IRBACC는 처음 유체 시간 단계의 시작시에 해석자가 유체를 무시하거나 정지된 질량으로 취급하도록 하게 해준다.

일반적으로 결합모델의 사용은 조절및 환경의 힘과 토크의 지정루틴인, RBCTRL 과 RBENVR의 특수 버젼을 필요로 한다. 배포된FLOW-3D버젼은 다양한 입력변수를 받아 들이는 특수 목적의 조절 루틴을 포함한다. 이 표준 루틴에 3가지선택이 가능하다. 이들의 선택은 다음과 같이IATYPE를 지정함으로써 가능하다.. :

IATYPE = 0: 힘이나 토크가 없음

IATYPE = 1: 스프링과 감쇠 모델로부터 계산되는 조절력

IATYPE = 2:스프링과 감쇠모델로부터 계산되는 조절토크

IATYPE = 3: 조절력과 토크성분의 표 형태  지정

IATYPE = 1 나2에 대한 스프링과 감쇠 계수는namelist MOTN에서RADKRB, RDMPRB 또는 THEKRB, TDMPRB로 지정된다. IATYPE = 1에대해 또한REQRB (평형각은IATYPE = 2에 대해서는 항상 0이다)를 통해 평형 반경을 지정할 수 있다.

표형태 입력선택은 기준계운동의 표형태 지정과 유사하다. 계산시 서브루틴rbctrl에 의해 읽혀질 수 있는 데이터 파일을 제공함으로써6개의 조절력과 토크의 데카르트 성분들을 위한 표형태로 지정할 수 있다. 이 파일의 디폴트 이름은 thrust.inp이다.

전환 상수(ACONV 와 OMCONV)가 추력 데이터 파일과 함께 사용될 수 있다. 파일로부터 읽혀진 힘의 성분은RBCTRL에서 사용하기 위해ACONV에 의해 곱해지며 한편 토크성분은OMCONV에 의해 곱해진다.

이 파일은 자유 포맷 구조를 갖는다. 처음 카드는motion에 의해 건너 뛰어야 할 줄의 수를 나타내는 정수를 가져야 한다. 이는 고유한 서식및 검사의 주석을 달게끔 해준다. 이 “두서(header)” 줄 뒤에 일련의 줄들이 따라야 한다. 각 줄은 특정 시간에서의 힘과 토크성분을 기술한다. 이 시간들은 엄격하게 증가하는 순서로 되어야 한다!

 

시간은 이 줄에서 첫 입력이어야 한다. 시간은FLOW-3D내의 시간과 일치하여야 한다. 재시작 계산에서 시간은 0이 아닌 선택된 값, TREST으로부터 시작한다. 다음 세 입력은 조절력의 성분들이다(x, y,z순서로). 그 다음 세 입력은 같은 순서로 조절 토크의 성분들이다. 다음 입력은 원하는 시간 단계의 크기이다. 이 값은 단지 양수이면 사용된다. 마지막 입력은 이 행이 부드럽게 변하는 힘과 토크 변화에 대한 보간점으로 사용되는 지 또는 충격 운동 지정으로써 사용되는 지를 결정한다. IATYPE > 3 이면 이 입력은 항상 이 행의 마지막에 있어야 한다. 이 데이터 입력이 양수이면 이는 충격 기간을 표시한다. 그렇지 않으면 이 행은 부드럽게 변하는 함수상의 보간점으로 간주된다.

위에서 암시된대로, 충격및 순조롭게 변하는 가속도를 결합할 수 있다. 이는 어떤 행들에서 충격기간에 대해 양의 값들 그리고 다른 행들에서는 0 또는 음의 값들을 지정함으로써 이루어진다. rbctrl은 이런 종류의 데이터 파일을 우선 부드럽게 변하는 데이터 점의 변수들을 현재의 시간으로 보간하고 현재의 시간이 충격 안에 있다면 충격값을 추가한다. 이는 급격한 변화를 가지는 가속도를 허용한다. 이는, 예를들면, 우주선에 탑재된 제어 시스템 작용의 결과로 발생할 수 있다. 시간 단계크기 조절이FLOW-3D가 충격가속도가 발생하기 전에 시간단계 크기를 감소시킬 필요를 예상할 수 있도록 주어질수있다. 이는 또한 FLOW-3D가 선택된 시간 단계 크기가 충격 간격보다 크기때문에 충격기간을 뜻하지 않게 건너뛰는 것을 방지한다. 자동 시간단계 선택이 가능하게 되면 시간 단계크기는 충격 중에 큰 값으로 회복되어 전체 계산 비용을 크게 감소시킨다.

Explicit Solution Algorithm and Stability Limitations 외재적해석알고리즘및 안정성 한계

외재적 결합 해법은  (1) 기준계의 이동에대한 전단계 해를 이용하여 유체유동방정식을 해석하고, (2) 이에 따른 유체의 힘과 토크를 평가하며, (3)강체 운동방정식을 표현하기위해 조절및 환경 효과를 추가하고, 그리고 (4)새 가속 변수들, 질량중심 위치 그리고 강체 방향성(자세)에 대한 방정식들을 해석 함으로써 진행된다. 새 기준계 변수들이 유체 해를 위한 다음 사이클에서 이용된다. 이 결합의 외재적 형태는 유체의 질량모멘트나 관성모멘트가 강체의 것보다 크면 불안정해 진다고 알려져 있다. (이는 단지 개략적 관계식이다. 유체의 질량모멘트나 관성모멘트가 다소 클 경우에도 불안정성이 발생하지 않는 경우들도 있다. 반대인 경우도 있을 것이라고 여겨진다.)

Implicit Solution Algorithms 내재적 해석알고리즘

두가지 내재적 해석 알고리즘이 존재한다. 첫째(IRBACC = 2)는 유체질량이 단계4에서 고정된 것으로 간주되는 것을 제외하고는 외재적방식과 유사하다. 즉 유체 힘은 마치 유체가 계산망에 고정된 것같이 기준계 가속도에 반응한다.

이의 변경은 외재적 해의 안정에 대한 제약을 없앤다. 이 알고리즘의 정확도는 “고정”의 근사가 해에 적절한 정도에 달려있다.

 

이 방식에 대한 개선은 계산된 기준계운동을 마지막 사이클의 결과의 추정으로 간주하는 것이다. 이때 유체 해는 추정된 기준계운동을 이용하여 같은 사이클에 대해 반복된다. 새로운 유체의 힘은 기준계운동을 재평가하기 위해 환경및 조절력과 결합된다. 가속도의 수렴기준이 만족될 때까지 반복이 계속된다. 이는 반복적 내재적 알고리즘이다(IRBACC = 3). 또 이의 정확도는 “고정” 근사의 유효성에 달려있다.

입력변수및 출력변수의 완전한 기술을 포함한 결합된 강체역학 모델의 추가정보를 위해References [Sic92], [Sic95]를 참고하라.

[FLOW-3D 이론] Numerical Approximations 수치근사

Numerical Approximations 수치근사

Overview 개요

FLOW-3D는 유한차분 (또는 유한 체적)근사를 이용하여 이전 절들에서 기술된 방정식들을 수치적으로 해석한다. 유동지역은 고정 직각 셀들로 이루어진 망으로 세분된다. 각 셀에서 모든 종속변수들의 지역 평균값들이 있다. 다음에서 설명되듯이 모든 변수들은 셀-면(엇갈린 망 배열)에 위치한 속도를 제외하고 셀의 중심에 위치한다.

굴곡진 물체, 벽경계 또는 다른 기하학적 형상은 유동에 열려있는 셀들의 면적및 체적율의 정의에 의해 망 내에 들어있다( FAVORTM 방식 [HS85]).

지배방정식에 대해 이산 수치근사를 구성하기 위해 유한체적들이 각 종속변수 위치를 둘러쌈으로써 정의된다. 각 유한 체적에 대해 표면유속, 표면응력 그리고 체적력이 주변 변수값들에 의해 계산될 수 있다. 이 양들은 운동방정식에 의해 표현된 보존법칙을 근사를 구성하기 위해 결합된다.

방정식의 대부분 항들은 다양한 내재적 선택이 또한 가능하지만 지역변수들의 현재시간 단계에서의 값들을 사용하여 즉 외재적으로 평가된다. 이는 대부분의 목적에서 단순하고 효율적인 계산방식을 제공하지만 계산적으로 안정되고 정확한 결과를 유지하기 위해 제한된 시간단계의 사용이 필요하다.

이 외재적 공식에의 한 중요한 예외는 압력에 의한 힘의 처리에 있다. 압력과 속도는 운동방정식에서의   시간-전진된 압력과 질량(연속)방정식에서의 시간-전진된 속도를 사용함으로써 내재적으로 결합되어 있다.  유한차분 방정식의 이런 내재적 공식은 저속의 압축 및 비압축성 유동문제에서 효과적이고 안정된 해석으로 이끈다.

내재적 압력-속도 공식은 반복법에 의해 풀어져야 하는 결합된 일련의 방정식들이된다. FLOW-3D에서는 3가지의 이런 기법이 주어진다. 가장 간단한 것은successive over-relaxation (SOR)방식이다. 더 내재적 방법들이 요구되는 경우에서는special alternating-direction, lineimplicit method (SADI)이 사용 가능하다. 뒤에 기술되겠지만 SADI기법은 해석할 문제의 특성에 따라 한,둘 또는 세 방향 모두에서 사용될 수 있다. 마지막으로Generalized Minimal Residual (GMRES) 방식이 또한 가능하다. GMRES솔버는 특히 커다란 소스항이 존재하거나 변하는 셀크기를 가지는 망에서SOR이나SADI방법보다 우월한 수렴 성질을 가지고 있다. GMRES솔버는 또한   메모리 공유 병열화에 더 적합하다. 이는 FLOW-3D의 기본(디폴트) 솔버이다.

FLOW-3D 에서 사용되는 기본 수치방식은 시간과 공간증분에 대해 공식 1차 정확도를 가진다. 유한 차분 망이 불균일할 때는 이 정확도를 유지하기 위해 특별한 주의가 필요하다. 2차 정확도의 선택 또한 가능하다. 어느 경우에나 경계조건은 모든 상황에서 최소한 1차정확도를 가진다. 예를들면 물체에 의해 부분적으로 점유된 셀에서FAVORTM방식은 셀내 경계조건의 1차보간과 대등하다. 그러나 유체/물체 경계면에서의 열전달 경계조건의 실행은 셀 크기에 대한 2차 정확도를 가진다. 다음 절에서 이 수치 기법은FAVORTM를 통한 유한 차분과 유한 체적근사의 특정한 예를 통해 더 정확하게 된다.

[FLOW-3D 이론] Auxiliary Model/Electro-mechanics 전기역학

Electric Field Model

To simulate physical processes such as movement of charged mass particles, particle and liquid dielectrophoresis, and electro-osmosis, an electric field distribution is needed. In FLOW-3D , the electric potential is solved for using the following equation.

(75)\nabla \cdot \left( {K\nabla \phi } \right) = - \frac{\rho _e}{\varepsilon _0}

with the electric field calculated by

(76){\mathbf{E}} = - \nabla \phi

where \rho _e, \varepsilon _0 are free charge density (i.e., electric charge per unit volume) and permittivity of vacuum or void (defined by elperm in FLOW-3D ) respectively while K is the spatially-varying dielectric constant.

Numerical solution of the Poisson equation Eq. (75) is done by an iterative solver using the GMRES method.

At open boundaries (i.e., boundaries through which fluid flow can occur), either an insulation condition (i.e., {\mathbf{n}} \cdot \nabla \phi = 0 where \textbf n is a normal vector on the boundary in question) or a specified value of the electric field can be specified. Solid objects can be assigned time-dependent potentials if they are conductors. Alternatively, an object can be a dielectric material with an assigned dielectric constant. In this case, the electric potential is computed within the object. Also, solid dielectrics may have non-zero conductivities that support free charges and currents.

Optionally, a charge density equation that includes charge convection, charge relaxation and charge sources associated with non-uniform electric properties is solved simultaneously with the electric potential. Please refer to Electric Fields in the Model Reference chapter and the Technical Notes listed at the end of this section for more information.

Most people use the SI system of units for electrostatic and electromagnetic problems. In the SI system the standard unit of charge is the coulomb and the unit of potential is the volt. Electric field intensities are measured in Newton/coulomb in the SI system. When using the CGS system of units, it is customary to express charge and potential in the “electrostatic” units of statvolt and statcoulomb, where:

1 coulomb (C) = 2.998 \times 109 statcoulomb (statC)

1 volt (V) = 3.336 \times 10-3 statvolt (statV)

Electric forces are then,

1 Newton = 1 C \cdot V \cdot m-1 = 105 dyne

1 dyne = 1 statC \cdot statV \cdot cm-1 = 10-5 Newton

If the fluid or solid regions are conducting then charges may develop in response to the applied field and to changes in conductivity and dielectric properties. If a thermal model is activated then Joule heating due to these currents is computed in all conducting fluids and solids. The heat source per unit open cell volume and time associated with passing electric current through conducting fluid, Q_{Joule}, is

(77)Q_{Joule} = F \cdot VF \cdot \kappa|\mathbf{E}|^{2}

where F and VF are the fluid and open volume fractions in cell and \kappa is the fluid electrical conductivity. For two-fluid problems, F is replaced with 1.0 and \kappa with the volume-averaged mixture value of the electrical conductivity. Similarly, the Joule heat source in solid is

(78)Q_{Joule solid} = (1 - VF) \cdot \kappa_{solid}|\mathbf{E}|^{2}

More details are given in Flow Science Technical Notes #52, #56, #69 and #70 on the users site at http://users.flow3d.com/technical-notes/.

 

Electro-osmosis

Generally, most substances such as silica and glass will acquire a surface electric charge when brought into contact with an aqueous (polar) medium (electrolyte solution), as shown in the figure below.

A layer called the electric double layer (EDL) is formed close to the charged surface in which there is an excess of counter-ions over co-ions to neutralize the surface charge. There are more counter-ions than co-ions in the region near the fluid/solid interface. The electric potential created due to the EDL is called the \xi -potential and is assumed to be imposed at the solid surface. The \xi -potential is a property of the solid-liquid pair and can be measured experimentally. The thickness of the EDL is indicated by the following parameter, called the Debye shielding distance or Debye length.

(79){\lambda _D} = {\left( {\frac{{\varepsilon RT}}{{2{F^2}{z^2}{c_0}}}} \right)^{\frac{1}{2}}}

where:

  • \varepsilon = liquid permittivity
  • R = gas constant
  • T = temperature
  • F = Faraday’s constant
  • c_0 = ion concentration
  • z = valence

More details on EDL and related physics can be found in [Pro94].

Electro-osmotic flow or electro-osmosis refers to the fluid motion that occurs when an electric field is applied to an electrolyte solution in the vicinity of a charged surface. The process can be described by the following equations

(80){\mathbf{F}} = {\rho _E}{\mathbf{E}}

(81){\rho _E} = - 2F{c_0}\sinh (\frac{{F\psi }}{{RT}})

(82){\nabla ^2}\psi = \frac{{\rho E}}{\varepsilon }

where \psi is the \xi -potential. An insulation-like boundary condition is imposed for the \xi -potential on all mesh boundaries with a symmetry condition. On solid (obstacle) surfaces, the \xi -potential is imposed. Here Boltzmann charge density distribution for liquid with single valence is assumed. But the user can easily provide other different charge density distributions in two simple functions named psi and dfdpfi released to users. How to do this customization is described in these two functions.

Please refer to Electro-osmosis in Model Reference to get related information on how to use the model.

 

Particle Movement and Fluid Flow Due to Electric Field

An electric charge or charge dipoles can be carried by molecules, small droplets, and particles which are called charge, or charge dipole carriers, or mass particles in FLOW-3D . If the charge of a mass particle, i, is e_i, the electric force (called coulomb force) acting on this particle is

(83){{\mathbf{f}}_i} = {e_i}{\mathbf{E}}

where {\mathbf{E}} is the electric field intensity and its calculation is described above. If a mass particle, i, carries charge dipoles and the corresponding dipole moment is p_i, the electric force (called polarization force) imposed on this mass particle is

(84){{\mathbf{f}}_i} = ({{\mathbf{p}}_i} \cdot \nabla ){\mathbf{E}}

If all energy losses in the carriers are neglected, the dipole moment is calculated by

(85){{\mathbf{p}}_i} = 4\pi {\varepsilon _0}r_p^3K_f\left( {\frac{{{K_p} - K_f}}{{{K_p} + 2K_f}}} \right){\mathbf{E}}

where K_p is the particle dielectric constant, K_f is the fluid dielectric constant, r_p is the particle radius and e_0 is the permittivity of vacuum. Then, the electric force imposed on mass particle, i, carrying the charge dipoles can be cast into

(86){{\mathbf{f}}_i} = 2\pi {\varepsilon _0}r_p^3K_f\left( {\frac{{{K_p} - K_f}}{{{K_p} + 2K_f}}} \right)\nabla {E^2}

Instead of paying attention to the carriers mentioned above, we consider fluid where these charges or charge dipole carriers are distributed. According to Newton’s law, the fluid will experience some body forces due to existence of charge or charge dipoles. The body force due to free charge is

(87){\mathbf{F}} = - {\rho _e}{\mathbf{E}}

where \rho _e is the free charge density. The body force due to charge dipoles is

(88){\mathbf{F}} = ({\mathbf{P}} \cdot \nabla ){\mathbf{E}}

where {\mathbf{P}} is the density of the dipole moment. For dilute dipoles in fluid, this density can be calculated by

(89){\mathbf{P}} = {\varepsilon _0}(K_f - 1){\mathbf{E}}

The body force is then calculated by

(90){\mathbf{F}} = \frac{1}{2}{\varepsilon _0}(K_f - 1)\nabla {E^2}

In FLOW-3D , users can provide a charge density distribution. This distributed charge density is defined using the FLOW-3D scalar variable having index IECHRG (a FLOW-3D input parameter).

Movement of particles carrying induced charge dipoles due to polarization is called particle Dielectrophoresis (DEP) while fluid flow due to existence of charge dipoles from polarization is called liquid Dielectrophoresis.

With the aforementioned equations included in FLOW-3D , users can simulate both particle and liquid Dielectrophoresis. Please refer to Dielectrophoresis to get related information on how to use the models.

[FLOW-3D 이론] Auxiliary Model/Electro-mechanics 전기역학

Electric Field Model 전기장 모델

전하를가진 질량입자의 운동, 입자 및 액체 유전이동 그리고 전기삼투 같은 물리적 과정을 모사하기위해 전장의 분포가 필요하다. FLOW-3D 에서 전위는 다음식을 이용하여 해석된다.

(75)\nabla \cdot \left( {K\nabla \phi } \right) = - \frac{\rho _e}{\varepsilon _0}

전기장는 다음에 의해 계산된다.

(76){\mathbf{E}} = - \nabla \phi

여기서 ρe, ε0 는 각기 자유전하밀도(즉, 단위 체적당 전하)와 진공이나 공간내의 유전율(FLOW-3D 에서 elperm으로 정의되는)이며 K 는 공간적으로 변화하는 유전상수이다.

Poisson식(10.75) 의 수치해석은 GMRES 방법을 사용하는 반복법 솔버에 의해 이루어진다.

개방된 경계 (즉, 유동이 발생할 수 있는 경계)에서 절연조건(즉, 이며 n 은 문제의 경계에서의 법선 벡터)이나 특정 전장값이 지정될 수 있다. 고체는 전도체이면 시간의존 포텐셜이 지정될 수 있다. 그렇지 않으면 물체는 지정된 유전장수를 가지는 유전물질일 수 있다. 이 경우 전위가 물체내부에서 계산된다. 또한 고체유전체는 자유전하와 전류를 가능케 하는0이 아닌 전도도를 가질 수 있다.

선택할 수 있는 가능성으로, 불균일한 전기 물성과 관련된 전하 대류, 전하 완화 그리고 전하 소스를 포함하는 전하밀도방정식이 동시에 전위와 함께 해석된다. 더 상세한 내용을 위해 모델 레퍼런스장에 있는 Electric Fields와 이 절 마지막 목록에 있는 Technical Notes 를 참조하라.

대부분의 사람들은 정전기 및 전자장 문제에 SI계 단위를 사용한다. SI 계에서 전하의 표준단위는 coulomb 이고 전위 단위는 volt이다. 전장 강도는 SI 계에서 Newton/coulomb으로 측정된다. CGS 단위 사용시 통상적으로 전하와 전위를 “ 정전기” 단위인 statvolt 와 statcoulomb로 표현하며, 여기서:

1 coulomb (C) = 2.998 \times 109 statcoulomb (statC)

1 volt (V) = 3.336 \times 10-3 statvolt (statV)

전기력은,

1 Newton = 1 C \cdot V \cdot m-1 = 105 dyne

1 dyne = 1 statC \cdot statV \cdot cm-1 = 10-5 Newton

유체나 고체지역이 전도성이면 전하가 가해진 장 및 전도도 및 유전 물성의 변화에 반응하여 전하가 발달할 수있다. 열모델이 활성화되면 이 전류에 따른 Joule 가열이 모든 전도성 유체 및 고체 내에서 계산된다. 전도성 유체를 통해 흐르는 전류에 의한 단위 개방 셀 체적 및 시간당 열소스 QJoule, 는

(77)Q_{Joule} = F \cdot VF \cdot \kappa|\mathbf{E}|^{2}

여기서 F V F 는 유체와 셀내의 개방 체적율이고 κ 는 유체 전기 전도도이다. 2유체 문제에서 F 는 1.0 으로 대체되고 κ 는 전기전도도의 체적평균 혼합값을 가진다. 유사하게, 고체내의 Joule 열소스는

(78)Q_{Joule solid} = (1 - VF) \cdot \kappa_{solid}|\mathbf{E}|^{2}

더 상세한 내용은 사용자 사이트 http://users.flow3d.com/tech-notes/default.asp 에 있는 Flow   Science  Technical Notes   #52,   #56,   #69, #70에 주어져 있다.

 

Electro-osmosis 전기삼투

일반적으로 실리카나 유리같은 대부분의 물질은 밑에서 보여진 바와 같이 수성(극성의) 매질(전해질)과 접촉할 때 표면전하를 가지게 될 것이다.

Surface electric charge when brought into contact with an aqueous (polar) medium (electrolyte solution)

전기 이중층(EDL) 이라고 불리는 층이 전하면에 가까이 형성되는데 이 표면전하를 중성화시키기 위해 같은 이온보다 더 많은 반대 이온이 과도하게 존재하게된다. 유체/고체면 가까이 지역에 같은 이온보다 더 많은 반대이온이 존재한다. EDL 에 의해 생성된 전장는 ξ−전위라고 불리며 고체표면에 부과된다고 가정된다. ξ−전위는 고체/액체 쌍의 물성이며 실험적으로 측정될 수 있다. EDL의 두께는Debye shielding distance 또는 Debye length라고 불리는 다음 변수에 의해 나타내진다.

(79){\lambda _D} = {\left( {\frac{{\varepsilon RT}}{{2{F^2}{z^2}{c_0}}}} \right)^{\frac{1}{2}}}

여기에서:

  • ε 는 액체유전율
  • R 은 가스상수
  • T 는 온도
  • F 는 Faraday 상수
  • c0는 이온 농도
  • z 는 원자가

EDL 과 관련 물리에 관한 상세내용은 [Pro94]에서 찾을 수 있다.

전기-삼투 유동 또는 전기-삼투는 전장이 전하를 띤 표면 근처에서 전해질 용액에 가해졌을 때 발생하는 유동을 뜻한다. 이 과정은 다음 식으로 기술될 수 있고

(80){\mathbf{F}} = {\rho _E}{\mathbf{E}}

(81){\rho _E} = - 2F{c_0}\sinh (\frac{{F\psi }}{{RT}})

(82){\nabla ^2}\psi = \frac{{\rho E}}{\varepsilon }

여기서 ψξ−전위이다. 절연 같은 경계조건이 대칭 조건을 가지는 모든 격자 경계상에ξ−전위를위해 가해진다. 고체(물체) 표면에서 ξ−전위가 가해진다. 여기서 단일 원자가의 액체에 대한 Boltzmann 전하밀도 분포가 가정된다. 그러나 사용자는 사용자에게 배포된 psi dfdpfi 두 간단한 함수로 쉽게 다른 전하밀도 분포를 줄 수가 있다. 어떻게 이 사용자 주문 변환을 하는지는 이 두 함수에 기술되어 있다.

이 모델을 어떻게 사용하는지에 대한 정보를 위해 모델 레퍼런스내 Electro-osmosis를 참조하라

 

Particle Movement and Fluid Flow Due to Electric Field 전장에 의한 입자운동 및 유체유동

전하 또는 전하 쌍극자는 FLOW-3D에서 전하, 또는 전하 쌍극자 운반자, 또는 질량입자라고 불리는 분자, 작은 방울, 또는 입자에 의해 운반될 수 있다. 질량입자 i 의 전하가 ei이면 이 입자에 작용하는 전기력(coulomb 힘이라고 불리는)은

(83){{\mathbf{f}}_i} = {e_i}{\mathbf{E}}

여기서 E 는 전장강도이고 이의 계산은 위에 기술된다. 질량입자 i 가 전하쌍극자를 운반하고 상응하는 쌍극자 모멘트가 pi라면 이 질량에 가해지는 전기력(편광력이라 불리는)은

(84){{\mathbf{f}}_i} = ({{\mathbf{p}}_i} \cdot \nabla ){\mathbf{E}}

이다.

운반자들에서의 모든 에너지 손실이 무시되면 쌍극자모멘트는 다음에 의해 계산되며

(85){{\mathbf{p}}_i} = 4\pi {\varepsilon _0}r_p^3K_f\left( {\frac{{{K_p} - K_f}}{{{K_p} + 2K_f}}} \right){\mathbf{E}}

여기서 Kp 는 입자 유전상수, Kf 는 유체 유전상수, rp 는 입자반경 그리고 e0 진공에서의 유전율이다. 이때 전하 쌍극자를 운반하는 질량입자 i 에 가해진 전기력는 다음으로 나타나게 된다.

(86){{\mathbf{f}}_i} = 2\pi {\varepsilon _0}r_p^3K_f\left( {\frac{{{K_p} - K_f}}{{{K_p} + 2K_f}}} \right)\nabla {E^2}

위에서 언급된 운반자에 집중하는 대신에 이 전하나 전하 쌍극자가 분포되어 있는 유체를 고려한다. 뉴튼 법칙에 의하면 유체는 전하나 전하 쌍극자의 존재로 인한 체적력을 받는다. 자유전하에 의한 체적력은

(87){\mathbf{F}} = - {\rho _e}{\mathbf{E}}

여기서 ρe 는 자유 전하밀도이다. 전하 쌍극자에 의한 체적력은

(88){\mathbf{F}} = ({\mathbf{P}} \cdot \nabla ){\mathbf{E}}

여기서 P 는 쌍극자모멘트의 밀도이다. 유체 내의 희석된 쌍극자에 대해 이 밀도는 다음에 의해 계산된다

(89){\mathbf{P}} = {\varepsilon _0}(K_f - 1){\mathbf{E}}

체적력은 이때 다음과 같이 계산된다.

(90){\mathbf{F}} = \frac{1}{2}{\varepsilon _0}(K_f - 1)\nabla {E^2}

FLOW-3D 에서 사용자는 전하밀도분포를 줄 수가 있다. 이 분포된 전하밀도는 색인 IECHRG(한 FLOW-3D입력변수)를 가지는 FLOW-3D 스칼라 변수를 이용하여 정의된다.

극성화에 의한 유도된 전하 쌍극자를 운반하는 입자의 운동은 입자 유전이동(DEP) 이라고 불리며 극성화에 의한 전하 쌍극자 때문에 발생하는유동은 액체 유전이동이라고 불린다.

FLOW-3D 에 포함된 앞서 언급된 방정식으로 사용자는 입자 및 액체 유전이동을 모사할 수 있다. 이 모델 사용 관련정보를 얻기 위해 Dielectrophoresis 를 참조하라.

[FLOW-3D 이론] Auxiliary Model/Drift-Flux Model 드리프트-플럭스 모델

다중 성분을 가지는, 예를 들면 유체/입자, 유체/기포, 그리고 성분들의 밀도가 다른 유체/유체 혼합물로 구성되있는 유체에서 성분들은 각기 다른 유동속도를 가지는 것으로 관측된다. 속도차이는 밀도차이가 불균일한 체적력을 초래하기 때문에 발생한다. 가끔 속도의 차이는 아주 뚜렷할 수도 있는데, 예를 들면 공기중에서 떨어지는 빗방울들, 또는 물에서 가라앉는 자갈 등의 경우다. 그러나 많은 조건에서 상대 속도는 한 성분의 다른 성분에 대한 “드리프트 “로 기술될 만큼 충분히 작다. 예를 들면 공기 중의 먼지와 물속의 토사이다.

”드리프트(drift)” 의 구별은 연속성분 내에서 이동하는 분산된 성분의 관성이 중요하냐와 관련이 있다. 상대속도의 관성이 무시되고 상대속도가 성분들 간의 구동적 힘(즉, 중력이나 압력구배)과 이에 반하는 항력간의 균형으로 축소된다면 이때 우리는 “표류-플럭스” 근사에 대해 이야기할 수 있다. 표류속도는 일차적으로 질량과 에너지의 이송에 영향을 미친다. 약간의 모멘텀도 또한 이송되지만 이 역시 매우 작고 FLOW-3D 의 표류 모델에서 무시되어 있다.

드리프트 모델의 발상은 성분 사이의 상대운동은 이산요소(즉 입자같은)보다는 연속체로 근사될 수 있다는 것이다. 이는 이산 요소의 운동 및 상호작용의 추적을 계산할 필요가 없으므로 계산효율를 증강시켜 준다.

FLOW-3D 에서 드리프트 속도가 사용될 수 있는 4가지 물리적 상황이 있다.

  • 혼합물이 각기 밀도가 ρ1 ρ2인 성분으로 되어있는 1유체 밀도변화 유동
  • 액체와 고체 혼합물이 각기 밀도가 ρ1 ρ2인 응고가 일어나는1 유체
  • 밀도 ρ1 ρ2를 가지는 2 비압축성 유체
  • 비압축성 성분을 가지는 압축성기체. 이 경우 압축성기체의 밀도는 상태방정식에 의해 주어지며 비압축성 물질은 밀도 ρ1를 가지고 항상 이는 기체밀도 보다 훨씬 크다고 가정된다.

표류 근사에서 상대속도의 공식화는 다음과 같이 진행된다. 유동이 2개의 이산 요소 또는 상인 상황에서, 하나는 연속상이고 다른 하나는 이산상이며, 이는 불연속적이고 연속인 상에의해 둘러싸여 있다고 가정한다. 이와 똑같은 유체계가 반대의 구성을 가질 수 있다: 많은 양의 디젤유를 오염시키는 작은 양의 물 같은 경우에 물이 이산 상이고; 역으로, 소량의 디젤유가 많은 물속에 펴져있다면 디젤유가 이산 상이다.

두 성분 유체의 비압축성 유동에 대해  \nabla \cdot \textbf u = 0 인 u = f1u1 + f2u2 를 정의한다. 혼합물을 구성하는 두 성분의 체적율은 f1 f2로 표시되며, 여기서

두 상이 비압축이라고 가정하면 연속상의 모멘텀균형은

(63){f_1} + {f_2} = 1

이며. 이산 상에 대해서는

(64)\frac{\partial {\mathbf{u}}_1}{\partial t} + {{\mathbf{u}}_1} \cdot \nabla {\mathbf{u}}_1 = - \frac{1}{\rho _1}\nabla P + {\mathbf{F}} + \frac{K}{f{\rho _1}}{\mathbf{u}}_r

반면에 dispersed phase 에 대해서는

(65)\frac{\partial {\mathbf{u}}_2}{\partial t} + {\mathbf{u}}_2 \cdot \nabla {\mathbf{u}}_2 = - \frac{1}{\rho _2}\nabla P + \mathbf{F} - \frac{K}{(1 - f){\rho _2}}{\mathbf{u}}_r

여기에서:

  • u1 와u2 는 각기 연속상과 이산상의 미세속도를 나타내고,
  • f 는 연속상의 체적율이다.
  • 미세속도는 작지만 유한한 체적의 유체에 대한 각 상의 속도를 뜻한다.
  • F 는 체적력이고,
  • K 는 두 상간의 상호작용을 연관시키는 항력계수이며,
  • ur 은 이산상과 연속상 사이의 상대속도 차이이다.

(66){\mathbf{u}}_r = {\mathbf{u}}_2 - {\mathbf{u}}_1

드리프트-플럭스 모델의 목적은 체적 가중속도 u¯에 대한 두상 의 운동을 계산하는 것이다. 체적가중 평균속도는

(67){\mathbf{\bar u}} = f{\mathbf{u}}_1 + (1 - f){\mathbf{u}}_2

이고, 체적-가중 평균속도가 질량-가중 평균속도 대신에 선택되는데, 이는 질량연속이 어떤 보정없이 자동적으로 만족되기 때문이다.

(68)\nabla \cdot {\mathbf{\bar u}} = 0

모멘텀 보존이 만족되듯이.

식 (10.64) 을 식 (10.65) 으로부터 빼면 상대속도의 식을 얻게되고 여기서 K 는 단위 체적당 항력계수이다.

(69)\frac{{\partial {\mathbf{u}}_r}}{\partial t} + {\mathbf{u}}_2 \cdot \nabla {\mathbf{u}}_2 - {\mathbf{u}}_1 \cdot \nabla {\mathbf{u}}_1 = \left( {\frac{1}{\rho _1} - \frac{1}{\rho _2}} \right)\nabla P - \left( {\frac{1}{{(1 - f){\rho _2}}} + \frac{1}{{f{\rho _1}}}} \right)K{{\mathbf{u}}_r}

목적은 상대속도 ur 을 결정하는 것이다. 식(10.69)을 그대로 사용하면 두 성분 유동을 위한 두 속도장을 구성할 것이다. 그러나 우리는 단순히 표류-플럭스 모델 근사를 했다. 즉 상대속도가 거의 정상상태이고 이류항들은 서로 상쇄한다(즉, 작은 상대속도 ur 에 대해)고 가정하는 것이다. 이 가정 하에 우리는 다음식을 얻는다

(70)\left( {\frac{1}{\rho _1} - \frac{1}{\rho _2}} \right)\nabla P = \left( {\frac{{f{\rho _1} + (1 - f){\rho _2}}}{{f(1 - f){\rho _1}{\rho _2}}}} \right)K{{\mathbf{u}}_r}

상대속도 ur 은 각 상의 미세속도에 의거하므로 이때 항력은 부유상의 체적율에 대한 어떤 정보를 가져야한다. 예를들면, 아주 작은 양의 부유상을 가지는 부유물은 성분 간에 아주 작은 모멘텀의 교환을 가져올 것이다.

부유상이 같은 크기인 입자들로 구성되어있고 단위체적당 n 개가있다고 추정하면, 이때에

(71){{\mathbf{u}}_r} = \left( {\frac{{{V_p}}}{{{K_p}}}} \right)\frac{{f({\rho _2} - {\rho _1})}}{{\bar \rho }}\nabla P

여기에서:

  • Vp = (1 − f)/n 는 입자 한개의 체적이며Kp 는 연속유체 내에 속력 |ur| 로 움직이는 단일 입자에 대한 항력계수이고,

(72)\bar \rho = f{\rho _1} + (1 - f){\rho _2}

위는 체적가중 밀도이다.

드리프트 계수와 상대 유동속도 간의 2차 의존도가 모델에서 사용된다. U 가 연속 유체 내에서 움직이는 입자의 상대속도의 크기라면, 이때

(73){K_p} = \frac{1}{2}{A_p}{\rho _c}\left( {{C_d}U + \frac{{12{\mu _c}}}{{{\rho _c}{R_p}}}} \right)

여기에서:

  • Cd 는 항력계수,
  • Rp 는 평균 입자 반경이며
  • Ap 는 구로 가정되는 입자의 단면적이다.

식(10.67)에서 U가 나타나있기 때문에 반복적 해법이 상대속도U를 구하기 위해 필요하나 반복법은 상당히 효율적이어서 많은 시간이 소요되지 않는다.

분산된 물질들의 체적율이 충분히 작지 않을 경우 성분간의 모멘텀 교환 계산을 위한 단일 입자 항력의 사용은 그렇게 적합하지는 않다. 이들 입자/입자 상호작용을 고려하기 위해 가장 흔히 사용되는 수정은 Richardson-Zaki 상관관계라고 불리는 실험적으로 결정된 관계식이다.

Richardson-Zaki상관관계는 입자 Reynolds 수, Re = (2ρcRpU/µc)에 의존한다. 이 상관식은 표류 속도를 계산된 속도 ur에 분산 성분 체적율에 지수 ζ 를 취한 값으로 대체한다.

(74){\mathbf{u}}_r^{eff} = {\mathbf{u}}_r\cdot max(0.5,f)^\zeta

지수 ζ 는 Richardson-Zaki 계수 rzmltζ0의 곱이며, 즉, ζ = rzmlt * ζ0 여기서

Re < 0.2 {\zeta _0} = 4.35
0.2 < Re < 1.0 {\zeta _0} = 4.35/{Re}^{0.03}
1.0 < Re < 500 {\zeta _0} = 4.45/{Re}^{0.1}
500 < Re {\zeta _0} = 2.39

조정된 또는 유효값은 모든 드리프트 플럭스를 계산하는데 이용된다. 이 상관식은 단지 2차 입자 저항 모델에만 사용된다.

[FLOW-3D 이론] Auxiliary Model/Defect Tracking 결함추적

Auxiliary Model/Defect Tracking 결함추적

주물의 기계적 물성은 산화막, 주조 충진 시 금속안에 갇힐 수 있는 기타 이물질과 공기와 관련된 금속의 “청결도“에 상당히 의존할 수 있다.

이 절에서는 표면난류에 기인하는 주조결함의 질적 예측을 하는 결함의 원천 및 추적 기법이 기술된다[BH98].
추적 모델은 또한 폼이 분해된 후 금속표면에 남아있는 잔류물로부터 발생하는 결함의 가능성을 예측하기 위해 로스트폼 모델과 결합되어 있다.

자유표면에서와 로스트폼 잔류로부터의 내포물들은 별도의 양들로 나타내진다.
자유표면 결함모델에서 양은 일정한 Defect generation rate DFTRF 율로 자유표면 상에 축적된다.
폼 잔류는 Residue generation rate DFTFOB 로 정의된 비례계수를 가지는 분해된 폼의 질량에 비례하는 양만큼 증가된다. 이송방정식은 식 (10.270) 에 주어진 것과 유사하게 수치적으로 이산 및 확산 항을 가지는 각 양들에 대해 해석된다.
공간에 대해2차적이고 단조성을 보존하는 기법이 증가된 정확도를 위해 이산항을 근사하는데 이용된다.
내포물의 확산은 보통 적어서 디폴트로 꺼져 있으나 Defect tracking Lost foam 모델을 위해 Molecular diffusion coefficient 와/또는 역Schmidt 수인 Turbulent diffusion coefficient multiplier 를 양의 값으로 정의함으로써 포함될 수 있다.

두 양 모두 오염 물질의 농도를 나타낸다.
오염으로 인한 결함의 가능성은 농도에 비례한다. 결함은 이 중 하나의 농도가 높게 나오는 지역에서 더 발생하기 쉽다.

지금 현재 부력이나 산화막의 강도 효과는 이 모델에 포함되어 있지 않다. 또한 몰드 벽에 접착하려는 산화막과 다공사형을 통해 배출된 폼 잔류물 같은 현상도 포함되어 있지 않다. 포스트폼 주조과정에서의 폼 분해에 의한 액체 생성물의 심지효과는 포함될 수 있다(자세한 내용은 Lost Foam Casting Process 참조).

모사가 끝날 때 내포물 농도의 공간적 분포는 결함이 있을 것 같은 위치를 알려준다.
대부분의 경우에 대다수의 내포물은 예측 되듯이 마지막으로 채워질 장소에서 발견된다. 다른 오염물은 금속전면에 모일 것이고 갈 곳이 없을 때까지 밀릴 것이다.
그러나 관심이 있는 곳은 두 전면이 만나서 금속내부의 표면 내포물이 갇히게되는, 또는 유체 재순환 지역이 발생하는 구석진 곳이다.
심지어 몰드 충진 후에 금속 내의 순환유동은 표면에서 생성된 오염물을 재 분포시킬 수 있다.

실제 결함이 있는 이런 형태의 모델의 예측을 상관시키시 위한 자세한 실험 비교없이 이 스칼라의 절대값에 어떤 의미를 두는 것은 불가능하다.
그러므로 반응 상수는 임의의 양수를 가질 수 있다. 그러나 심지어 가장 간단한 실험으로부터의 정성적결과도 주조물 내의 결함생성과 이들의 최후 분포에 대해 영향을미치는 과정에 대한 상당히 유용한 통찰을 준다.

[FLOW-3D 이론] Auxiliary Model/Core Gas Generation and Flow in Sand Cores and Moulds 모래 코어 및 사형 내 코어가스생성 및 유동

가끔 모래안의 모래 점결제의 열분해에의해 생성된 가스로부터 발생하는 가스 결함은 사형 주조 부품에서 발견된다. 계산 모델은 코어 내가스의 생성과 이송의 예측을 가능케 한다. 가스는 불량한 코어 배출 상황및 커다란 온도구배에 따른 코어지역을 통한 이송등을 고려하기위해 압축성으로 간주된다. 이 모델은 실제 몰드형상 및 코어 배출위치를 감안하며 임의의 코어 표면 상위치로부터 용탕내로 유입되는 가스의 양을 예측하는데 이용될 수 있다.

코어는 프로그램 내에서 다공이 없는 고체물질로 간주되므로 용탕이 내부로 유입될 수 없다. 코어가스 모델은 코어 물체를 마치 가스가 생성되어 하기와 같이 기술된 특별 경계 조건을 사용하는 경계에서 유출입할 수있는 다공질 물질로 간주한다.

생성된 코어가스는 이상기체로 간주되고 기체상수 Rcg 를 가지는 고정된 구성을 갖는다. 가스상수는 고정 체적 기구로부터 모아져서 가스압력이 측정되는 실험으로부터 추론될 수 있다. 측정된 전체 얻어진 표준체적 Vstd 과 점결제의 초기질량, mb 로부터 가스상수는 다음과 같다.

(53){R_{\rm{cg}}} = \frac{p_{\rm{std}}}{T_{\rm{std}}}\frac{V_{\rm{std}}}{m_b}

코어 가스의 미시적 속도, \vec u_{cg}, 는 다공질 유동 방정식에 의해 지배된다:

(54){\phi }{\vec u_{\rm{cg}}} = - \frac{K}{{\mu}_{\rm{cg}} }\nabla {p_{\rm{cg}}}

여기서

  • φ 는 코어 다공도
  • K 는 모래 투과도
  • µcg 는 코어가스 점도이며
  • pcg 는 코어가스 압력

거시적 관성항은 일반적인 코어가스 유동에서 작다고 알려져 있으므로 식 (10.54)에 포함되어 있지않다. 투과도 K 는 다음에 따라 항력계수 Fd 로 표현되어 있다.

(55){K = \frac{{\phi } {\mu }_{\rm{cg }} } {{\rho }_{\rm{cg}} F_ {\rm{d}}}}

그리고 항력계수 Fd 는 1차 Darcy 항과 2차 Forchheimer 항을 갖는다:

(56){F_ {\rm{d}} } = a\frac{{{\mu _{\rm{cg}}}}}{{{\rho _{\rm{cg}}}}}\frac{{{{(1 - \phi )}^2}}}{{{\phi ^2}}} + b\frac{{1 - \phi }}{\phi }{\vec u_{\rm{cg}}},

여기서:

  • \rho _{\rm{cg}} 는 코어가스 밀도이며
  • a, b 는 모래 특정의 선형 및2차 유동 손실계수이다.

코어가스의 밀도는 질량 이송방정식에 의해 지배된다:

(57)\frac{{\partial {\rho _{\rm{cg}}}}}{{\partial t}} + \nabla \cdot \left( {{\rho _{\rm{cg}}}{{\vec u}_{\rm{cg}}}} \right) = - \frac{{d{\rho _b}}}{{dt}}

여기서 ρcg ρb 는 미시적 코어가스 밀도 및 거시적 코어 점결제 밀도이다. 코어가스는 압축성이므로 심지어 가스 소스가 없더라도 코어가 가열될 때 코어내의 초기 공기의 열팽창 및 유동이 있을 수 있다.

코어가스 밀도는 더구나 이상기체 법칙에 의해 지배된다:

(58){p_{\rm{cg}}} = {R_{\rm{cg}}}{\rho _{\rm{cg}}}T

여기서 가스의 온도는 지역 코어의 온도와 같다고 가정된다. 이는 가스에 비해 고체 코어재질의 더 많은 열량때문에 합당한 근사이다.

코어와 가스의 온도가 같고 코어의 온도 T 는 이미 열전달 해석에 의해 계산이 되었으므로 이상기체 법칙이 식(10.57) 에 의해 정의된 밀도를 가지는 가스압력을 계산하는데 이용될 수 있다. 고체 점결제의 가스로의 전환은 Arrhenius 식에 의해 기술되며,

(59)\frac{{d{\rho _b}}}{{dt}} = - {\rho _b}{C_b}\exp \left( { - \frac{{{E_b}}}{{RT}}} \right),

여기서:

  • ρb 는 고체점결제밀도
  • Cb 는 경험적으로 얻어지는 반응율상수
  • Eb 는 성분결합에너지이며
  • R 은 일반가스상수이고
  • T 는 온도이다.

가스소스는 계산요소의 이산화크기에 따른 비물리적인 진동 거동을 나타낼 수 있다. 이를 완화하기 위해 단순한 세부분할 기법이 각 계산셀 내에서 이용된다. 고체점결제 밀도는 분할셀 교점에서 저장되며 분할셀 온도는 주 망의 교점으로부터 온도장의 선형 보간에 의해 얻어진다.

코어가스 모델은 몰드 충진모델과 동시에 사용되기 때문에 코어내 가스유동이 계산시간 단계 크기(외재적 이류 근사를위한Courant안정조건에 의해) 를 충진에 필요한 크기보다 작게 제한할 만큼 빠를 수가 있기 때문에 더 긴 계산시간을 초래한다. 이런 것에 대비하여 코어가스 모델에는 보조-시간-단계 기법이 포함되어있다. 이는 모델이 코어외부 표면의 경계조건을 통하여 금속충진과 결합되어 있기 때문에 가능하다(하기참조). 코어내 가스유동이 계산시간 단계크기가 충진의 크기보다 작다고 판명되면 코어가스계산의 시간단계크기는 하나의 안정적인 값으로 감소되고 모델은 금속 내에 고정된 열과 유동조건으로 대여섯 개의 분할시간 단계에서 해석된다. 분할시간 단계의 크기는 주 시간 단계 크기를 정수로 나눈 분수이므로 가스계산은 가스해를 옳은 시간까지 전진시키기위해 이 정수만큼의 시간 단계들에 대해 반복된다.

코어재질 경계에서 코어가스의 유출입이 있을 수 있다. 이 교환은 코어가스 모델의 경계조건으로 간주된다. 코어경계를 통한 가스의 통로는 코어 외부에 무엇이 있느냐에 달려있다. 예를 들어, 코어표면이 공기에 노출되어 있으면 가스는 압력 차이에 따르는 방향으로 경계를 통해 흐를 수 있다.

코어표면에 액체 금속이 있을 때 코어의 압력이 액체의 압력보다 크면 가스가 나올 수 있으나 금속은 코어 내로 들어갈 수 없다. 금속이 이미 코어표면에서 응고하면 가스는 그 표면을 통과할 수가 없다.

다른 경계조건이 즉, 코어표면이 그 몰드의 또다른 고체부분과 접촉할 경우 같은 코어 프린트 표면에서 발생한다. 이런 표면 위치에서 (가스)배출이 될 수 있도록 몰드 내에 채널이 만들어지지 않는 한 가스는 정상적으로 흐르지 않는다. 코어가스모델은 프린트표면에서 배출할 수있는 선택을 할 수가 있다.