[FLOW-3D 이론] Numerical Approximations 수치근사 – Diffusion Process, Heat Conduction and Heat Transfer 확산과정, 열전도및 열전달

Numerical Approximations 수치근사

Diffusion Process, Heat Conduction and Heat Transfer  확산과정, 열전도및 열전달

Diffusion Process  확산과정

내부 점성 전단(즉, 고체 경계로부터 떨어진 곳에서의 전단응력)과 유체 분율의 확산, 유체 밀도, 유체 에너지밀도, 난류에너지 그리고 난류 소산들은 모두 유사한 방식으로FLOW-3D에서 처리된다. 각 과정에서의 확산계수는 동적 점도µ에 비례하며 이는 FLOW-3D에 포함된 난류 모델과 사용자에 의해 포함된 다른 모델의 결과로 공간및 시간에 대해 변할 수 있다.  이 변화를 고려하기 위해 모든 확산 과정들은 더 통상적인 Laplacian 형태 보다 유속의 발산으로 공식화 된다.

확산과정을 위해 사용된 외재적 차분 기법은 상당히 간단하다. 점도와 확산된 양은 전시간 단계의 배열로부터 평가된다. 이는Stability Considerations에서 토론된 안정성 요구로 연결된다. 이러한 제약은 점성모멘텀 확산항의 내재적 처리가 선호되는 어떤 경우들에는 상당히 심각하다(Implicit Viscous Stress참조). 스칼라양 확산의 공간 차분 또한 상당히 단순하다. 망 셀 표면을 통과하는 유속량은 확산계수와 셀 표면적율의 평균곱에 면의 양 편에 있는 셀 중심 양들의 1차 차분을 곱하여 평가된다.

점성 응력의 공간 차분은 필연적으로 벡터 성격으로 인해 더 복잡하다. 간단히, 지역전단율을 얻기 위해 1차차분을 속도성분에 적용하고 점도와 셀 표면적율의 평균 곱을 곱한다. 그후 이 결과는 순수 점성응력을 근사화하도록 차분화된다.  모든 양들은 이 계산에서 외재적으로 평가된다.

Heat Conduction and Heat Transfer  열전도와 열전달

열전도및 열전달 항들은 유체 에너지 밀도 방정식(10.21)및 고체 온도 방정식(10.38) 둘 다에서 나타난다. 이항들은FLOW-3D에의해 유사하게 처리되고 관련된 시간단계 안정성 제약을 없애기 위해 선택할 수 있는 두가지 내재적 공식을 포함한다. 이 근사는 이 열과정과 관련된 시간 규모와 유체유동 자체, 특히 정상및 준정상 상태에서의 시간규모의 큰 차이 때문에 포함되어있다(Implicit Viscous Stress를 보라).

사용된 방법을보여주기위해 식(10.21)의 재공식화된 상사형

  (10.369)

을 고려하며, 여기서

  • TDIF는 유체내 열전도를 나타내고,
  • HADT는 유체와 인근 고체와의 열전달을 나타내며,
  • X는 식(10.21)의 나머지 항들을 보여준다.

지금 식(10.369)을 근사하기 위해FLOW-3D에 의해 사용된 차분관계를 고려한다. 통상적으로 에너지 유속량에 대해서  전진 시간차분및 1차공간차분을 이용하면, 아래와 같이 되며,

  (10.370)

여기서,

  • T는 유체온도를나타내고,
  • TW는 고체온도,
  • h 는 열전달계수,
  • WA는 경계면 면적,
  • A는 셀면 면적,
  • k는 평균열전도계수이며,
  • δ는 적절한 공간 증분이다.

식(10.370) 오른쪽 모든 항은 외재적 근사가 열유속량에 대해 사용되는 것을 뜻하는 “전” 시간 단계tn 에서 평가된다. 외재적 알고리즘은 시간-단계 크기 안정성 제약을 뜻한다(Stability Considerations참조). 외재적 근사는 같은 차분이 어느 셀 표면의 양쪽에서 사용되기 때문에 에너지를 보존한다.

Implicit Viscous Stress and Heat Flux Approximations  내재적 점성 응력및 열유속량근사

점성및 열확산에 대한 외재적 근사의 수치 안정성을 위해 필요한 최소 시간-단계 크기는 전체 망에 대해 추정되므로 이 시간-단계 크기 제약은 가끔 너무 심해서 효율적인 수치해석을 할 수가 없는게 확실하다, 이는 비록 변수들이 이 지역에서 변하지 않을 수있더라도 시간-단계 크기가 가장 작은 셀에서의 안정성 한계에 의해 조절되는 아주 찌그러진 망에서의 경우일 수도 있다.
두 개의 내재적 솔버가 이 문제를 완화시키기 위해 개발되었다. 첫번째는 모멘텀 방정식에서의 내재적 항들을 위한 Jacobi 반복법이다. 둘째는ADI방식에 근거한다.

이 방식들을 설명하기 위해 앞에서의 모멘텀및 연속 방정식(10.315), (10.316), (10.317) 및 (10.352)들로부터 면적및 체적분율 생략하고 2차원 데카르트 형상을 가정함으로써 단순화 할 것이다. 일정한 점도를 가지는 비압축성 유체에 대해 유동방정식은 다음과 같이 균일(단순화하기 위해)한 격자에서 차분화된다.

   (10.371)

    (10.372)

여기서,

  • i j는x 와 y방향에서의 셀 색인이고,
  • δx δy셀 간격이며,
  • FUX, FUY, FVX, 와 FVY는 Momentum Equation Approximations 에서 기술된 모멘텀 이류항들이다.
  • 벽 전단 응력은계수 wx wy,를 가지는   속도에 대해 선형으로 가정되고 fx fy 유체에 작용하는 모든 다른 힘을 포함한다.

상첨자n을가진 항들은 시간tn에서의 해를 이용하여 명시적으로 계산된다. 상첨자n + 1는 내재적으로 시간tn+1 = tn + δt에서 평가되는 항들을 표기한다.  식(10.371) 와 (10.372)은 서로 결합되어 있고 이웃 셀들에서의 변수로 연결되어 있으므로 한 반복법이 모든0 < i < Nx + 1, 0 < j < Ny + 1에 대해ui,j, vi,j pi,j해를 얻기 위해 이용되어야 한다.

두 반복 순환이Jacobi반복법에서 정의된다: 모멘텀 방적식에서의 점성항을 위한 하나와 연속방정식의 압력을 위한 다른 하나이다. 두번째 반복 순환은Pressure Solution Algorithm에서 기술된 것과 유사하다. 완전한 압력 수렴은 점성 반복 후에야 이루어진다. 다른 말로 각 점성 반복은 일반적으로 많은 압력반복을 필요로 한다(주어진 시간단계에서 모든 점성반복에 대한 평균압력 반복수는 진단 출력에서 계산된다).각 점성반복kv + 1에서 모멘텀방정식은 대각 내재적항ui,j vi,j들에 대해 해석된다(kv = 1에 대해서는 전시간 단계로부터의 해석이 이용된다):

  (10.373)

여기서u* v*는 식(10.371)의 모든 외재적 항들을 포함한다.
식(10.373)으로부터kv + 1단계의 속도를 식(10.372)으로 치환하면 압력방정식을 얻게된다. Pressure Solution Algorithm에서 기술된 표준SOR 방식과 유사하게 각 유한 체적내에서 발산이0이 되는 조건을 만족시키도록 압력이 조절된다.

  (10.374)
여기서,

  • S는 식(10.372)의 좌측이고,
  • ω는 입력변수OMEGA에의해 정의되는 완화인자이며,
  • kp는 압력 반복수이다.

압력과 속도성분들의 상응하는 조절은Dx Dy가 식에 존재한다는 점에서 식(10.356)에서 주어진 것들과는 다르다:
  (10.375)

  (10.376)

주어진 점성 반복kv에대해 압력 반복이 왼료되면 다음 점성 반복kv + 1 에 대해 수렴된 해가 식(10.373)을 풀도록 사용된다.
해는 점성과 압력반복 둘다가 수렴될 때 수렴이된다. 압력반복EPSI을위한 수렴 기준은 외재적 점성 알고리즘에 대한 경우와 같다. EPSI는 일반적으로 시간-단계의 크기의 함수이다. EPSI의 추가 승수EPSADJ가 수렴을 더 작게 하거나(EPSADJ < 1.0) 느슨하게(EPSADJ > 1.0)하는 데 이용될 수 있다. EPSADJ의 디폴트 값은 1.0이다.

연속적 감속완화(SUR) 또는 Jacobi 점성 반복은 연속 반복 사이의 정규화된 속도의 변화가0.01× EPSI를 넘지 않을 때 수렴되었다고 여겨진다. 추가 승수EPSVIS는 수렴을 더 완화하거나 작게하는 데 이용될  수 있다. EPSVIS의 디폴트 값은 1.0이다.
감속 완화인자 OMEGVS는 점성 반복의 수렴을 증진시키기 위해 식 (10.373)에서 사용된다. SUR반복법에 대한OMEGVS의 디폴트 값은 0.9이다. 이 방법은IMPVIS=1일 때 활성화된다.

선 내재적(ADI)점성 솔버에서 모멘텀과 압력방정식은 같은 반복내에서 함께 해석된다; 압력과 속도 반복간의 구별이 없다. 반복k+1에서 전의 반복k로부터 속도를 이용하여δp를 평가하기 위해 격자를    한번 통과한다:

  (10.377)
여기서 분모는 식(10.376)으로부터 평가되고 Dx Dy는 식(10.373)에서 정의된다. 식(10.377)의ω는 1.0으로 디폴트된다. 1보다 큰 값을 사용하는 것은 권장되지 않는다.

그때 모멘텀방정식은 3중 대각행렬 솔버를 이용하여 식(10.377)으로부터의 가장 최근 압력을 이용하여 해석된다. 우선, y-미분을 평가하기 위해k– th단계의 속도성분을 이용하여 x-방향에서 한 통과가 이루어진다:

  (10.378)

식(10.378)은 상첨자k + 1를 가지는 항들에 대해 풀어진다.  x-방향 미분을평가하기위해.  x-방향 통과 후에 얻어진 속도성분들을 이용하여y-방향으로 한 통과를 한다. 이때z-방향에서의 통과가 3차원 문제에서 더해진다.  또한 가속완화가 이 방법에서 사용되나 수치적 잡음을 최소화하기 위해 OMEGVS는 현재 1.0으로 디폴트가 된다.

반복과정을 위한 첫 추측으로 시간단계n에서의 속도들이 이용된다. 해법은 연속방정식의 잔여가 격자 내 모든 곳에서EPSI보다 작아지고 정규화된 모멘텀방정식 잔여가1.0-8 × EPSVIS보다 작을 때 수렴된다. EPSVIS는 디폴트값 1.0을 가지는 입력변수이다. 모멘텀방정식 잔여(x-방향에서)는 와  식(10.378)의 우측항의 계수 중의 최대값에 의해 정규화된다.

ADI반복법은IMPVIS = 2를 지정함으로써 활성화된다. 점성 저항이 우세한 힘일 경우SUR방법보다 더 빨리 수렴하고 더 정확한 해를 준다. 이는 모사가 점성 안정 한계보다   몇 차 정도 큰 시간 단계에서 진행할 때 반응고 금속, 폴리머 그리고 비결정 물질같은 아주 높은 점도 유동에 적용된다.

SUR방식은 압력에 의해 구동되는 유동에서ADI솔버보다 더 효율적인데 이 경우 점성응력이 압력에 부수적이고 시간단계 크기가 단지 점성 안정 한계보다 한 두차 정도 크기가 크다.

ADI반복법은 유체및 물체상에서 동시에 해석되는 열에너지 방정식의 해를 얻는데 이용된다. 이 반복 알고리즘에서 이용되는 상향 완화 과정은 수렴율을 증가시킨다. 상향 완화의 정도는 최적 수렴을 위해 지정된 디폴트 값을 가지는 OMEGHT에 의해 조절된다. 온도에 의해 정규화된 에너지방정식의 최대 잔여가EPSIHT = 1.0e-04 × EPSHTC보다 작을때 반복은 수렴되었다고 간주된다. EPSHTC의 디폴트값은 1.0이다.

열및 점성 해의 수렴기준은 자동적으로 계산 중에 지정되며 시간 단계마다 변할 수 있다.
난류 모델들은 내재적으로 처리되지 않는 다른 종류의 확산과정을 포함하므로 내재적 점성 선택은 난류모델들과 함께 사용될 수 없다.

내재적해석 방법은 열전도및 전달 과정에 의해 δt에 강제되는 제약을 없앨 수 있다. 그러나 다음의 언급은 내재적및 외재적 근사의 차이를 인식하는 데 유용할 수 있다. 수치적으로 안정된 일시적 외재적 해를 확실케 하는 최대 시간-단계 크기는 열파동이 공간에서 한 개 이상의 셀을 전파하지 않는 시간 주기로 설명될 수 있다. 더 큰 시간-단계가 내재적 알고리즘에서 사용될 때 이 파동들은 한 시간 단계에 더 멀리까지 전파될 수가있다. 이 기간 동안에 대 여섯개의 열파동들이 유동내 같은 위치를 통과하므로 이 위치에서의 유체 에너지및 온도는 이 같은 기간 동안에 복잡한 과도적 거동을 겪을 수 있다. 그러나 에너지방정식에서의 순간적 편미분은 한시간-단계에서 실질적으로 선형 에너지변화를 가정하는 1차 차분기법으로 근사화 되어있다. 그러므로 내재적기법에서의 큰 시간-단계를 사용하는 것은 해를 평준화(또는 확산화)하고 순간해의 정확도의 확실한 손실로 나타난다. FLOW-3D에서 내재적 해를 위한 시간 단계의 크기는 반드시 내재적 해의 정확성 및 효율을 유지하도록 자동적으로 조절된다. 시간 단계의 크기는 반복수IDTHT를 넘지 않도록 조절된다(디폴트는 20이다).

[FLOW-3D 이론] Numerical Approximations 수치근사 – Scalar Advection 스칼라 이류

Numerical Approximations 수치근사

Scalar Advection 스칼라 이류

압축성 연속방정식(10.1), 유체분율 방정식(10.19), 내부에너지 방정식(10.21), 그리고 난류에너지 와 소산 방정식 (10.270) 와(10.275) 이 모든 식들은 벡터형태에서 공통적 형태를 지니는데,

  (10.362)

이는 확산과 소스 기여를 포함할 수 있는 우측 항RHS를 가지는 양Φ에 대한 스칼라 이류 방정식이다. 다섯개 모든 방정식은 망의 셀중심에 위치한 같은 이류양Φ를 가지는 유한 차분 알고리즘에 의해 표현된다. 밑에 기술된 대로 경계면을 뚜렷하게 하기위한 알고리즘이 특정 유체 형태의 양들에 대한 유속량 항∇ · (uAΦ)에 대해 이용된다.

밀도, 에너지, 유체분율 그리고 난류 이송방정식과 별도로FLOW-3D는 임의의 스칼라 양들에 대해 식(10.362)와 유사한 추가의 보존방정식을 가진다.이 양들은 유동에서의 피동적 오염물의 진화를 기술하거나 이송방정식을 필요로 하는 새 물리적 모델, 예를들면, 화학 반응 과 시간의존 유변학모델을 개발하는데 이용될 수 있다.

각기 다른 유체와 관련된 모든 이류 양들은 이류된 유체의 양δF에 비례한다는 것이 중요하다(Time Advancement of Fluid Configuration를보라). 2 유체 문제에서의 난류와 임의의 스칼라양들은 이들이 두 유체경계면을 지나서 연속적으로 분포되어 있다고 가정하고 있기 때문에 이 범주에 속하지 않는다.

Time Advancement of Fluid Configuration시간 증가에따른 유체형상

유체분율F의 스칼라 이류는 식(10.362)에의해 지배된다. 식(10.362)이 한 계산셀 상에서 적분될 때 한 셀에서의F의 변화는 셀 면들을 통과하는F유속량의 합으로 치환된다. 그러나 유체경계의 뚜렷한 선명도를 보존하기 위해 특별한 주의를 기울여야 한다. 여기에서 사용된 방법은 공여-수용 유속량 형태를 사용한다. 이의 본질적 개념은 근사적 경계형상을 설정하고 유속량을 계산하는데 이 형상을 사용하기 위해 유속량 경계의 상류및 하류에서F의 정보를 사용한다는 것이다.

여기서 사용되고 있는 Volume-of-Fluid (VOF) 기법에서 사용되기 위해 개발된 기본 방식은 시간 단계δt기간 동안에 셀면을 통하여 x방향으로 유속되어질F의 양을 고려함으로써 이해될 수있다[HN81]. 단위 단면적당 이 면을 통과하는 체적 유속은L = uAxδt이며 여기서u는 면에 수직한 속도이고Ax는 유동에 열려있는 면적이다. u 의 부호는 공여 와 수용 셀을 결정한다, 즉 각기 체적을 잃거나 얻는 셀들 예를 들면, u가 양수이면 상류 또는 좌측의 셀은 공여셀이고 하류나 우측의 셀은 수용셀이다. 한 시간 단계에 셀면을 통과하는 F의 유속량은 δF와 막히지 않은 단면적(δy, δz)의 곱으로 정의되며, 여기서

   (10.363)

이며,  여기에서

  (10.364)

이다. 단일 첨자는 수용(A) 또는 공여(D) 셀을 표기한다. 이중 첨자AD는 밑에 설명되 듯이 유동방향에 상대적인 경계면의 향배에 따른 A또는D 중의 하나를 뜻한다. FDMFD의 최대값이며 공여셀 상류의 이웃 셀에서의F값이다.

간단히 식 (10.363)의MIN기능은 주어져야 될 것보다 더 많은 F 유속량이 발생하는 것을 방지하며MAX기능은 유속되어야 할 공간(1-F)의 양이 실제 이용가능 공간 양을 능가하면  추가F유속량CF을 고려한다. 밑의 일련의 그림은 식(10.363)의 도식 설명을 제공한다. 공여및 수용 셀들은 수직 셀 표면을 통과하는 유속량에 대해 그림 (a) 에서 정의된다. AD = D일때 유속량은 일반 공여셀 값,

  (10.365)

이며,   여기서 공여셀에서의F값은 유체에 노출된 가능한 유동면적 부분을 정의하는데 이용된다(하기 그림(b)참조). 수치 안정성으로 인해  |L|이V δx보다 작아야하므로 이 경우에 공여셀을 비우는 것이 불가능하다.

그림 10.14 F 이류에 사용된 자유표면 형태의 예제. 공여-수용 배열이 점선이 이류되는 전체 체적의 좌측경계를 가리키는 곳(a)에서 보여진다. (b-d)에서 보여지는 음영 지역은 실제로 이류된F의 양이다.

AD = A 일때 수용셀 내의F값은F가 이동하는 유동면적의 부분을 정의하는 데 이용된다. 위 그림(c)에서 공여 셀내의 모든F 유체는 점선과 유속 경계사이의 모든 것이 수용셀로 이동하므로 유속화된다. 이는 식(10.363)에서의MIN테스트의 예제이다. 위의 그림(d) 에서는FA|L|양 보다 더 많은 유체F가 유속되어야 하므로 이는MAX테스트의 예이다.

수용 또는 공여셀이 사용되는 지의 여부는 평균 표면의 향배에 달려있다. 표면이 자체에 수직한 방향으로 이류될 때 수용셀이 이용된다(셀표식 NF에 의해 정의되는 것같이); 그렇지 않으면 공여셀 값이 사용된다. 그러나 수용셀의F가 없거나 공여셀의 상류셀이 비어 있으면 수용셀F값은 표면의 방향에 상관없이 유속량을 결정하기 위해 사용된다. 이는 임의의 유체F가 하류의 빈 셀에 들어가기 전에 공여셀이 거의 가득차 있어야 한다는 것을 뜻한다.

표면의 방향성에 대해 테스트하는 이유는 수용 셀이 항상 유속량을 계산하는데 이용되면 표면파의 부정확한 고 경사도가 발생하기 때문이다. 예를 들면, 양의 x방향으로 움직이는 작은파를 가지는 수평표면을 고려해 보자.  F의 하류(수용) 값에 의거한 유속량은 결국 계단 형태의 불연속성을 가지는 고 경사도의 파를 이루게 된다. 사실상 수용셀 방식은F의 역 확산(즉, 음의 계수를 가지는 확산같은 이송)을 소개하므로 수치적으로 불안정하다. 그러나 불안정성은 유속량 정의에 이용되는MIN 과 MAX테스트때문에 무제한의 값으로 증가하지는 않는다. 반면에 표면이 자체에 수직하게 이류하면F의 단계함수 형태를 유지하기 위한 급격한 경사도는 정확히 원하는 바이다.

일단 유속량이 위의 방식에 의해 계산이 되면 공여셀에서 차감되고 수용셀에서 추가되는 F유체의 양을 얻기 위해 이는 유속 경계 면적에 의해 곱하여 진다. 이런 식으로F에 의해 정의된 유체의 양이 보존된다. 이류과정이 망 내의 모든 셀 경계에서 반복될 때 결과적F값은 식(10.19)을 만족시키는 시간 전진의 값에 상응하며 뚜렷한 경계면을 그럼에도 유지한다.

실제에서 유체운동이 표면의 반복된 뭉침과 파괴를 일으킬 만큼 과격하고 특히 이런 과정이 한 셀보다 작은 크기의 규모로 발생한다면 경계면의 선명도는시간이 지나면서 다소 악화된다. 특별한 경우에 대한 다중의 테스트가 경계면 추적능력을 더 향상시키기 위해 상기에서 기술된 기본F-이류의 루틴에 포함되어 있다. 또한 자유표면 문제에서 내부유체 지역에서의 부분적 공간(즉 1보다작은F값)을 폐쇄할 수 있는 방식이 추가되어 있다. 이런 지역들은 가끔 유체표면이 다른 표면과 또는 고상물체와 부딪힐 때 나타난다. 이 폐쇄방식은 또한 유체가 물체에 부딪힐 때 압력의 급격한 증감의 완화를 돕는다. 입력 표식(IFPK)은 이 선택을 활성화 또는 비활성화 시키는데 이용 되는데 이는 단지 공간의 압축성이 없는, RCSQV = 0.0, 자유표면 문제에서만 이용될 수 있다. 유체가 붕괴하는 극단적인 경우에 자유표면을 보존하기 위해 추가로 경계면 선명화가 더해질 필요가 있을 수 있다. 이 변수IFPK 가 보통 약간 더 뾰죽한 자유표면을 발생시키므로 일반적으로는 권장되지 않지만, 이 선택을 하는데 역시 이용된다.

Other Scalar Advection   다른 스칼라 이류

F이류를 위해 이용된 알고리즘은 압축 유동의ρ ρI로 확장된다. 경계면 위치, 그러므로 유속이 될 각 유체의 양은 위에 기술된 바와 같이 결정된다. 이 때 ρ (또는 ρI)의 유속량은 쉽게

  (10.366)

나 동등하게

   (10.367)

로 계산되며, 여기서Φ1 와 Φ2는 2 유체및η = δF/|L|와 연관되어 있다. 양쪽 모두의 유체에 분포되어있는2 유체 문제에서의 난류량에 대해, 유속량은 단순히Φ|L|와 같다. Second-Order Monotonicity Preserving Method에서 기술된 2차 단조 보존 이류방식은 또한 유체분율, F, 밀도, r, 및/또는 에너지, I ,의 이류에도 적용될 수 있다.

Determining Surface Normals and Cell Flags  표면 법선 및 셀 표식의 결정

자유표면 경계조건및F함수의 이류의 응용에서 표면에 근사적 법선방향을 지정하는 것이 필요하다. 한 표면 셀에서의 근사 법선은 표면의 내부로 수직한 방향에 가장 가까운 이웃 셀을 확인하는 정수값(NF배열에 지정되는)에 의해 기록된다. 이 방향은 빈 이웃셀로부터 멀어지는 방향을 가리켜야 한다. 표면셀이 하나 이상의 빈 이웃을 가지면 이 때 선택된 방향은 반대편 이웃 셀에서 가장 큰F값을 가지는 방향이다. 표면장력이 요구되면 더 정확한 표면 법선을 계산하는 것이 필요하다. 이를 위해 경계는 방향에 의존하는 단일 함수X(y,z), 또는Y (x,z), 또는 Z(x,y)로써 지역적으로 나타내질 수 있다. 예를들면, 근사 법선이z방향을 가리키면 그 때 경계는Z(x,y)에 의해 표시되고 다섯개 열(i, j), (i ± 1, j) 그리고 (i, j ± 1)에 대해 Z값을 계산한다. 이 양들로부터 δZ/δx δZ/δy및 표면장력을 위한 이의 2차미분에 대한 유한 차분근사가 계산될 수있다.

k의 단계에서 Zi,j 대해 사용된 근사는 세개의 셀 열 상에서의 합이며,

   (10.368)

여기서δz에의 한 ‘*’ 상첨자는 이 양이 표시된 k-단계에 대해δz이거나 0이라는 것을 뜻한다. 예를들면, δzk*−1에 대해 다섯 열 중에 어느 하나에서k-1 와 k단계 사이의 유동면적이 0일 경우에 0의 값 이 이용된다. 이 조치는 물체 및 배플 근처에서 합리적인 표면 구배와 곡률 값을 얻는데 필요하다.

Surface Location within a Cell   셀내 표면의 위치

일단 내부 법선방향이 결정되면 표면위치는 세 셀 열에서 적절한 높이까지 확장되는 이 좌표 방향에 수직한 평평한 표면에 의해 정의된다. 실제로는 미분계산에서 사용된 것으로부터 다소 수정된 열의 높이를 사용하는 것이 필요하다. 내부법선 방향의 이웃셀의 중심으로부터 측정되는 높이는 표면셀 유체에 이웃 셀 유체 높이의 반을 더한 것과 같다. 이 조치는 이웃셀이 거의 비어있더라도 표면이 이웃셀 중심보다 높게 위치하는 것을 확실하게 해준다.

Cleanup of Misty Fluid Regions   분무형 유체지역의 제거

어떤 응용에서는 해의 좋은 정확도를 얻기 위해 부드러운 자유표면 형태를 유지하는 것이 중요하다. 예를들면, 제대로 해결이되지 못하거나 들쑥날쑥한 자유표면은 표면장력의 평가시 수치적 ”잡음”을 일으킬 수 있다. 교대로 이는 더 들쑥날쑥한 자유표면 형태를 초래할 수 있다. 이런 변동은 궁극적으로 작은 양의 유체가 여러 개의 인접 셀들에 분포되는 “분무”형 유체지역의 형성으로 나타날 수있다. 이런 지역에서 자유표면이 만족스럽게 해결되지 않으므로 표면장력 해의 정확도는 급격히 악화된다.

“분무”지역이 수치해석의 질을 저하시킬 수있는 또 다른 예는 과도한 튀김과 자유표면의 부서짐이 있는 유동에서 발생한다. 이 모든 지역내의 작은 유체방울들의 움직임을 기술하려는 시도는 비효율적일 수있고 아마 불 필요할 것이다.

이런 상황을 방지하기 위해 인위적으로 “분무”형 유체지역(즉, 유체분포를 제거함으로써)을 정화하는 방편이FLOW-3D에서 주어진다. 이 수정은 변수FCLEAN에 의해 조절된다. 한 셀내의 유체는 그 안에서 그리고 그 모든 이웃셀에서의F값이FCLEAN밑으로 떨어지면 폐기된다. 일반적인FCLEAN값은 0과1사이이다. 현저한 경계를 갖는 2유체 문제에서 이과정은 해의 대칭성을 유지하기 위해 두 유체에 다 적용된다.

이 제거 알고리즘의 사용은 단지 유동이 주로 표면에서의 조건에 의해 조정되는 극심한 자유표면 변형이 있는 유동에 대해 권장된다. 일반적으로 이의 적용은 전체 유동체적에서 단지 작은 에러(1%보다 작은)를 유발하며 정확하고 효과적인 해를 준다.

Bookkeeping Adjustments  부기조절

상기방식으로 결정된 새 F값은 가끔 약간 0보다작거나 약간 1보다 큰 값을 갖는다. 그러므로 스칼라 이류계산을 마친 후에 0보다 작은F값을 0으로 그리고 1보다 큰F값을 1로 되돌리기위해 망을 통하여 한 경로를 다시 한다.

F 의 추가조절은 값이 0이나 1에 가까울 때 이루어진다. 한 셀이εF보다작거나1 − εF 여기서εF = 10−6 보다 큰 F값을 가지면 이 때F는각기 0이나 1로 재지정된다. F가 0으로 재지정될 때 모든 인접한 꽉 차있는 셀들은 표면 셀들이 된다.

Cleanup of Misty Fluid Regions에서 된 조절을 포함하여 전체유체체적을  허용된 범위로 유지하기 위해F 의 다양한 조절을 통해  제거및 추가된 전체 유체 체적은, 절대량일 뿐만아니라 모사시 영역을 통과한  전체 유체 1의 체적의 백분율인, cumulative volume error으로 기록된다. 이는 해 요약 파일의 긴 프린트에서 VCHGT로 뿐만 아니라flsgrf파일에서General History 데이터 카탈로그에 쓰여진다. 전체 유체 체적VL또한 인쇄된다. 일반적으로 축적 체적 에러는 전체 유체 1 체적의 1퍼센트보다 작아야한다.

전체 에러에 추가로 축적 체적 에러는 모든 셀에서 계산되고 공간양으로flsgrf파일에서 기록된다. 모사기간 동안에 이 체적은 조건이 허락할 때, 즉F값이 허용된 범위 안에 있고 셀의 상태(비어있던가, 표면이던가, 가득차 있던가)가 이 결과로 변하지 않으면 유체에 추가(이 체적이 음이면 제거)될 수있다. 이는 전반적 유체 체적 보존을 향상시키는 데 도움이 된다.

 

[FLOW-3D 이론] Numerical Approximations 수치근사 – 압력 솔루션 알고리즘

Numerical Approximations 수치근사

Pressure Solution Algorithm 압력 해 알고리즘

질량보존의 수치 처리는 압축성과 비압축성에서 상당히 다르다. 그러나 어느 경우든지 적합한 질량 방정식은 셀내의 압력을 결정하고 속도를 갱신하는 알고리즘에 이르게한다. 압축 유동이나 제한적 압축 유동(내재적 압력 해법을 사용해야 하는)에서 연속방정식(10.6) 또는 (10.8)은 셀에서의 압력과 속도의 타원조건으로 직접 해석될 수있다.  압축유동에서는 연속방정식(10.1)이 포물선 방정식으로 즉, 시간에 대해 전진하는 알고리즘에 의해 해석된다. 이 때 압력은 상태방정식 밀도가 갱신되는 셀 밀도와 같게함으로써 결정된다.  이 경우 시간단계 크기가 음파 전달에관한 안정성을 확실하게할 만큼 충분히 작다면 속도는 갱신될 필요가 없다. 압축유동에 대해 내재적 선택을 사용하면 더 큰 시간단계를 허용하나 충격파나 저밀도 파형에는 덜 정확할 수 있다.

때때로 인위적 제한 압축성을 유체에 추가하는 것이 해에 상당한 에러를 일으키지 않고 수렴을 증진시킬 수 있다. FLOW-3D에서 이는IMP = 2(디폴트는 1)로 지정함으로써 자동적으로 이루어진다. 상세한 모델 내용을 위해 http://users.flow3d.com/tech-notes/default.asp의 사용자 주소에서Flow Science Technical Note #55를 참조하라.

Incompressible SOR Method  비압축 SOR방식

식(10.315)으로부터 계산된 속도는 제한적 압축 연속 방정식(10.8)에 대한 다음의 이산화 근사식을 만족시켜야하며,

(10.352)

여기서XCi는 셀i의 중심의x위치이다. 원통좌표계에서는CYL = 1.0 및 Ri = XIM1/XCi 이며XIM1는 망내 마지막 실재 셀의 바깥 가장자리의 반경(x위치)이다. 데카르트 형상에서는 모든i 에 대해Ri = 1.0 이고CYL = 0.0이다. 항RSOR 은 셀내의 유체 체적 소스를 뜻한다. 압축성계수1/(ρc)2는 2유체 문제의 공식에 의해 계산된다.

 (10.353)

1유체 문제에서는 단지 두번째 항만 존재한다. 첨자l v는 각기 유체 1과 2를 뜻한다. 제한적 압축성은 두가지 목적, 물리적및 수치적, 으로 사용될 수 있다. 이 두 입력변수RCSQL = 1/(︀ρc2)︀ 와 RCSQV = 1/(︀ρc2)︀를 지정함으로써 유체의 제한적인 물리적 압축성을 모델할 수 있다. 또한 RCSQL 와 RCSQV를 적절히 지정함으로써 물체 경계상에서의 자유표면의 붕괴에 의해 종종 발생하는 수치적 압력 파동의 효과를 완화시킬 수 있다.

속도가 식 (10.352)을 만족시키기 위해 압력 그러므로 유체가 차지하고 있는 셀 내의 속도를 조절하는 것이 필요하다. 이는 둘 중 하나의 방식으로 행해진다. 가장 간단한 방법은 successive over-relaxation(연속가속완화) (SOR)반복 과정이다. 계산망을 망내의 첫번쩨 비 경계셀에서 시작하여 하나씩 쓸어나간다. 쓸림은 먼저i에 대해 시행되고 다음에j그리고 마지막으로k값에 대해서 되어진다. 계산은 단지 유체를 포함는 유체가 빈 이웃이 없는 셀들에 대해 실행된다.  셀(i, j, k)에서의 속도가 식(10.352)을 만족시키기 위해 필요한 압력변화는 아래와 같으며, 여기서 S는 식(10.352)의 좌측이다.

(10.354)

식(10.354)은 단순히S = 0를 이루기 위해 필요한p값을 생성하는 완화과정의Newton형태이다. 각 셀에서S를 평가하기 위해 사용된 속도 값은 반복과정에서 사용 가능한 가장 최신의 값이다.  식(10.354)으로부터의 결과를 이용하여 셀압력의 세 추정치는

(10.355)

이며, 셀의 면들에 위치한 속도들의 새 추정치는 다음 식과 같다.

      (10.355)

여기서 여기에 나타나는 속도들 또한 반복중에 가장 최신의 값들이다. 반복과정을 시작하기 위해 식(10.315)으로부터의 새 추정 속도들은 전 시간단계로부터 남아있는 압력값과 함께 이용된다. 물론 면적이0인 곳에서의 속도는 이 단계에 수정되지 않는다.

자유표면을 가지는 셀들에서, 즉 유체가있지만 하나 또는 더 많은 빈 이웃 셀들을 가지는 셀에서는 다른 과정이 이용된다. 이러한 셀들에서 요구되는 경계조건은 압력이 표면에서 지정된 값, 즉 표면에서의 ps 이다. 표면압력은 이웃한 void영역 압력, PR, 과 표면장력 압력, PS,의 합과 같도록 지정되며

(10.357)

여기서n은 인접한 빈 공간의 색인이다. PS의 평가는Surface Tension with Wall Adhesion에서 기술된다. 표면압력, ps 은 셀내의 정압분포를 가정하여 표면셀의 중심에서의 압력pi,j,k 으로 외삽하여 해석에 적용된다.정압변화는NF에의해 정의된 바와같은 표면에 수직한 방향에서의 순수가속도에 의존한다. 이 표면 셀 압력은 압력 반복동안에 변하지않으며 고정 경계값으로 취급된다. 이런 방법으로 셀내의 자유 표면의 실제 위치가 확실하게 고려된다.

고압, 단열 기포가 존재할 때[Hir92], [BC94] 수치 불안정성이 공간지역 압력의 외재적 근사로 인해 발생할 수 있다고 알려져 있다(단열 기포모델은Variable Pressure (Adiabatic) Void Region에서 기술되어 있다). 이러한 어려움을 없애기 위해FLOW-3D에 내재적 기포모델이 추가되어 있다. 이의 목적은 기포압력 변화가 한 사이클의 마지막에서 계산되도록 기대하고 이를 통상적인 압력-속도 반복과정에 포함하는 것이다. 내재적 기포 모델은 기포의 “강성도”가 너무 크지 않다면 잘 작동한다. 다른 말로, 강성 기포는 기대되고 그리고 실제압력 변화가 한시간 단계내에서 너무 다른 기포를 뜻한다. 이런 강성 기포들이 발생하면 해석은 매 사이클 마다 큰 압력변화와속도와 다른 양 들에서 이에 상응하는 커다란 변동을 가질 수 있다(상세 내용을 위해 Ref. [Hir92] 참조하라).

2 유체문제에서 모든 셀들은 유체로 가득 차 있다고 간주된다; 즉, 압력과 속도는 반복하는 동안에 모든 유체 셀 내에서 조절된다. 표면장력이 작용되면 압력PSi,j, 은 두 유체중에 하나에만 작용해서 표면장력으로 의한 경계면을 통과하는 압력에서의 불연속성이 유지된다.

완전한 반복은 식(10.354), (10.355) 및 (10.356)에따라 모든 유체가 가득찬 셀내의 압력및 속도를 조절하는 것으로 이루어져 있다. 반복시의 수렴은 모든 셀들이 어떤 일정 작은 수인EPSI ·VFi,j,k 보다 작은S값을 가질 때 이루어진다.

EPSI의 값은 자동적으로FLOW-3D에 의해 각 시간 사이클에서 시간 단계 크기의 함수로 계산된다. 이 알고리즘은 입력변수EPSADJ의 값이 양수이면 원용된다. 선택적으로 한 EPSI의 상수가 한 계산 과정에서 사용 가능하다.

어떤 경우에는 반복의 수렴이 식(10.355)에서의δp를 완화인자OMEGA로 곱함으로써 가속화될 수 있다.

OMEGA 는 1.7 또는 1.8이 최적값이다. 어떤경우에도 는 2.0이 넘어서는 안되는데 이는 이럴 경우 불안정한 반복이 발생하기 때문이다. 압축성유동에서OMEGA는 1.0으로 지정된다. 또한 시간 단계 크기가 상당히 대류 안정성 한계보다 작을 때 비압축성 유동에 대해OMEGA는 1.0으로 사용하는것이 권장된다. 이는 해에서의 잠재적 압력 잡음을 감소시킬 수 있다.

Incompressible Line Implicit SADI Method 비압축성 선 내재적 SADI 방식

앞에 언급된 압력을 계산하기 위한SOR반복법은 간단하고 많은 문제에서 잘 작동한다. 그러나SOR방식의 수렴이 상당히 느려지는 경우들이 있다. 예를들면, 한방향으로의 셀 크기가 다른 방향으로 보다 훨씬 큰 망은SOR압력 완화가 횡방향의 작은 셀크기에 의해 제한되므로 큰 셀방향으로 느린 수렴성을 보여줄 것이다.

이런 형태의 더딘 수렴에 대한 보완은 더 작은 셀크기의 방향에서 더 내재적인 해석 방식을 사용하는 것이다. 이런 목적으로 수정된Alternating-Direction-Implicit (SADI)방식이 개발되었다. SADI는 망 셀의 한 i, j, 또는 k열을 따라 압력에 대한 표준 3중 대각해법에 근거한다.   이 해법은 주기적 경계를 포함하는 모든 경계조건에 적용 가능하다.

주기적 경계가 원통좌표계에서 방위각의 방향에 사용될 때IADIY=1로 지정함으로써-단지 이 방향으로만- ADI압력 해법을 사용하는 것이 권장된다. 이는 가끔 발생할 수 있는 압력과 속도해에서의 수치 잡음을 제거하는데 도움이 될 것이다. 다른 방법으로는OMEGA= 1.0을 지정하여 상향 완화를 잠금으로써 잡음을 감소시키는데 도움이될 수 있다.

SADI방식이 z방향으로사용되면 반복은 모든 i j색인을 거치고 각(i, j)쌍에 대해k-색인 방향에서 압력에 대해 내재적으로 해석하는 것으로 이루어진다. 이웃 열들에서 필요한 압력 값들은 표준ADI에서는 반드시 항상 그렇게 실행되지는 않지만 최신의 반복값을 항상 취하는데 이는 수렴을 향상시킨다.

SADI방식은 방향의 어떤 조합으로도 사용될 수있다: 어느 하나 또는 둘 또는 셋 모두. 이는 더 비용이 드는 내재적 소해가 단지 전체 반복의 수렴을 향상시키는데 필요한 방향에서만으로 제한될 수 있다는 것을 뜻한다.

단지 한 또는 두 방향으로만 내재적으로 처리될 때 열간의 상향 완화는SOR상향 완화에서 사용되는 변수OMEGA에 의해 조절된다. 일반적으로SADI는 이 변수에 그렇게 민감하지 않으며 디폴트 값OMEGA=1.7은 보통 만족스럽다.  SADI 가 세 방향 모두에 사용될 때 이 경우 상향 완화가 최대값에서 고정되므로 변수OMEGA 는 영향을 미치지 않는다.

Compressible Solution Method 압축성 해 방식

압축성 유동에서 셀 압력은 연속방정식 밀도를 상태방정식으로부터 결정되는 밀도와 동일시함으로써 결정된다. 이 방정식에는 SOR 와 SADI 방식 둘 다 이용 가능하다. 두 알고리즘에서 식(10.354)에서 사용되는 반복 함수는 아래와 같이 정의된다.

  (10.358)

인자∇ · (uAΦ)는 시간n+1에서의 속도에 의거하여 셀의 압축및 팽창을 조절하며 여기서

 (10.359)

C1은 유체 1에서의 음속이다. 외재적 해석 방식에서 외재적 모멘텀 방정식으로부터의 속도는 이며, 이는 가 셀압력에 대한 반복에 독립적이라는 것에 주목한다. 상태 방정식밀도는

 (10.360)

로 정의된다.

상태방정식 밀도는 압축이나 팽창에 의한 시간 정도 n에서 n+1사이에 발생하는 에너지변화에 대해 조절되지 않는 것에 주목한다. 반복이 매시간 단계에서 되풀이 되므로 에너지를 갱신하지 않는데서 비롯되는 오차가 기껏해야 한시간 단계 쳐지고 시간단계 안정성을 이루는 목적을 위해서도 의미가 없다. 식(10.354)에서 사용된 양δS/δp은 각 사이클에서 한 번씩 평가되고 저장되는 항DSDPU = δDTi,j,k/δpi,j,k을 필요로 한다.

SADI 해법은 비압축성에서와 같이 진행한다. 반복함수는 아직 3중 대각 시스템(횡방향에서의p대한 최신의 반복값을 유지하며)으로 처리될 수 있다. 모든 셀에 대해 어느 경우에도 수렵에 도달한다.

  (10.361)

SOR해석 알고리즘에서 상태방정식에서의 밀도에 대한 압력의 의존도가 비선형(코드의 사용자 수정에 의해 가능)일 경우에 유용한 선택이 주어진다. 다음셀로 전진하기 전에 식(10.361)의 조건을 만족하기 위해 셀내의 압력을 변화시키도록내부 반복이 실행된다. 사용자는 내부 반복의 최대수(IITMX)와  완화인자 OMEGA를 지정해도 된다. IITMX > 1이면OMEGA = 1.0이 권장된다.

지정된 속도 및 지정된 압력경계 조건은 계산영역 내에 이용 가능한 반복법중에 어느 방법으로도 계산하기 어려운 균일한 압력변화를 형성할 수가 있다.  이런 경우에 전반적으로 균일한 압력조절을 주기 위해 추가 알고리즘이 압축 해 과정에 주어진다. 이 선택은 변수IPUN를 1로 지정함으로써 활성화된다. If IPUN = 0이면 균일 압력변화는 평가되지 않는다. 이 균일 압력조절은 셀들에 대해 Si,j,k와 셀 체적의 곱을 합하고 이 결과를δS/δp와 셀 체적의 곱의 셀들에 대한 합으로 나눔으로써 계산된다.

GMRES Pressure-Velocity Solvers GMRES 압력-속도 해법

새로운 압력-속도해법이FLOW-3D [AMS90], [BBC+94] [Saa96] 에서 실행되고 있다. GMRES는  일반화된 최소 잔류 방법을 뜻한다. GMRES솔버에 추가하여 또 새로운 선택적 알고리즘- 일반화된 짝 구배(GCG)알고리즘-이 새 GMRES솔버에서 점성항을 위해 실행되고 있다. 이 새 솔버는 많은 범주의 문제들에 대해 아주 정확하고 효과적이다. 좋은 수렴성, 대칭 및 속도성을 갖는다; 그러나SOR 이나 SADI방법보다 더 많은 메모리를 사용한다. GMRES 솔버는  과소 또는 상향 완화를 사용하지 않는다.

사용된 방법들에 대한 상세한 내용은http://users.flow3d.com/tech-notes/default.asp에 있는 사용자 주소상의 Flow  Science TN68에서 찾을 수 있다.

[FLOW-3D 이론] Auxiliary Model/Fan and Impeller Model 팬과 임펠러모델

팬과 임펠러모델

FLOW-3D 에서 정의된 팬과 임펠러 모델은 날개의 회전율이 유체가 정상상태에 이를 때까지 많은 회전이 필요할 때 사용될 수있다.

이 모델은 회전과 축속도 성분을 유도한다. 팬이나 임펠러는 구역을 정의하나 실제 물체의 막힘효과가 없는 “phantom” 물체의 형태로 정의된다. 일반적으로 이런 물체들은 회전 날개에 의해 휩쓸어지는 외경 R , 내경 r, 두께 L 인 직원통으로 가정된다.

형상 이외에 팬이나 임펠러의 성능을 결정하는 나머지 변수들은 회전율 Sd, 날개가 얼마나 효율적으로 유체에 운동을 가하는 지를 결정하는 조절계수 Ad, 그리고 유도된 축방향 유동량을 조절하는 계수 Bd 들이다.

팬이나 임펠러의 성능은 상세한 날개의 크기와 형태 그리고 날개의 수에 의존하므로 경험식으로부터 Ad Bd 값을 결정하는 것이 최선이다. 이 장치들의 제조자들은 가끔 이 값들을 장치 통과시 압력 저하대 이를 통과하는 평균 유량의 그림인 소위 “성능곡선”으로 특성화한다.

Typical performance curve (solid line) and |f3d| approximation (dashed)그림 10.4 전형적 성능곡선(실선)과 FLOW-3D 근사치(점선)

FLOW-3D 에서 사용된 모델을위한 성능곡선은 회전 모멘텀소스를 장치 두께를 통과시의 등가 압력저하와 전체 단면을 통과하는 평균유량을 연관시켜 유도될 수 있다. 이 결과는:

(135)\Delta p = \rho L{A_d}\left( {\frac{2}{3}{S_d}{B_d}R\left( {1 - \frac{{{r^3}}}{{{R^3}}}} \right) - \frac{Q}{{\pi {R^2}}}} \right)

이 식에서 ρ 는 유체밀도이며 Q 는 순수 유량이다. 이 관계식은 다음에 의해 주어지는 위 그림의 y-절편 ∆ρ0 와 x-절편 Q0 를 갖는 선형 성능 곡선을 준다:

(136)\Delta {p_0} = \rho L\left( {\frac{{{Q_0}}}{{\pi {R^2}}}} \right){A_d}, \quad {Q_0} = \frac{2}{3}\pi \left( {{R^3} - {r^3}} \right){S_d}{B_d}

이 관계 및 주어진 회전율 OSPIN = Sd를 이용하여 OADRG = Ad 와 OBDRG = Bd 변수들이 원하는 성능 곡선에의 선형근사를 주도록 계산 될 수 있다.

[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.