Adjoint Optimization을 이용한 차량용 후처리시스템의 입구형상 최적화

1. 서 론

날로 강화되는 배기규제에 대응하기 위하여 배출되는 유해가스를 줄이기 위해 촉매 변환기내의 배기 유동분포개선 및 활성화 시간 단축을 위한 많은 노력들을 기울이고 있다.¹⁾

특히 촉매변환기에 입구부 형상을 변경시켜서 촉매담체로 유입되는 유량을 균일하게 하는 형상최적화는 여러가지 방법으로 진행되고 있다.

하지만 그 방법들은 인자수가 많아지면 조합수가 기하급수적으로 증가하여, 얼마나 설계인자를 적게하여 조합수를 최소화시킬 수 있느냐가 관건이다.

그래서 전역적인 방법론은 설계인자가 적은 경우에는 유용한 방법론이나. 설계인자수가 많아지는 경우는 실제 제품설계 과정에서 적용하기는 어려움이 많다.

하지만, 지배방정식의 범함수를 이루는 Adjoint 변수에 대한 Duality 문제로 변환하는 기울기 기반의 Adjoint 최적화는 원방정식 1회와 Adjoint방정식 1회를 구하면 국부적인 부분에서의 설계민감도도 손쉽게 구할 수 있으므로, 최적화수행회수를 크게 줄일 수 있다.

그래서, 본 논문에서는 전체적 형상에 대해 영향을 미치는 설계인자 3개를 전역적인 방법으로 이용하여 최적화를 수행한 후에, 많은 설계인자를 가지는 국부적인 형상변경에 대한 최적화는 Adjoint 방법을 이용하여 추가 수행 방식을 이용하였다.

그 결과, 기존의 방식에서는 파악하지 못한 국부적이 형상개선안을 파악할 수 있었고, 전역적인 최적화이후 추가적인 배기계 유속균일도의 성능향상을 얻을 수 있었다.

2. 최적화 방법

이 장에서는 본 논문에서 사용하는 두가지 최적화 기법으로 실험계획방법과 Adjont방법에 대한 내용을 소개한다.

전체적인 최적화 프로세스는

• 1) 1단계로 실험계획방법을 이용하여 컨셉설계상의 모델 중 전역적으로 대략적인 적접 조합을 선택후,
• 2) 2단계로 그 설계인자 조합으로 이루어진 형상에서의 Adjoint 최적화를 수행한다.

실험계획법(DOE)은 직교하는 설계변수의 범위를 지정하여 그 범위내의 평가결과의 최대 혹은 최소값을 산출하여 인자별 최적값을 산출하는 방법으로 전역적 최적화 방식중 가장 간단한 방식이다.

물론, 2차식의 근사모델을 이용한 RSM(반응표면법)이나 주변 관측값의 통계학적 상관강도를 이용한 근사식의 선형결합으로 계산하는 Kriging이 있으나, 좀 더 정확한 만큼 더 많은 데이터와 계산과정이 필요하게 된다.

2단계로 이루어지는 Adjoint 최적화를 통한 국부적인 최적화 과정이 있기 때문에 1단계의 전역적 탐색을 최대한 간단한 방식으로 채택하였다.

특히, 노이즈 인자를 가정하여 노이즈 인자에 Robust한 최적화가 가능하다는 장점이 있어서 1단계에서 적용하게 되었다.

이 최적화의 목적함수는 1/2단계 모두 촉매담체 입구에서의 속도균일도로 선정하였다.

이 논문에서는 1단계인 실험계획법에서는 간단하게 해석모델, 인자수준, 1차적으로 진행시킨 최적화 평가결과만을 간단히 소개하고, 다음 단계인 Adjoint 최적화 단계를 주로 다루려고 한다.

1단계인 실험계획법의 상세한 내용은 “다꾸치 방법을 이용한 디젤엔진용 후처리시스템의 입구부 형상 최적화”¹⁾에 제시되 있다.

2.1 실험계획법(다꾸치) 방법

설계 변수와 최종 제품의 유속균일도간의 관계 규명 및 최적화 작업에 사용될 실험 방법은 우선 실험횟수를 최소화하여 시간과 경제적 부담을 줄여야 하며, 실험에서 결정되는 설계변수가 잡음에 대하여 강건해야 하는 요구사항을 만족시켜야 한다.

이러한 실험 요구 조건을 만족시키기 위한 대표적인 실험 계획법인 다꾸찌 기법을 사용하였다. 다꾸찌 기법에서는 교란인자의 영향으로부터 품질을 강건하게 유지하기 위해서 SN비(Signal to Noise ratio)를 사용한다. SN비는 신호 대 잡음의 비율을 의미하며, 망대, 망소, 망목 특성으로 최적화가 가능하다.

2.1.1 컨셉모델

개념설계모델은 총 3가지 설계인자로 되어있다.

• (1) Passline Radius : Turbo-Charger 출구부에서 직선으로 돌출되어 “ㄱ” 형태로 꺾이는 배관형태의 곡률반경
• (2) Sphere Diameter : 유동이 분리되는 형태를 가지는 Sphere 지름
• (3) Forming Angle : 유동이 분리되는 지점의 Forming형상의 각도

Fig. 1

The concept model of exhaust downpipe

2.1.2 설계인자수준

컨셉모델의 설계인자별 수준은 외부 타부품과의 간섭을 피할 수 있는 공간적 제한사항을 고려하여야 하며, 그 한계치를 고려한 수준은 다음과 같다.

Table 1

The level of design parameter

2.1.3 노이즈 인자수준

제어하지 못하는 인자를 노이즈 인자를 설정하지만, 실험계획법에서는 조건변동에 Robust(강건한) 설계를 하기위해서, 해당 조건을 노이즈 인자로 선정하기도 한다.

이는 기울기기반의 최적화 기법과 차별되는 실험계획법의 적용 특징이기도 하다. 여기서는 엔진운행조건에 강건한 형상최적화를 목표로 하기 때문에, 엔진운행조건에 따라 달라지는 질량유입량, 배기가스온도, TC 회전수와 그에 따라 달라지는 출구쪽 압력, 촉매측 유동저항계수를 노이즈 인자로 선정하고 그 수준은 다음과 같다.

Table 2

The level of noise parameter

2.1.4 실험계획 탐색결과 및 S/N비 주효과도

인자별 주효과는 아래와 같이 관측치의 인자집단간 편차를 나타낸 것이다.

실험계획을 산정후, 다꾸치 방법을 통해 산출한 실험계획법에 따라서 관측을 한 다음, 주효과도를 분석하면, 최대값을 가지는 인자수준을 알 수 있다.

그 결과로 설계인자 1, 2, 3순으로 3수준, 3수준, 2수준으로 산출되었다.

Fig. 2

Main effect plot for SN ratio

2.2 Adjoint 방법

형상최적화를 반복적인 방법으로 접근하기 위해서는 설계인자에 대한 유동장의 민감도를 계산하는 것이 필요하다. 제어변수가 3~4개 이하인 일반적인 경우는 유한차분법을 이용하여 쉽게 구할 수가 있으나, 설계인자수가 증가할수록 지수에 비례적으로 해석 수행수가 증가하는 단점이 있다. Adjoint방법은 이 단점을 보안하는 대체방법으로 해석을 하고자 하는 식과 거기에 상응하는 Adjoint 변수로 구성된 범함수의 설계인자에 대한 도함수를 구하여 각각의 설계인자에 해당하는 민감도를 구할 수 있다.

이 Adjoint 방식의 목적함수 J(U, α)는 유동변수 U와 설계변수 α의 함수형태이다. 이 설계변수에 대한 목적함수의 설계민감도는 다음과 같다.

$\[ \frac{dJ}{d\alpha }=\frac{\partial J}{\partial U}\frac{dU}{d\alpha }+\frac{\partial J}{\partial \alpha } \]$

(1)

이 유동변수가 Navier-stoke 방정식 N(U, α)를 만족시키는 제한조건을 Adjoint 변수 λ를 도입하면, 여기서,

$\[ \frac{dJ}{d\alpha }=\frac{\partial J}{\partial U}\frac{dU}{d\alpha }+\frac{\partial J}{\partial \alpha }+{\lambda }^T\frac{dN}{d\alpha } \]$

(2)

$\[ \frac{dN}{d\alpha }=\frac{\partial J}{\partial U}\frac{dU}{d\alpha }+\frac{\partial N}{\partial \alpha }=0 \]$

(3)

이다.

따라서, 이 식을 재구성하면,

$\[ \frac{dJ}{d\alpha }=\left(\frac{\partial J}{\partial U}-{\lambda }^T\frac{\partial N}{\partial U}\right)\frac{dU}{d\alpha }+\left(\frac{\partial J}{\partial \alpha }-{\lambda }^T\frac{dN}{\partial \alpha }\right) \]$

(4)

만약, $$ \frac{\partial J}{\partial U}-{\lambda }^T\frac{\partial N}{\partial U}=0$$을 만족하는 λ를 구한다면, 다시 말해

$\[ \frac{\partial N}{\partial U}{\lambda }^T=\frac{\partial J}{\partial U} \]$

(5)

위의 식을 Adjoint 방정식으로 정의하고 이를 만족하는 Adjoint 변수 λ를 구한다면, 목적함수에 대한 설계민감도는 다음과 같다.

$\[ \frac{dJ}{d\alpha }=\frac{\partial J}{\partial \alpha }-{\lambda }^T\frac{\partial N}{\partial \alpha } \]$

(6)

여기서, 일반적 유동문제인 경우, $$ \frac{\partial J}{\partial \alpha }=0$$이고, 표면에서의 $$ \frac{\partial N}{\partial \alpha }$$는 축약형태로 산출할 수 있고, λ는 Adjoint 방정식을 통해서 구해지기 때문에 목적함수에 대한 설계민감도를 쉽게 구할 수 있다.

이 과정에 가장 시간이 많이 소요되는 과정은 Adjoint 방정식을 풀어서 λ를 구하는 과정이나, 이는 원식을 한번 푸는 것과 거의 동일한 시간만을 소요하기 때문에 설계인자가 증가하더라도 소요되는 시간의 증가는 상대적으로 거의 없다.³⁾

2.2.1 Steady 조건하의 Navier-Stork 방정식과 그에 대응하는 Adjoint 방정식

1) 원 NS 방정식

유동관련 최적화 문제에서는 유동의 흐름은 Navier-Stokes(이하 NS)방정식에 제한된다. 여기서 Steady상태의 NS 방정식은 다음과 같다.

$\[ \rho u_j{\partial }_j{u}_i+{\partial }_i\left(p+\frac{2}{3}k\right)-{\partial }_j\left(\left(\mathrm{\mu }+{\mathrm{\mu }}_t\right)\left({\partial }_j{u}_i+{\partial }_i{u}_j\right)\right)=0 \]$

(7)

$\[ {\partial }_j\left(\text{ρ}{u}_j\right)=0 \]$

(8)

여기서, ρ, u, p는 각각, 유체의 밀도, 속도, 압력을 의미하고, 표기법 ∂_ju_i는 $$ \frac{\partial u_i}{\partial x_j}$$를 나타낸다.

본 논문에서 사용하는 난류모델 k-ϵ model 방정식은 다음과 같다.

$\[ \text{ρ}{u}_i{\partial }_{i}k-{\partial }_i\left(D_k{\partial }_{i}k\right)-\text{ρ}G\text{+ρϵ}=0 \]$

(9)

$\[ \text{ρ}{u}_i{\partial }_i\text{ϵ}-{\partial }_i\left(D_{\text{ϵ}}{\partial }_i\text{ϵ}\right)-\text{ρ}{c}_{1}G\frac{\text{ϵ}}{k}\text{+ρ}{c}_2\frac{{\text{ϵ}}^2}{k}=0 \]$

(10)

여기서 D_k , D_ϵ은 특정 상수이며, G 는 난류에너지생성항이다.

특히, 원방정식에서 난류모델을 사용하는 경우에 Adjoint 방정식에도 포함시켜야 High Reynold 조건에서 정확성이 높다고 알려져 있다.

2) 연속 Adjoint NS 방정식

여기서, u*, p*, k*, ϵ*는 각각, 유체의 속도, 압력, 난류에너지, 난류소산율의 Adjoint 변수를 의미한다. 연속 Adjoint NS 방정식에 대한 상세한 유도과정은 관련 문헌을 참조하였다.^4,5)

$\[ \begin{array}{c}-\mathrm{\rho }{u}_j{\partial }_j{u}_i^{*}+u_i^{*}{\partial }_j\mathrm{\rho }{u}_i-{\partial }_j\left(\mathrm{\rho }\left(v+v_t\right)\left({\partial }_j{u}_i^{*}+{\partial }_i{u}_j^{*}\right)\right)\hfill \\ +{\partial }_j{p}^{*}+\mathrm{\rho }{k}^{*}{\partial }_{i}k+\mathrm{\rho }{\text{ϵ}}^{*}{\partial }_i\text{ϵ}+{2\partial }_j\left(\left(k^{*}+{\text{ϵ}}^{*}{c}_1\frac{\text{ϵ}}{k}\right)\mathrm{\rho }{v}_t\right.\left({\partial }_j{u}_i\right.+\hfill \end{array} \]$

(11)

$\[ {\partial }_i\text{ρ}{u}_i^{*}=\frac{\partial J_{\mathrm{\Omega }}}{\partial p} \]$

(12)

난류모델에 대한 Adjoint NS 방정식은 다음과 같다.³⁾

$\[ \begin{array}{c}-\rho u_j{\partial }_j{k}^{*}-{\partial }_j\left(D_k{\partial }_j{k}^{*}\right)+\rho {\text{ϵ}}^{*}{c}_{1}G\frac{\text{ϵ}}{k^2}\hfill \\ -\text{ρ}{\text{ϵ}}^{*}{c}_2\frac{{\text{ϵ}}^2}{k^2}-\rho v^{*}\left(c_{\mu k}^{*}\right)+\frac{2}{3}{\partial }_i\rho u_i^{*}=\frac{\partial J_{\mathrm{\Omega }}}{\partial k}\hfill \end{array} \]$

(13)

$\[ \begin{array}{c}-\text{ρ}{u}_j{\partial }_j{\text{ϵ}}^{*}-{\partial }_j\left(D_{\text{ϵ}}{\partial }_j{\text{ϵ}}^{*}\right)-\text{ρ}{c}_{1}G\frac{{\text{ϵ}}^{*}}{k}\hfill \\ +2\text{ρ}{c}_2\frac{\text{ϵ}}{k}{\text{ϵ}}^{*}+\text{ρ}{k}^{*}+\text{ρ}{\text{v}}^{*}\left(c_{\mu \text{ϵ}}^{*}\right)=\frac{\partial J_{\mathrm{\Omega }}}{\partial \text{ϵ}}\hfill \end{array} \]$

(14)

$\[ \begin{array}{cc}{v}^{*}=\hfill & -\rho \left({\partial }_j{u}_i^{*}\left({\partial }_j{u}_i^{*}+{\partial }_i{u}_j^{*}\right)\right)+\rho k^{*}{G}_k\hfill \\ \hfill & +\rho {\text{ϵ}}^{*}{c}_1{G}_k\frac{\text{ϵ}}{k}-\rho \frac{1}{{\sigma }_k} {\partial }_{i}k{\partial }_i{k}^{*}-\rho \frac{1}{{\sigma }_{\text{ϵ}}}{\partial }_i\text{ϵ}{\partial }_i{\text{ϵ}}^{*}\hfill \end{array} \]$

(15)

3) Adjoint 경계조건

연속 Adjoint 방정식은 다음과 같은 경계에서의 Adjoint 방정식을 포함하게 된다.⁴⁾

$\[ \delta u\left(\rho u_j{u}_i^{*}{n}_j-p^{*}{n}_i+\rho \left(v+v_t\right)\left({\partial }_j{u}_i^{*}\right)n_j-2\left(k^{*}+{\text{ϵ}}^{*}{c}_1\frac{ϵ}{k}\right)\text{ρ}{v}_t\left({\partial }_j{u}_i+{\partial }_i{u}_j\right)n_j-\frac{\partial J_{\text{T}}}{\partial u}\right)=0 \]$

(16)

$\[ \delta p\left(-n_i\text{ρ}{u}_i^{*}-\frac{\partial J_{\text{T}}}{\partial p}\right)=0 \]$

(17)

$\[ \delta k\left(\text{ρ}{u}_j{k}^{*}{n}_j+\left(D_k{n}_j{\partial }_j{k}^{*}\right)-\frac{\partial J_{\text{T}}}{\partial k}\right)=0 \]$

(18)

$\[ \delta \text{ϵ}\left(\text{ρ}{u}_j{\text{ϵ}}^{*}{n}_j+\left(D_{ϵ}{n}_j{\partial }_j{\text{ϵ}}^{*}\right)-\frac{\partial J_{\text{T}}}{\partial \text{ϵ}}\right)=0 \]$

(19)

$\[ \delta {\text{v}}_t\left(\text{-ρ}\left(n_{i}k{\partial }_i{k}^{*}\right)-\text{ρ}\left(n_i\text{ϵ}{\partial }_i{\text{ϵ}}^{*}\right)\right)=0 \]$

(20)

여기서 δ, n_j는 변분과 표면수직벡터를 나타낸다.

위 경계조건하의 Adjoint 방정식을 만족시키기 위해서 다음과 같은 Adjoint 경계조건을 채택하였다. 특히 원 NS방정식의 경계조건에 따라서 각 상태량의 변분값이 0인 경우와 아닌 경우로 나누어져서 경계에서의 Adjoint 변수의 값 설정이 달라질 수 있으므로, 원방정식의 경계조건과 그에 따른 Adjoint 방정식의 경계조건을 동시에 표시하였다.

또한, 유속균일도의 목적함수에는 압력, 난류에너지와 난류소산율에 대한 항이 없는 경우로 가정하여, $$ \frac{\partial J_T}{\partial p},\frac{\partial J_T}{\partial \text{ϵ}},\frac{\partial J_T}{\partial k}$$을 0으로 가정하여 경계조건을 최대한 단순화 시켰다.

Table 3

The boundary condition of adjoint velocity

Table 4

The boundary condition of adjoint pressure

Table 5

The boundary condition of adjoint turbulent kinetic energy

Table 6

The boundary condition of adjoint dissipation ratio

3) 목적함수

Adjoint 최적화에 사용되는 목적함수는 촉매담체입구 유동의 균일 정도를 파악하기 위한 유속균일도 지수로 유사한 형태로 다음과 같은 식을 사용한다.⁶⁾

여기서, U는 경계면 Γ의 평균유속을 의미하며, 경계면 Γ는 배기계 촉매전면이자 Down pipe의 출구면을 의미한다.

$\[ J=\frac{c}{2}\int {\left(u-\overline{U}\right)}^{2}d\mathrm{\Gamma }=\frac{1}{2}\int \frac{{\left(u-\overline{U}\right)}^2}{\overline{U}} d\mathrm{\Gamma } \]$

(21)

3. 배기계 Downpipe 최적화 결과

1단계의 최적화방법은 기존의 상용 유동해석코드의 해석 평가결과를 이용하여 최적화를 수행할 수 있는 반면에, Adjoint 최적화는 Adjoint 방정식을 풀기 위한 전용 Solver가 필요하다.

물론 상용프로그램도 Adjoint 전용 프로그램이 나와있긴 하지만, 아직 수렴성문제 및 다양한 문제적용 한계가 있는 부분이 있어서, 여러가지 수치적 기법을 수정할 수 없는 한계가 있다.

Adjoint solver는 공개된 유동해석 Library인 Open FOAM를 사용하여 작성 되었다. 이 라이블러리는 편미분형태의 지배방정식의 형태를 쉽게 Code화시키는 것이 용이하기도 하고, 많은 Adjoint관련 공개코드가 많아서 Code 구현이 용이하여 선택하였다.⁷⁾

본 논문에서는 원 NS 방정식은 기존의 상용코드로 풀고 수렴된 속도, 압력, 난류에너지, 난류소산율, 밀도등의 정보를 바탕으로 연속 Adjoint 방정식을 푸는 Adjoint solver를 직접 작성하였다.

물론, 원 NS 방정식도 같은 Code를 사용할 수도 있지만 제품개발시에는 특정회사의 인증된 상용 Code의 해석 평가결과만을 인정하는 경우가 많아서, 상용 Code를 이용하여 원 NS 방정식의 해를 구하고, Adjoint 방정식을 이용하여 그에 대응하는 Adjoint 변수를 산출하여 민감도를 산출하였다.

3.1 1단계 실험계획법 결과

결과분석에서 나온 Pipe 곡률 3수준, Sphere지름 3수준, 턱의 각도 2수준의 조합은 이미 실험 직교표에 L9으로 포함되어 있어서 따로, 최적 조합으로 해석을 재검증할 필요는 없다.

따라서 L9에 해당하는 모델과 기존의 양산모델과 그 결과를 Fig. 3과 Fig. 4에서 비교하였다.

이전 모델의 경우, 배기관의 형상에 의해 배기가스가 촉매담체 후방부 외곽쪽에 집중되는 것에 반해, 최적 조합에서는 생성된 형상에 의하여 최외곽으로 배기가스가 집중되는 것이 완화가 되는 것을 알 수 있었다.

Fig. 3

The shape and uniformity of 1st optimization

Fig. 4

The shape and uniformity of 1st optimization

3.2 2단계 Adjoint 최적화 결과

2단계 Adjoint 최적화를 이용한 형상변경은 Fig. 5와 같이 Adjoint solver를 통해서 설계민감도를 구한 다음, 설계민감도에 따라서 Surface를 변형시키면 된다.

Fig. 5

The design sensitivity of flow uniformity

Fig. 5에서는 Adjoint solver를 통해서 산출된 설계민감도를 나타냈는데 1단계 실험계획법에서는 고려하지 못한 입구측 측면 형상에 국부적으로 설계민감도가 가장 크게 나타났다.

산출된 설계민감도를 기반으로 형상을 개선한 결과 Fig. 6과 같이 1단계 최적화에서 향상시킨 유속균일도가 0.80에서 0.85까지 추가 개선된 것을 나타났다.

Fig. 6

The shape and uniformity of 2nd optimization

물론 Iteration 과정을 거쳐서, 최종적인 최적형상을 산출하여야 하겠지만, 실제 개발과정에서는 설계민감도를 참조하여 양산성이 가능한 형상으로 모델을 재구성하여야 하기 때문에 1번의 반복계산만을 필요로 했다.

4. 결 론

실제 제품개발에서 성능을 육성하는 과정 즉, 최적화 과정은 반드시 필요하다.

본 논문에 소개된 Adjoint 최적화 방법론은 한번의 Adjoint 미분방정식의 해를 구함으로 형상개선점을 파악할 수 있어서 시간 효율적인 부분에서 기존의 전역적인 방식이나 차분적인 방법에 비해 매우 뛰어남을 알 수 있었다.

그를 통해 실제 개발과정에서 디젤엔진용 Downpipe에 형상 변경으로 유속균일도를 0.8에서 0.85까지 추가적으로 개선할 수 있었다.

하지만, 다음과 같은 사용상 제한점이 있었다.

• 1) 불연속적이거나 이산적인 수준(예: 1개, 2개)
• 2) 노이즈인자 반영 불가(Robust 최적화 불가)
• 3) 전용 S/W가 필요로 함.

하지만, 위의 사항 1), 2)은 본 논문과 같이 실험계획법을 1단계로 사전처리 함으로써 보완적으로 해결이 가능하였다.

References

• S. C. Ko, J. H. Jung, K. C. Oh, H. G. Kim, C. H. Lee, and C. B. Lee, “A study on the NOx Reduction Rate and Fluid Mixing Characteristics Depending on the Mixer Configurations in Urea-SCR Systems”, KSAE Annual Conference Proceedings, p615-620, (2011).
• J. H. Jung, J. H. Kim, and S. H. Kim, “The Inlet Shape Optimization for Aftertreatment System for Diesel Engine with Taguchi Method”, Transactions of KSAE, 20(5), p145-151, (2012). [https://doi.org/10.7467/KSAE.2012.20.5.145]
• M. B. Giles, and N. A. Pierce, “An Introduction to the Adjoint Approach to Design”, Flow, Turbulence and Combustion, 65(3-4), p393-415, (2000). [https://doi.org/10.1023/A:1011430410075]
• H. Jasak, M. Popovac, and H. Rusche, “Implementation and Numerical Stabilisation of Adjoint Flow and Turbulence Model in Open FOAM,”, 6th European Conference on Computational Fluid Dynamics (ECFD VI), (2014).
• E. Helgason, Development of Adjoint-based Optimization Methods for Ducted Flows in Vehicles, Ph. D. Dissertation, Chalmers University of Technology, Gothenburg, (2015).
• C. Othmer, “A Continuous Adjoint Formulation for the Computation of Topological and Surface Sensitivities of Ducted Flows”, International Journal for Numerical Methods in Fluids, 58(8), p861-877, (2008). [https://doi.org/10.1002/fld.1770]
• H. Jasak, Finite Volume Discretisation in OpenFOAM, Wikki Ltd, http://www.tfd.chalmers.se/~hani/kurser/OS_CFD_2015/HrvojeJasak/DiscretisationBestPractice.pdf (2012).