A high-order finite-volume solver for supersonic flows

Abstract

Nowadays, Computational Fluid Dynamics (CFD) is a powerful tool in engineering used in various industries such as automotive, aerospace and nuclear power. More than ever the growing computational power of modern computer systems allows for realistic modelization of physics. Most of the open-source codes, however, offer a second-order approximation of the physical model in both space and time. The goal of this thesis is to extend this order of approximation to what is defined as high-order discretization in both space and time by developing a two-dimensional finite-volume solver. This is especially challenging when modeling supersonic flows, which shall be addressed in this study. To tackle this task, we employed the numerical methods described in the following. Curvilinear meshes are utilized since an accurate representation of the domain and its boundaries, i.e. the object under investigation, are required. High-order approximation in space is guaranteed by a Central Essentially Non-Oscillatory (CENO) scheme, which combines a piece-wise linear reconstruction and a k-exact reconstruction in region with and without discontinuities, respectively. The usage of multi-step methods such as Runge-Kutta methods allow for a high-order approximation in time. The algorithm to evaluate convective fluxes is based on the family of Advection Upstream Splitting (AUSM) schemes, which use an upwind reconstruction. A central stencil is used to evaluate viscous fluxes instead. When using high-order schemes, discontinuities induce numerical problems, such as oscillations in the solution. To avoid the oscillations, the CENO scheme reverts to a piece-wise linear reconstruction in regions with discontinuities. However, this introduces a loss of accuracy. The CENO algorithm is capable of confining this loss of accuracy to the cells closest to the discontinuity. In order to reduce this accuracy loss Adaptive Mesh Refinement (AMR) is used. This algorithm refines the mesh near the discontinuity, confining the loss of accuracy to a smaller portion of the domain. In this study, a combination of the CENO scheme and the AUSM schemes is used to model several problems in different compressibility regimes, with a focus on supersonic flows. The scope of this thesis is to analyze the capabilities and the limitations of the proposed combination. In comparison to traditional implementations, which can be found in literature, our implementation does not impose a limit on the refinement ratio of neighboring cells while utilizing AMR. Due to the high computational expenses of a high-order scheme in conjunction with AMR, our solver benefits from a shared memory parallelization. Another advantage over traditional implementations is that our solver requires one layer of ghost cells less for the transfer of information between adjacent blocks. The validation of the solver is performed in different steps. We assess the order of accuracy of the CENO scheme by interpolating a smooth function, in this case the spherical cosine function. Then we validate the algorithm to compute the inviscid fluxes by modeling a Sod shock tube. Finally, the Boundary Conditions (BCs) for the inviscid solver and its order of accuracy are validated by modeling a convected vortex in a supersonic uniform flow. The curvilinear mesh is validated by modeling the flow around a NACA0012 airfoil. The computation of the viscous fluxes is validated by modeling a viscous boundary layer developing on a flat plate. The BCs for viscous flows and the curvilinear implementation are validated by modeling the flow around a cylinder and a NACA0012 airfoil. The AUSM schemes are tested for shock robustness by modeling an inviscid hypersonic cylinder at a Mach number of 20 and a viscous hypersonic cylinder at a Mach number of 8.03. Then, we validate our AMR implementation by modeling a two-dimensional Riemann problem. All the validation results agree well with either numerical or experimental results available in literature. The performance of the code, in terms of computational time required by the different orders of approximation and the parallel efficiency, is assessed. For the former a supersonic vortex convection served as an example, while the latter used a two-dimensional Riemann problem. We obtained a linear speed-up until 12 cores. The highest speedup value obtained is 20 with 32 cores. Furthermore, the solver is used to model three different supersonic applications: the interaction between a vortex and a normal shock, the double Mach reflection and the diffraction of a shock on a wedge. The first application resembles a strong interaction between a vortex and a steady shock wave for two different vortex strengths. In both cases our results perfectly match the ones obtained by a Weighted Essentially Non-Oscillatory (WENO) scheme documented in literature. Both schemes are approximating the solution with the same order of accuracy in both, time and space. The second application, the double Mach reflection, is a challenging problem for high-order solvers because the shock and its reflections interact strongly. For this application, all AUSM-schemes under investigation fail to obtain a stable result. The main form of instability encountered is the Carbuncle phenomenon. Our implementation overcomes this problem by combining the AUSM+M scheme with the formulation of the speed of sound of the AUSM+up scheme. This combination is capable of modeling this problem without instabilities. Our results are in agreement with those obtained with a WENO scheme. Both, the reference solutions and our results, use the same order of accuracy in both, time and space. Finally, the third example is the diffraction of a shock past a delta wedge. In this configuration the shock is diffracted and forms three different main structures: two triple points, a vortex at the trailing edge of the wedge and a reflected shock traveling upwards. Our results agree well with both, numerical and experimental results available in literature. Here, a formation of a vortex-let is observed along the vortex slip-line. This vorticity generation under inviscid flow condition is studied and we conclude that the stretching of vorticity due to compressibility is the reason. The same formation is observed when the angle of attack of the wedge is increased in the range of 0-30. In general, the AUSM+up2 scheme performed best in terms of accuracy for all problems tested here. However, for configurations, in which the Carbuncle phenomenon may appear, the combination of the AUSM+M scheme and the computation of the speed of sound formula of the AUSM+up scheme is preferable for stability reasons. During our computations, we observe a small undershooting right behind shocks on curved boundaries. This is imputable to the curvilinear approximation of the boundaries, which is only second-order accurate. Our experience shows that the smoothness indicator formula in its original version, fails to label uniform flow regions as smooth. We solve the issue by introducing a threshold for the numerator of the formula. When the numerator is lower than the threshold, the cell is labeled as smooth. A value higher than 10^-7 for the threshold might force the solver to apply high-order reconstruction across shocks, and therefore will not apply the piece-wise linear reconstruction which prevents oscillations. We observe that the CENO scheme might cause unphysical states in both inviscid and viscous regime. By reconstructing the conservative variables instead of the primitive ones, we are able to prevent unphysical states for inviscid flows. For the viscous flows, temporarily reverting to first-order reconstruction in the cells where the temperature is computed as negative, prevents unphysical states. This technique is solely required during the first iterations of the solver, when the flow is started impulsively. In this study the CENO, the AUSM and the AMR methods are combined and applied successfully to supersonic problems. When modeling supersonic flow with high-order accuracy in space, one should prefer the combination of the AUSM schemes and the CENO scheme. While the CENO scheme is simpler than the WENO scheme used in comparison, we show that it yields results of comparable accuracy. Although it was beyond the scope of this study, the AUSM can be extended to real gas modeling which constitutes another advantage of this approach.

Hesaplamalı Akışkanlar Dinamiği günümüzde; arabalar, uçaklar, motorlar, nükleer reaktörler vb. gibi her türlü kullanım için yeni araçlar tasarlamakta kullanılan güçlü bir araçtır. Modern iş istasyonlarının sürekli artan gücü, her zamankinden daha gerçekçi fizik modellemesine olanak tanımaktadır. Çoğunlukla açık kaynak kodlarının hem uzayda hem de zamanda modellenen fiziğin ikinci mertebeden bir yaklaşımını sunduğu fark edilebilir. Bu tezin amacı, iki boyutlu bir sonlu hacim çözücüsü geliştirerek bu yaklaşım mertebesini hem uzay hem de zamanda yüksek mertebe ayrıklaştırmaya genişletmektir. Bu yöntem özellikle süpersonik akışlar modellenmek istendiğinde zordur. İlk adım, hesaplama alanı ve araştırılan nesnenin uygun şekilde oluşturulmasıdır, bu nedenle eğrisel bir ağ üzerinde çalışan bir çözücü şarttır. Uzayda yüksek mertebe, süreksizliklerin olduğu bölgede k-tam yeniden oluşturma (reconstruction) ve parçalı doğrusal (piece-wise linear) yeniden oluşturmayı birleştiren bir Merkezi Esasen Salınımsız (Central Essentially Non-Oscillatory CENO) şema tarafından güvence altına alınmaktadır. Runge-Kutta yöntemi gibi çok adımlı bir yöntemin kullanılması, çözücünün zamanda yüksek mertebe yaklaşıma sahip olmasını sağlamaktadır. Konvektif akıları hesaplamak için kullanılan algoritma, akıları yeniden oluşturmak için akım önü şemasını (upwind scheme) kullanan Adveksiyon Yukarı-akış Bölme Yöntemi (Adevction Upstream Splitting Method AUSM) şema ailelerini temel almaktadır. Viskoz akıları hesaplamak içinse merkezi bir şablon kullanılmıştır. Süreksizliklerin olduğu bölgelerde doğruluk kaybını azaltmak için, bu bölgelerdeki ağı iyileştirmeye yönelik bir adaptif ağ iyileştirme (Adaptive Mesh Refinement AMR) algoritması kullanılmıştır. AMR ve CENO şemalarının kombinasyonu çözücünün, süreksizliklerin olmadığı bölgede yüksek mertebeden yaklaşımı korumasına ve süreksizliklerin olduğu bölgede birinci mertebeden yaklaşımı sınırlandırmasına izin verir. AMR uygulamamız birbiriyle temas halinde olan bölgeler arasındaki iyileştirme oranına bir sınır getirmemektedir. Çözücü ayrıca, bloklar arasında bilgi aktarmak için bir tane daha az hayalet hücre katmanı kullanmaktadır. Bu tezin kapsamı, süpersonik akışlara odaklanarak farklı sıkıştırılabilirlik rejimlerinde çeşitli problemleri simüle ederken CENO şemasının kabiliyetlerini ve sınırlarını analiz etmektir. Çözücünün doğrulanması farklı adımlarda gerçekleştirilmiştir. CENO şemasının doğruluk mertebesi, düzgün bir fonksiyon olan küresel kosinüs fonksiyonunun interpolasyonu yapılarak değerlendirilmiştir. Ardından, bir Sod şok tüpü simüle edilerek viskoz olmayan akıların hesaplandığı algoritma doğrulanmıştır. Son olarak, viskoz olmayan çözücü için sınır koşulları ve doğruluk mertebesi, süpersonik üniform bir akış içerisinde taşınan bir girdap simüle edilerek değerlendirilmiştir. Eğrisel ağ, bir NACA0012 kanat profili etrafındaki akışı simüle ederek doğrulanmıştır. Viskoz akıların hesaplanması, düz bir plaka üzerinde gelişen viskoz bir sınır tabaka simüle edilerek doğrulanmıştır. Viskoz akışların sınır koşulları ve çözücünün eğrisel ağ üzerindeki performansı, bir silindir ve bir NACA0012 kanat profili etrafındaki akışların modellenmesiyle doğrulanmıştır. AUSM şemaları, bir silindir etrafında 20 Mach ve 8.03 Mach sayılarındaki hipersonik viskoz akışlar modellenerek, çözücünün geniş bir akış koşulu aralığında şokları yakalama kabiliyeti bakımından test edilmiştir. Son olarak, iki boyutlu bir Riemann problemi simüle edilerek AMR uygulaması doğrulanmıştır. Tüm doğrulama sonuçları, literatürde mevcut olan sayısal veya deneysel sonuçlarla karşılaştırılmıştır. Kodun performansı, farklı yaklaşım mertebelerinin gerektirdiği hesaplama süreleri ve paralel verimlilik açısından değerlendirilmiştir. Birincisi, süpersonik bir girdap taşınımını simüle ederek hesaplanırken, ikincisi iki boyutlu bir Riemann problemini simüle ederek hesaplanmıştır. 12 çekirdeğe kadar doğrusal bir hızlanma elde edilmiştir. En yüksek hızlanma değeri 20 olup 32 çekirdek kullanılarak elde edilmiştir. Sonuncu fakat bir o kadar da önemli olarak, çözücü üç farklı uygulamayı modellemek için kullanılmıştır. Tüm uygulamalar süpersonik akışlara odaklanmaktadır. Bunlardan ilki, iki farklı durum için bir girdap ile durağan şok dalgası arasındaki güçlü etkileşimdir. İki durum için de elde ettiğimiz sonuçlar literatürde Ağırlıklı ENO (Weighted Essentially Non-Oscillatory WENO) şeması kullanılarak elde edilen mevcut sonuçlarla iyi bir uyum içindedir. İki şemanın da yaklaşık çözümü uzayda ve zamanda aynı doğruluk mertebesindedir. İkinci uygulama yüksek mertebe çözücüler için zorlu bir problem olan çift Mach yansımasıdır (double Mach reflection). Bu uygulama için incelenen tüm AUSM şemaları, kararsızlık içermeyen bir sonuç elde etmekte başarısız olmuştur. Karşılaşılan ana kararsızlık oluşumu Carbuncle fenomenidir. Bu durum için bizim çözümümüz, AUSM+M şemasını, AUSM+up şemasının ses hızı formülasyonu ile birleştirmektir. Bu kombinasyon, problemi kararsızlıklar olmadan modelleme yeteneğine sahiptir. Elde ettiğimiz sonuçlar literatürde WENO şeması ile elde edilen mevcut sonuçlarla iyi bir uyum içindedir. Hem referans hem de sonuçlarımız, zaman ve uzayda aynı doğruluk mertebesini kullanmaktadır. Son uygulama, bir delta kamasını geçen şokun kırınımıdır. Sonuçlarımız literatürde mevcut olan hem sayısal hem de deneysel sonuçlarla uyumludur. Bu problemde girdap kayma hattı boyunca girdapçık oluşumu gözlemlenmektedir. Viskoz olmayan akış koşulu altında gözlemlenen bu girdaplılık (vorticity) oluşumunu inceledik ve bu oluşumun sebebinin sıkıştırılabilirlik nedeniyle girdap gerilmesi olduğu sonucuna vardık. Aynı oluşum, kamanın hücum açısı 0-30 aralığında arttırıldığında da gözlemlenmiştir. Genel olarak, deneyimlerimiz AUSM+up2 şemasının burada simüle edilen tüm problemler için en iyi doğruluğa sahip olduğunu söylemektedir, ancak Carbuncle fenomeninin ortaya çıkabileceği konfigürasyon için, AUSM+M şeması ve AUSM+up şemasının ses hızı formülasyonunun kombinasyonu tercih edilebilir. Hesaplamalarımız sırasında, kavisli sınırlarda şokların hemen arkasında küçük bir düşüş gözlemliyoruz. Bu durum, sınırların eğrisel yaklaşımına atfedilebilir. Buradaki çözüm, kodu spline veya eşdeğer bir teknik aracılığıyla sınırların yüksek mertebe temsiliyle çalışacak şekilde genişletmektir. Deneyimlerimiz, orijinal versiyonundaki süreklilik göstergesi formülünün, üniform akış bölgelerini düzgün olarak etiketlemediğini göstermektedir. Bu sorunu formülün payı için bir eşik değeri getirerek çözdük. Pay, bir hücrenin eşik değerinden düşük olduğunda, hücre sürekli ismiyle etiketlenmektedir. Eşik için 10^-7'dan yüksek bir değer, çözücüyü şok boyunca k-tam yeniden oluşturma uygulamaya zorlayabilir ve bariz nedenlerden dolayı bundan kaçınılmalıdır. CENO şemasının hem viskoz hem de viskoz olmayan rejimde fiziksel olmayan durumlara neden olabileceğini gözlemledik. Birincil (primitive) değişkenler yerine korunumlu (conservative) değişkenleri yeniden oluşturarak viskoz olmayan akışlar için fiziksel olmayan durumu önleyebiliyoruz. Viskoz akışlar için çözümümüz, sıcaklığın negatif olarak hesaplandığı hücrelerde geçici olarak birinci mertebe çözüme dönmektir. Bu yalnızca, akış aniden başlatıldığında çözücünün ilk iterasyonları sırasında gerçekleşmektedir. Sonuç olarak, AUSM şemaları ile CENO şemasının kombinasyonunun süpersonik akışları simüle etmek için iyi bir aday olduğunu söyleyebiliriz.