Tez No	İndirme	Tez Künye	Durumu
421178		A parallel monolithic approach for the numerical simulation of fluid-structure interaction problems / Akışkan-yapı etkileşimi problemlerinin sayısal simülasyonu için paralel monolitik bir yöntem Yazar:ALİ EKEN Danışman: YRD. DOÇ. DR. HAYRİ ACAR ; DOÇ. DR. MEHMET ŞAHİN Yer Bilgisi: İstanbul Teknik Üniversitesi / Fen Bilimleri Enstitüsü / Uçak ve Uzay Mühendisliği Ana Bilim Dalı Konu:Havacılık Mühendisliği = Aeronautical Engineering ; Uçak Mühendisliği = Aircraft Engineering Dizin:Akışkanlar dinamiği = Fluid dynamics ; Akışkanlar mekaniği = Fluid mechanics	Onaylandı Doktora İngilizce 2016 144 s.
Bu çalışmada akışkan-yapı etkileşimi (FSI) problemlerinin paralel tam bağlaşık bir çözüm yaklaşımıyla simülasyonuna yönelik özgün bir sayısal algoritma geliştirilmiştir. Problemin akışkan kısmi için daimi olmayan, sıkıştırılamaz Navier-Stokes denklemleri, Arbitrary Lagrangian-Eulerian (ALE) formda kullanılmıştır. Akışkan için ALE tabanlı bu hareket denklemleri, yapısal olmayan bir sonlu hacimler yaklaşımı ile ayrıklaştırılmıştır. Bu ayrıklaştırmada birincil değişkenler kenar-merkezli bir şema ile konumlandırılmıştır. Bu konumlandırmada, hız vektörü komponentleri her hücre yüzünün orta noktasında tanımlanmışken, basınç ise eleman merkezinde tanımlanmaktadır. Birincil değişkenlerin bu şekilde konumlanması kararlı bir sayısal şema oluşmasını sağlamaktadır ve basınç-hız bağlaştırması için fazladan düzenlemeler yapılmasına gerek duyulmamaktadır. Bu sonlu hacimler yaklaşımının en çekici tarafı, üniform Kartezyen çözüm ağlarında, klasik MAC (Marker and Cell) şemasında olduğu gibi, Poisson denklemi için klasik yedi-nokta Laplace operatörüne yol açmasıdır. Bu durum, çoğunlukla akışkan alt problemi hesaplama kaynaklarının büyük kısmına ihtiyaç duyduğundan, özellikle büyük ölçekli akışkan-yapı etkileşimi problemlerinin çözümünde çözüm veriminin artırılması acısından çok önemlidir. ALE tabanlı sonlu hacimler yaklaşımında, doğruluk ve kararlılık acısından dikkate alınması gereken başka bir nokta da kullanılacak çözüm ağı hareketi yöntemidir. ALE yönteminde çözüm ağı akışkan ve kati sınırları arasındaki ara yüzü takip eder ve hareket denklemleri hareketli bir çözüm ağında ayrıklaştırılır. Bu durum da zaman integrasyon şemasının doğruluğu ve kararlılığı açısından, çözüm ağı hareketine özel şartlar empoze edilmesini gerektirir. Bu şartlar, ayrık geometrik korunun yasasının (DGCL) uygulanması ile sağlanmaktadır. Bu çalışmada, mesh hareketinden kaynaklı akılar geometrik korunum yasası ayrıklaştırma seviyesinde sağlanacak şekilde hesaplanmaktadır. Süreklilik denkleminin her eleman içinde sağlanması için de özel bir dikkat sarf edilmiştir. Her elemen için süreklilik denkleminin toplamı çözüm alanı sınırlarına tam bir şekilde indirgenebilir ki bu da özellik global kütle korunumu acısından çok önemlidir. Akışkan bölgesinde zaman integrasyonu için ikinci mertebeden geri farklar formülü kullanılmaktadır. Çözüm ağı deformasyonu algoritması için verimli cebirsel bir yöntem uygulanmıştır. Bu yöntemde, akışkan iç bölgesindeki çözüm noktaları, akışkan-yapı ara yüzündeki en yakın çözüm noktalarının eksponensiyel bir deplasman fonksiyonuna göre deforme edilmektedir. Uygulanan bu cebirsel yöntemin en temel avantajı, oldukça spars cebirsel bir denkleme yol açmasıdır ki bu da bütün algoritmanın verimi acısından çok önemlidir. Yapı bölgesinin deformasyonu Saint Venant-Kirchhoff malzeme modelinin uygunluk denklemlerine dayanmaktadır. Bu yöntem yapının büyük elastik deplasman gösterdiği geometrik doğrusal olmayan problemlere uygulanabilmektedir. Yapı bölgesi hareket denklemlerinin ayrıklaştırması Lagrangian bir çerçevede klasik Galerkin sonlu hacimler yöntemine dayanmaktadır. Üç boyutlu sonlu elemanlar ayrıklaştırmasında 8-nodlu, eş-parametreli 6-yuzlu elemanlar kullanılmakta iken, iki boyutlu ayrıklaştırmalarda 4-nodlu dörtgen elemanlar kullanılmıştır. Yapısal katılık matrisi, kütle matrisi ve kuvvet vektörüne ait integraller 2-nokta Gauss tümlevi (Gauss quadrature) yöntemi ile hesaplanmaktadır. Yapı bölgesi denklemlerinin zaman integrasyonunda ise genelleştirilmiş-α yöntemi uygulanmıştır. Bu yöntem, yapının hız ve yer değiştirmeleri için Newmark tipi yaklaşımlar kullanan, tek-adim implicit, ikince mertebe bir integrasyon yöntemidir. Bu yöntemde, sayısal algoritmanın yüksek frekans sönümleme karakteri, uygun genelleştirilmiş-α parametreleri seçilerek kolaca kontrol edilebilir. Yapısal dinamik denklemlerinin doğrusallaştırılması, Newton tipi bir algoritmaya dayanmaktadır. Bu algoritmada, denklem sisteminin Jacobian matrisi her zaman adımında tam hesaplanmamaktadır. Bu yaklaşım ALE tabanlı akışkan bölgesi denklemleri için de uygulanmıştır. Aksi halde tam Newton yöntemini, tam bağlaşık denklem sisteminde fazladan sıfır olmayan bloklar oluşmasına neden olmakta ve bu da hafıza gereksinimlerini özellikle 3-boyutlu hesaplamalarda oldukça artırmaktadır. Akışkan ve yapı bölgelerine ait denklemlerin çözümü tam bağlaşık bir yaklaşıma dayanmaktadır. Bu yaklaşımda, akışkan ve yapı denklemleri tek bir denklem sistemi oluşacak şekilde inşa edilmektedir ve bu denklem sistemi her zaman adımında tam bağlaşık şekilde çözülmektedir. Akışkan ve yapı bölgeleri arasında bağlaştırma, akışkan-yapı ara yüzü boyunca birbirine uyumlu akışkan ve yapı çözüm ağları kullanılmasıyla basitleştirilmiştir. Hali hazırdaki monolitik FSI çözücü, ortaya çıkan tam bağlaşık denklem sisteminin çözümünde, ön-koşulu bir Krylov alt uzay metodu kullanmaktadır. Hızın diverjansinin sıfır olması koşulu nedeniyle oluşan sıfır blok diyagonal, bütün sitem için verimli ön-koşullandırıcıların uygulanmasını zorlaştırmaktadır. Mevcut yöntemde ise, orijinal sistemdeki sıfır blok yerine ölçeklenmiş bir ayrık Laplacian oluşturan bir üst üçgen sağ ön-koşullandırıcı uygulanmaktadır. Bu ön-koşullandırma, matris-matris çarpımları nedeniyle sıfır olmayan eleman sayısında belirgin bir artışa neden olduğu için, ön-koşullandırıcının sıfır olmayan bloğu hesaplama açısından daha az pahalı bir matris ise değiştirilmektedir. Momentum denkleminde basınç gradyanlarına olan katkı, hız vektörlerinin ayrıklaştırıldığı ortak eleman yüzeyini paylasan sağ ve soldaki elemanlardan olduğundan, kullanılan bu matris sadece bu katkılardan kaynaklanan terimleri içermektedir. Bu yaklaşım iteratif çözücünün yakınsama karakterini çok belirgin şekilde etkilemese de, özellikle 3-boyutlu hesaplamalarda, hesaplama zamanı ve hafıza gereksinimlerinde ciddi azamlar sağlamaktadır. Mevcut tek seviye iteratif çözüm yaklaşımı, sistemin simetrik olmayan doğası gereği, kısıtlı aditif Schwarz ön-koşullu esnek GMRES(m) (restricted additive Schwarz preconditioned flexible GMRES) Krylov alt uzay algoritmasına dayanmaktadır. Sistemin her alt bölgesinde blok ILU faktorizasyonu uygulanmıştır. Cebirsel denklemlerin doğasından kaynaklanan doğrusal olmama durumu nedeniyle, yeterince tatmin edici bir yakınsama kriterine ulaşıncaya kadar her zaman adımında alt-iterasyonlar uygulanmıştır. Mevcut ön-koşullu Krylov alt uzay algoritması, matris-matris çarpımları ve kısıtlı aditif Schwarz ön-koşullandırıcı uygulaması için PETSc (Portable, Extensible Toolkit for Scientific Computation) kütüphanesi kullanılmıştır. Bu kütüphane, doğrusal ve doğrusal olmayan denklem çözücülerinin paralel impilementasyonu için oluşturulmuş veri yapıları ve rutinleri içermektedir. Bütün akışkan-yapı ara yüzünün yapısal olmayan çözüm ağı, METIS kütüphanesi kullanılarak alt parçalara ayrılmıştır. Bu kütüphane, yapısal olmayan grafik ve çözüm ağlarının paralel programlamaya yönelik parçalanması için geliştirilmiş programlar içeren bir kütüphanedir. Geliştirilen bu FSI çözücünün doğruluğunu test etmek ve önerilen algoritmanın ölçeklenme karakterini incelemek amacıyla, mevcut metot literatürde sıklıkla adres edilen birçok FSI test problemine uygulanmıştır. İlk doğrulama örneği oldukça popüler 2-boyutlu bir FSI test problemidir. Problem rijit dairesel bir engel arkasına yerleştirilmiş elastik bir çubuk ile etkileşime giren Newtonian bir akıştan oluşmaktadır. Elastik çubuk akis bölgesinin alt ve üst duvarları arasında asimetrik yerleştirilmiştir ve bu akışkan yapı etkileşimi senaryosunda, engelden kopan ve ilerleyen girdapların elastik çizim üzerinde indüklediği periyodik titreşimler oluşmaktadır. Bu test probleminde hem daimi hem de daimi olmayan durumlar dikkate alınmıştır. Geliştirilen algoritmanın çözüm ağı yakınsama ve ölçeklenme karakterinin ortaya konulması amacıyla, üç fark çözüm ağı çözünürlüğü dikkate alınmıştır. Çözüm ağı çözünürlüğünün, işlemci sayısının, ILU(k) ön-koşullandırıcı seviyesinin ve kısıtlı aditif Schwarz ön-koşullandırıcıdaki üst üste binme miktarının performansa etkileri ortaya konulmuştur. Bu test probleminde elde edilen sonuçlar literatürdeki birçok çözümle karşılaştırılmış ve mevcut çözücünün doğruluğu ispat edilmiştir. İkinci FSI test probleminde yine literatürde sıklıkla adres edilen 3-boyutlu bir konfigürasyon dikkate alınmıştır. Bu problem konfigürasyonu basitçe elastik arterlerden kan akisini simüle etmektedir. Problemde sıkıştırılamaz viskoz bir akis esnek dairesel bir tüp ile çevrelenmiştir ve empoze edilen sinir koşulları ile elastik tüpte zamana bağlı radyal ve eksenel deplasmanlar oluşturacak şekilde ilerleyen bir dalga çözümüne ulaşılmaktadır. Bu 3-boyutlu test problemi için de bir ölçeklenme testi yapılmıştır. Bu test için iki farklı çözüm ağı çözünürlüğü dikkate alin mistir. Radyal deplasmanlar için hesaplanan sonuçların literatür ile iyi uyum içinde olduğu gösterilmiştir. Takip eden 3-boyutlu FSI test problemlerinden bir diğeri dikdörtgensel bir kanal içine yerleştirilmiş elastik bir cisimden oluşmakta ve daimi bir akis çözümü vermektedir. Cisim üzerinde bir kontrol noktasının deplasman çözümü verilmiş ve literatür ile olan uyumu gösterilmiştir. Dördüncü test problemi, dikdörtgensel bir engel arkasına yerleştirilmiş bir bayrağın girdap indüklü titreşimlerini modelleyen 3-boyutlu bir konfigürasyondur. Bu problem için mevcut hesaplamalar oldukça yüksek çözünürlüklü bir çözüm ağında yapılmıştır. Bu test problem bilgisayar gücü açısından oldukça zorlayıcı olsa da mevcut FSI algoritması belirgin bir performans kaybı yasamadan benzer ölçeklenme özellikleri göstermiştir. Son test problemi, geliştirilen FSI algoritmasının, özellikle akışkan bölgesi yapı bölgesi ile tamamen çevrelendiği durumlarda, kütle korunum kabiliyetinin ortaya konulması amacıyla tasarlanmıştır. Bu amaçla tasarlanan problem, paralel rijit duvarlar arasına simetrik olarak yerleştirilmiş iki boyutlu dairesel elastik bir yüzükten oluşmaktadır. Bu dairesel yapı bir akışkan ile çevrelenmiş olmakla birlikte kendisi de bir akışkan bölgesini çevrelemektedir. Bu konfigürasyon basitçe kırmızı kan hücrelerinin düşük Reynolds sayılı akışta deformasyonunu simüle etmektedir. Bu son test problemi, geliştirilen algoritmanın, kullanılan uyumlu FSI ara yüz şartı ile makine hassasiyetinde kütle korunumunu sağladığını göstermektedir. FSI çözücüsünün doğrulanması amacıyla yapılan sayısal deneylerden sonra, mevcut algoritma kardiyovaskular akışkan-yapı etkileşiminde sıklıkla karşılaşın gerçekçi bir problemin çözümünde kullanılmıştır. Bu akışkan-yapı etkileşimi problemi, damar çatallanma tepesinde anevrizma ihtiva eden bir beyin arterinde impulsif olarak başlatılan bir FSI problemidir. Kan Newtonian bir akışkan olarak, damar duvarı ise Saint Venant Kirchhoff modeline uygun bir elastik malzeme olarak modellenmiştir. Başlangıçta kullanılan tamamen hegzahedral konformal çözüm ağı oktree metodu kullanılarak oluşturulmuştur. Akışkan hız alanı, kan basıncı, duvar kayma gerilmeleri gibi birçok hemodinamik büyüklüğün yanında, zamana bağlı damar duvar deplasmanları hesaplanmıştır. İyi ölçeklenme karakteri bu problemde de elde edilmiştir. Bu çalışmada geliştirilen FSI çözücüsünün geliştirme ve test aşamalarında kullanılan yöntemler özetlenmiştir. Mevcut yöntemin avantaj ve dezavantajlarına değinilmiş ve muhtemel gelecek uygulamalarından bahsedilmiştir.
In this study, a novel numerical algorithm has been developed for the simulation of fluid-structure interaction problems in a parallel fully coupled solution approach. For the fluid part of the problem, the Arbitrary Lagrangian-Eulerian (ALE) form of the incompressible, unsteady Navier-Stokes equations are used to model the fluid motion. The ALE based governing equations of the fluid domain are discretized using an unstructured finite volume approach, where the primitive variables of the fluid are defined based on a side-centered arrangement. In this side-centered arrangement, the velocity vector components are defined at the mid-point of each cell face, while the pressure is defined at the element centroid. This arrangement of the primitive variables leads to a stable numerical scheme and it does not require any \textit{ad-hoc} modifications in order to enhance the pressure-velocity coupling. The most appealing feature of the present finite volume approach is that it leads to the classical seven-point Laplace operator on uniform Cartesian meshes for the pressure Poisson equation as in the classical MAC scheme, which is very important for the efficient solution of the large-scale FSI problems in 3D, since the fluid subproblem requires the majority of the computational resources in the coupled system for the most of cases. Another point to be considered in terms of accuracy and stability of the ALE based finite volume approach is the mesh movement strategy. In the ALE method, the mesh follows the interface between the fluid and solid boundary and the governing equations are discretized on a moving mesh, which requires the imposition of special conditions on the mesh movement for the accuracy and stability of the time integration scheme. This condition is imposed by the enforcement of the discrete geometric conservation law (DGCL). In the present study, the fluxes due to mesh motion are computed in a way that the geometric conservation law is satisfied in a discrete level. A special attention is also given in order to satisfy the continuity equation exactly within each element. The summation of the discrete continuity equations can be exactly reduced to the domain boundary, which is important for the global mass conservation. For time integration in the fluid domain second-order backward difference formula is used. The mesh deformation algorithm within the fluid domain is based on an efficient algebraic method, where the interior fluid nodes are deformed based on a exponential function of the displacement of the nearest node on the common fluid-structure interface. The main advantage of the present algebraic method is that it leads to a highly sparse algebraic equation, which is very important for the efficiency of the overall algorithm. The deformation of the solid domain is governed by the constitutive laws for the Saint Venant-Kirchhoff material model, which is applicable to geometrically nonlinear analysis of the structure with large elastic displacements. The discretization of the governing equations of the solid domain is accomplished by the classical Galerkin finite element method in a Lagrangian frame. For three-dimensional finite element discretization, eight-node iso-parametric hexahedral elements are used to model the solid domain, while four-node quadrilateral finite elements are used for two-dimensional discretizations. The numerical integrals related to the structural stiffness matrix, mass matrix and force vector are calculated using two-point Gauss quadrature method. The time integration of equations for the solid domain is based on the generalized$-\alpha$ method, which is a second order accurate, one-step implicit time integration method based on Newmark type approximations for the displacements and velocities of the structure. In this method, the high frequency damping character of the numerical algorithm can be easily controlled by selecting the appropriate generalized$-\alpha$ parameters. The linearization of the structural dynamics equations is based on a quasi-Newton type algorithm, where the Jacobian matrix of the equation system is not exactly calculated at each time level. This approach is also applied for the ALE algorithm for the fluid domain, where the exact Newtons' method leads to extra non-zero blocks in the whole coupled system of equations, which significantly increases the memory demand in three-dimensional computations. The solution of the equations resulting from the discreatization of the fluid and structure domains is based on a monolithic approach, where the fluid and solid equations are assembled into a single equation system and solved in a fully coupled way at each time level. The coupling between the fluid and structure domains is simplified by creating fluid and structure meshes which are conformal along the fluid-structure interface. The present monolithic FSI solver uses a preconditioned Krylov subspace method to solve the resulting fully coupled system. Because of the zero-block diagonal resulting from the divergence-free constraint, it is rather difficult to construct robust preconditioners for the whole coupled system. In the present approach, we use an upper triangular right preconditioner which results in a scaled discrete Laplacian instead of the zero block in the original system. Since this preconditioning leads to a significant increase in the number of non-zero elements due to the matrix-matrix multiplications, we replace the non-zero block of the preconditioner with a computationally less expensive matrix, which contains the contributions from the right and left element only, since the largest contribution for the pressure gradients in the momentum equations comes from the right and left elements that share the common face where the components of the velocity vector are discretized. Although, this approximation does not change the convergence rate of the iterative solver significantly, it leads to a significant reduction in the computing time and memory requirement in 3D. The present one-level iterative solution approach is based on the restricted additive Schwarz method with the flexible GMRES(m) Krylov subspace algorithm due to the non-symmetric nature of the system. A block-incomplete LU (ILU) factorization is used within each partitioned sub-domains. In order to handle the nonlinear nature of the algebraic equations, sub-iterations are performed for each time step until a satisfactory convergence criteria is reached. The implementation of the present preconditioned Krylov subspace algorithm, matrix-matrix multiplications and the restricted additive Schwarz preconditioner are carried out using the PETSc (Portable, Extensible Toolkit for Scientific Computation) library developed at the Argonne National Laboratories, which is a collection of data structures and routines for parallel implementation of linear and nonlinear equation solvers. The unstructured computational mesh for the whole fluid and solid domains is decomposed into a set of sub-domains using the METIS library, which is a set of programs developed for partitioning unstructured graphs, meshes, and producing fill reducing orderings for sparse matrices. In order to validate the accuracy of the developed FSI solver and examine the scaling properties of the proposed algorithm, the present method is initially applied to several FSI benchmark problems that are frequently addressed in the literature. The first validation case is a rather popular two-dimensional FSI test problem. The problem consists of a Newtonian fluid flow interacting with an elastic bar behind a fixed rigid circular obstacle which is placed asymmetrically between parallel lateral walls. In this fluid-structure interaction scenario, shedding vortices induce time dependent displacements and periodic oscillations of the structure. Both steady and unsteady flow solutions are considered for this test case. In order to demonstrate the mesh convergence and scaling properties of the present FSI algorithm, three different mesh resolutions are considered. The effects of the mesh resolution, the number of processors, the level of ILU(k) preconditioner and the amount of overlap for the restricted additive Schwarz preconditioner are presented. The obtained results for this two-dimensional FSI problem are compared to several solutions from the literature and the accuracy of the present algorithm is validated. In the second FSI benchmark case a three-dimensional problem configuration is considered, which is also a well addressed FSI problem in the literature. This problem setting simply simulates the blood flow through elastic arteries. The test configuration consists of an incompressible viscous flow through a flexible circular tube. The considered boundary conditions for this FSI problem produces a wave propagation through the elastic tube with resultant time dependent radial and axial deformations of the structure. A scaling test is also conducted for this 3D benchmark problem for two different mesh resolutions. It is shown that the computed results for radial displacements of the structure are relatively in good agrement with the results in the literature. In the subsequent validation case, a steady FSI benchmark is considered, where an elastic solid is immersed in a rectangular channel, and the deflection of a test point on the elastic structure is taken in account for comparison to the result given in the literature. The computed results for this FSI test case are also found to be in good agrement with the literature. The fourth test case represents a three-dimensional configuration for the vortex-induced vibrations of a flag attached behind a rectangular rigid body in an incompressible viscous flow. The present calculations for this test case are performed on a relatively fine mesh. Although this test case is rather demanding in terms of the required computer power, the FSI algorithm achieves similar scaling without a significant performance loss. The final benchmark problem is motivated to demonstrate the exact mass conservation property of the present FSI algorithm when the fluid domain is fully enclosed in a solid domain. For this purpose, a two dimensional flexible circular ring is placed symmetrically between two parallel rigid plates. This circular ring, which is immersed in a fluid domain, completely encloses a part of the fluid domain itself. This configuration mimics the deformation of a red blood cell inside a capillary wall at a very low Reynolds number. This final test demonstrates that the developed algorithm with its compatible FSI interface condition achieves mass conservation at machine precision. Subsequent to the numerical experiments and FSI solver validation, the present algorithm is applied to a realistic fluid-structure interaction problem frequently encountered in cardiovascular FSI. This fluid-structure interaction problem corresponds to an impulsively accelerated blood flow within a cerebral artery with aneurysm located at the apex of the bifurcation of a branch. The blood is assumed to be a Newtonian fluid, and the artery wall is modeled as a Saint Venant Kirchhoff material. The octree method is used to generate the initial all hexahedral conforming coarse mesh. Various hemodynamic quantities of interest like fluid velocities, blood pressure and wall shear stresses (WSS) are computed as well as the time dependent artery wall deformations. Good scaling character has also been obtained for this FSI application. Finally, the methods used to develop and test the present FSI solver are summarized. Advantages and the drawbacks of the present solver are addressed with possible future applications.