어떤 물질은 전체응력 상태에서 압력과 점도외에 추가의 응력이 존재한다: 탄성응력. 점성응력은 물질의 변형율의 함수인 반면에 탄성응력은 물질의 지역 전체 변형의 함수이다. 이런 물질은 탄성응력이 항복 한계까지 증가하며 그 때에 점성유체 같이 거동하는 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 내의 다른 온도 의존 양에서와 마찬가지로 구간 선형 보간이 두 데이터 점 사이에서 가정된다.