ISTANBUL TECHNICAL UNIVERSITY ⋆ GRADUATE SCHOOL OF SCIENCE ENGINEERING AND TECHNOLOGY
A PARALLEL MONOLITHIC APPROACH FOR THE NUMERICAL SIMULATION OF FLUID-STRUCTURE INTERACTION PROBLEMS
Ph.D. THESIS Ali EKEN
Department of Aeronautics and Astronautics Engineering Aeronautics and Astronautics Engineering Program
ISTANBUL TECHNICAL UNIVERSITY ⋆ GRADUATE SCHOOL OF SCIENCE ENGINEERING AND TECHNOLOGY
A PARALLEL MONOLITHIC APPROACH FOR THE NUMERICAL SIMULATION OF FLUID-STRUCTURE INTERACTION PROBLEMS
Ph.D. THESIS Ali EKEN (511082107)
Department of Aeronautics and Astronautics Engineering Aeronautics and Astronautics Engineering Program
Thesis Advisor: Assist. Prof. Dr. Hayri ACAR Co-Advisor: Assoc. Prof. Dr. Mehmet ¸SAH˙IN
˙ISTANBUL TEKN˙IK ÜN˙IVERS˙ITES˙I ⋆ FEN B˙IL˙IMLER˙I ENST˙ITÜSÜ
AKI ¸SKAN-YAPI ETK˙ILE ¸S˙IM˙I PROBLEMLER˙IN˙IN SAYISAL S˙IMÜLASYONU ˙IÇ˙IN
PARALEL MONOL˙IT˙IK B˙IR YÖNTEM
DOKTORA TEZ˙I Ali EKEN (511082107)
Uçak ve Uzay Mühendisli˘gi Anabilim Dalı Uçak ve Uzay Mühendisli˘gi Programı
Tez Danı¸smanı: Assist. Prof. Dr. Hayri ACAR E¸s Danı¸sman: Assoc. Prof. Dr. Mehmet ¸SAH˙IN
Ali EKEN, a Ph.D. student of ITU Graduate School of Science Engineering and Tech-nology 511082107 successfully defended the thesis entitled “A PARALLEL MONO-LITHIC APPROACH FOR THE NUMERICAL SIMULATION OF FLUID-STRUCTURE INTERACTION PROBLEMS”, which he prepared after fulfilling the requirements specified in the associated legislations, before the jury whose signatures are below.
Thesis Advisor : Assist. Prof. Dr. Hayri ACAR ... Istanbul Technical University
Co-advisor : Assoc. Prof. Dr. Mehmet ¸SAH˙IN ... Istanbul Technical University
Jury Members : Prof. Dr. Fırat O. ED˙IS ... Istanbul Technical University
Prof. Dr. Muammer KALYON ... Istanbul Commerce University
Prof. Dr. N. L. Ok¸san ÇET˙INER YILDIRIM ... Istanbul Technical University
Assoc. Prof. Dr. Metin MURADO ˘GLU ... Koç University
Assoc. Prof. Dr. Kerem PEKKAN ... Koç University
Date of Submission : 26 November 2015 Date of Defense : 14 January 2016
FOREWORD
I would like to thank my thesis advisors Assoc. Prof. Dr. Mehmet ¸Sahin and Assist. Prof Dr. Hayri Acar for their guidance and endless support throughout this thesis study. I also wish to express my gratitude to the jury members for sharing their scientific background and for their constructive comments on this dissertation.
I would also like to acknowledge the financial support of the Turkish National Scientific and Technical Research Council (TUBITAK) through project number 112M107. This acknowledgement extends to the Faculty of Aeronautics and Astronautics at ITU for the use of the Chimera machine, to the National Center for High Performance Computing of Turkey (UYBHM) and to TUBITAK ULAKBIM High Performance and Grid Computing Center for the use of their computer resources. I wish to thank my friends at Trisonic Lab. Last but not least; I would like to thank my family for their endless support throughout my life.
January 2016 Ali EKEN
TABLE OF CONTENTS Page FOREWORD... ix TABLE OF CONTENTS... xi ABBREVIATIONS ... xiii SYMBOLS ... xv
LIST OF TABLES ...xvii
LIST OF FIGURES ... xix
SUMMARY ... xxi
ÖZET ... xxv
1. INTRODUCTION ... 1
1.1 Fluid Structure Interaction and Applications ... 1
1.2 Literature Review ... 3
1.3 Outline ... 7
2. MATHEMATICAL MODELLING ... 9
2.1 Kinematics of a Deformable Body ... 9
2.2 Lagrangian, Eulerian and ALE Descriptions... 11
2.3 Structure ... 14
2.3.1 Deformation gradient... 14
2.3.2 Volume change and compresibility... 15
2.3.3 Strain... 16 2.3.4 Stress measures... 17 2.3.5 Constitutive equation ... 19 2.4 Fluid... 21 2.4.1 Kinematics ... 21 2.4.2 Constitutive equation ... 22
2.4.3 Conservation law and governing equations ... 23
2.4.3.1 Continuity equation ... 24
2.4.3.2 Momentum equations ... 25
2.4.3.3 ALE form of Navier-Stokes equations ... 26
3. DISCRETIZATION OF GOVERNING EQUATIONS... 27
3.1 Discretization of Structure Equations... 27
3.1.1 Weak formulation and finite element discretization ... 28
3.1.1.1 Stiffness matrix ... 29
3.1.1.2 Mass matrix ... 32
3.1.1.3 Surface traction ... 33
3.1.2 Time integration... 34
3.1.2.1 Generalized−α method ... 34
3.2 Discretization of Fluid Equations ... 36
3.2.1 Geometric conservation law and swept volume calculations ... 37
3.2.1.1 DGCL for first-order time integration ... 38
3.2.1.2 DGCL for second-order time integration ... 38
3.2.2 Three-dimensional numerical discretization ... 39
3.2.3 Mesh deformation algorithm ... 44
4. FLUID-STRUCTURE COUPLING AND THE MONOLITHIC SOLVER . 47 4.1 Interface Conditions ... 47
4.2 Discretization of fluid-structure interface equations ... 47
4.3 Coupled Equations and Monolithic Solver... 49
5. FSI BENCHMARKS ... 53
5.1 Test Case I: An Elastic Bar Behind a Rigid Cylinder ... 53
5.2 Test Case II: A Pressure-pulse Propagating in a Flexible Tube ... 59
5.3 Test Case III: A Three-Dimensional Elastic Solid in a Channel Flow ... 64
5.4 Test Case IV: Vortex-induced Vibrations of a Flag ... 67
5.5 Test Case V: Elastic Circular Ring Moving in a Channel ... 71
6. UNSTEADY BLOOD FLOW IN A CEREBRAL ARTERY WITH ANEURYSM... 77
6.1 Computational Methods in Hemodynamics ... 77
6.2 Unsteady Blood Flow In A Cerebral Artery With Aneurysm ... 79
7. CONCLUSIONS AND RECOMMENDATIONS ON FUTURE WORK ... 93
REFERENCES... 97
APPENDICES ... 105
APPENDIX A.1 ... 107
ABBREVIATIONS
ALE : Arbitrary Lagrangian-Eulerian AMG : Algebraic Multi-Grid
ASM : Additive Schwarz Method
BDF2 : Backward Difference Formula-Second order CFD : Computational Fluid Dynamics
CSD : Computational Structural Dynamics DOF : Degrees of Freedom
DGCL : Discrete Geometric Conservation Law FGMRES : Flexible Generalized Minimal Residuals FSI : Fluid-Structure Interaction
GCL : Geometric Conservation Law
GHz : GigaHertz
GMRES : Generalized Minimal Residuals ILU : Incomplete Lower-Upper LCS : Lagrangian Coherent Structures
M : Mesh, i.e. M1, M2
MAC : Marker and Cell MAV : Micro Air Vehicles
PETSc : Parallel Extensible Toolkit for Scientific Computing RBF : Radial Basis Functions
WSS : Wall Shear Stress
SYMBOLS
Ani : Area vectors of the textitith triangular face b : Body forces
C : Cauchy-Green deformation tensor d : Displacement vector
E : Green-Lagrange strain tensor EI : Lagrangian cartesian basis
F : Deformation gradient tensor
J : Deformation gradient determinant (Jacobian) K, M, F : Stiffness matrix, mass matrix, load vector n, N : Normal vector
P : First Piola-Kirchhoff stress tensor p : Fluid pressure
S : Second Piola-Kirchhoff stress tensor t, T : Tangent vector
t : Time
ux, uy, uz : Components of the velocity vector in x, y and z directions
V : Volume
Vn : Total dual control volumes at time levels n v(x) : Test function
X : Lagrangian position vector xi, yi, zi : Global coordinates of node i β,γ : Newark parameters
ε : Strain rate tensor
ξ,η,ζ : Natural coordinates
λ : Second viscosity coefficient
µ : Dynamic viscosity
ρ : Solid material density
ρ0 : Solid material density per unit undeformed volume
σs : Cauchy stress tensor
τ : Shear stress tensor
φ(x,t) : Mapping function Ω0 : Undeformed domain
Ω : Domain, solid or fluid ()T : Transpose
LIST OF TABLES
Page
Table 5.1 : Parameter settings for test case I ... 54
Table 5.2 : Displacements and forces for FSI-1(Re = 20) of test case I... 55
Table 5.3 : Displacements, forces and frequencies for FSI-3(Re = 200) of test case I. ... 58
Table 5.4 : The convergence and scaling properties for FSI-3(Re = 200) (Intel Xeon 5550, 2.67GHz). ... 60
Table 5.5 : Parameter settings for test case II... 60
Table 5.6 : The convergence and scaling properties for pressure pulse (Intel Xeon 5550, 2.67GHz). ... 65
Table 5.7 : Parameter settings for test case III ... 66
Table 5.8 : Displacement at point A for test case III ... 67
Table 5.9 : Parameter settings for test case IV... 69
Table 5.10 : The convergence and scaling properties for 3D Flag (Intel Xeon 5550, 2.67GHz)... 71
Table 5.11 : Parameter settings for test case V. ... 72
Table 6.1 : Material properties for cerebral aneurysm. ... 81
Table 6.2 : The convergence and scaling properties for the cerebral aneurism test case (Intel Xeon 5550). ... 82
LIST OF FIGURES
Page Figure 1.1 : Simple description of FSI. ... 1 Figure 1.2 : Mathematical FSI domain. ... 2 Figure 1.3 : Applications of fluid-structure interaction... 3 Figure 2.1 : General motion of a deformable body. ... 9 Figure 2.2 : Mesh descriptions. ... 13 Figure 2.3 : Deformation of continuum. ... 15 Figure 2.4 : Traction vectors in reference and current configurations. ... 18 Figure 2.5 : Description of an arbitrary control volume. ... 23 Figure 3.1 : Hexahedral element geometry. ... 28 Figure 3.2 : Three-dimensional unstructured mesh with a dual control volume. . 40 Figure 3.3 : Three-dimensional volume swept by the area A125. ... 43
Figure 5.1 : The geometric description of the first validation case. ... 54 Figure 5.2 : The computational coarse mesh M1 for the fluid and solid domains. 55 Figure 5.3 : The time variation of vertical displacement for an elastic bar
behind a rigid cylinder at Re= 200. ... 56 Figure 5.4 : The computed u−velocity vector component contours with the
streamlines for an elastic bar behind a rigid cylinder ... 57 Figure 5.5 : The partial view of the computational unstructured fine mesh M2
for both the fluid and solid domains. ... 61 Figure 5.6 : The radial and axial components of displacement on the inner
wall halfway along the pipe on meshes M1 and M2 ... 62 Figure 5.7 : The computed radial displacement contours for the propagating
pressure-pulse problem ... 63 Figure 5.8 : The computed solid radial displacement contours and the fluid
pressure contours ... 64 Figure 5.9 : The computed u−velociy profiles at t = 0.0069 ... 64 Figure 5.10 : The geometry of the FSI problem for a three-dimensional elastic
solid in a channel flow... 66 Figure 5.11 : The computational mesh for a three-dimensional elastic solid in a
rectangular channel flow. ... 67 Figure 5.12 : The computed streamtraces on the solid walls for an elastic bar
confined in a rectangular channel at Re= 40. ... 68 Figure 5.13 : The computed u−velocity vector component isosurfaces with the
streamtraces for an elastic bar confined in a rectangular channel at Re= 40. ... 69 Figure 5.14 : The geometry for test case IV. ... 70 Figure 5.15 : The computational mesh for test case IV... 72 Figure 5.16 : The computed q-criteria behind the flag... 73
Figure 5.17 : Vertical displacements at control points for test case IV. ... 73 Figure 5.18 : The deformed geometry with the vertical displacements contours
at several different time levels. ... 74 Figure 5.19 : The computational mesh for test case V... 75 Figure 5.20 : The time variation of the flexible ring location ... 75 Figure 5.21 : The computed x−velocity component contours with the deformed
ring at t= 2.4 s. ... 76 Figure 5.22 : The time variation of the total fluid area within the flexible ring ... 76 Figure 6.1 : The geometry of the cerebral aneurysm ... 80 Figure 6.2 : The computational all hexahedral conformal mesh for fluid and
solid domains ... 83 Figure 6.3 : The time variations of the mass flow rate at the inflow, outflow
sections and the inflow pressure ... 84 Figure 6.4 : The computed velocity vector magnitude contours at several
planes normal to y−axis... 85 Figure 6.5 : The wall shear stress lines and the vertical velocity component
contours... 86 Figure 6.6 : The time variation of the computed velocity vector magnitude
contours... 87 Figure 6.7 : The time variation of the computed stream-traces ... 88 Figure 6.8 : The time variation of the computed pressure iso-surfaces ... 89 Figure 6.9 : The time variation of the computed displacement vector
magnitude contours... 90 Figure 6.10 : The time variation of the computed wall shear stress contours... 91 Figure 6.11 : The time variation of total displacements at various points on the
A PARALLEL MONOLITHIC APPROACH FOR THE NUMERICAL SIMULATION OF FLUID-STRUCTURE INTERACTION PROBLEMS
SUMMARY
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 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−α 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−α 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 algo-rithm 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.
AKI ¸SKAN-YAPI ETK˙ILE ¸S˙IM˙I PROBLEMLER˙IN˙IN SAYISAL S˙IMÜLASYONU ˙IÇ˙IN
PARALEL MONOL˙IT˙IK B˙IR YÖNTEM ÖZET
Bu çalı¸smada akı¸skan-yapı etkile¸simi (FSI) problemlerinin paralel tam ba˘gla¸sık bir çözüm yakla¸sımıyla simülasyonuna yönelik özgün bir sayısal algoritma geli¸stir-ilmi¸stir. Problemin akı¸skan kısmi için daimi olmayan, sıkı¸stırılamaz Naviar-Stokes denklemleri, Arbitrary Lagrangian-Eulerian (ALE) formda kullanılmı¸stır. Akı¸skan için ALE tabanlı bu hareket denklemleri, yapısal olmayan bir sonlu hacimler yakla¸sımı ile ayrıkla¸stırılmı¸stır. Bu ayrıkla¸stırmada birincil de˘gi¸skenler kenar-merkezli bir sema ile konumlandırılmı¸stır. Bu konumlandırmada, hız vektörü komponentleri her hücre yüzünün orta noktasında tanımlanmı¸sken, basınç ise eleman merkezinde tanımlanmaktadır. Birincil de˘gi¸skenlerin bu ¸sekilde konumlanması kararlı bir sayısal sema olu¸smasını sa˘glamaktadır ve basınç-hız ba˘gla¸stırması için fazladan düzenlemeler yapılmasına gerek duyulmamaktadır. Bu sonlu hacimler yakla¸sımının en çekici tarafı, üniform Kartezyen çözüm a˘glarında, klasik MAC (Marker and Cell) semasında oldu˘gu gibi, Poison denklemi için klasik yedi-nokta Laplarca operatörüne yol açmasıdır. Bu durum, ço˘gunlukla akı¸skan alt problemi hesaplama kaynaklarının büyük kısmına ihtiyaç duydu˘gundan, özellikle büyük ölçekli akı¸skan-yapı etkile¸simi problemlerinin çözümünde çözüm veriminin artırılması acısından çok önemlidir. ALE tabanlı sonlu hacimler yakla¸sımında, do˘gruluk ve kararlılık acısından dikkate alınması gereken ba¸ska bir nokta da kullanılacak çözüm a˘gı hareketi yöntemidir. ALE yönteminde çözüm a˘gı akı¸skan ve kati sınırları arasındaki ara yüzü takip eder ve hareket denklemleri hareketli bir çözüm a˘gında ayrıkla¸stırılır. Bu durum da zaman integrasyon semasının do˘grulu˘gu ve kararlılı˘gı açısından, çözüm a˘gı hareketine özel ¸sartlar empoze edilmesini gerektirir. Bu ¸sartlar, ayrık geometrik korunun yasasının (DGCL) uygulanması ile sa˘glanmaktadır. Bu çalı¸smada, mesh hareketinden kaynaklı akılar geometrik korunum yasası ayrıkla¸stırma seviyesinde sa˘glanacak ¸sekilde hesaplanmaktadır. Süreklilik denkleminin her eleman içinde sa˘glanması için de özel bir dikkat sarf edilmi¸stir. Her elemen için süreklilik denkleminin toplamı çözüm alanı sınırlarına tam bir ¸sekilde indirgenebilir ki bu da özellik global kütle korunumu acısından çok önemlidir. Akı¸skan bölgesinde zaman integrasyonu için ikinci mertebeden geri farklar formülü kullanılmaktadır. Çözüm a˘gı deformasyonu algoritması için verimli cebirsel bir yöntem uygulanmı¸stır. Bu yöntemde, akı¸skan iç bölgesindeki çözüm noktaları, akı¸skan-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˘gi geometrik do˘grusal olmayan problemlere uygulanabilmektedir. Yapı bölgesi hareket denklemlerinin ayrıkla¸stırması Lagrangian bir çerçevede klasik Galerkin sonlu hacimler yöntemine dayanmaktadır. Üç boyutlu sonlu elemanlar ayrıkla¸stırmasında 8-nodlu, es-parametreli 6-yuzlu elemanlar kullanılmakta iken, iki boyutlu ayrıkla¸stırmalarda 4-nodlu dörtgen elemanlar kullanılmı¸stı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¸stirilmi¸s-α yöntemi uygulanmı¸stır. Bu yöntem, yapının hız ve yer de˘gi¸stirmeleri için Newmark tipi yakla¸sı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¸stirilmi¸s-α parametreleri seçilerek kolaca kontrol edilebilir. Yapısal dinamik denklemlerinin do˘grusalla¸stırılması, Newton tipi bir algoritmaya dayanmaktadır. Bu algoritmada, denklem sisteminin Jacobian matrisi her zaman adımında tam hesaplanmamaktadır. Bu yakla¸sım ALE tabanlı akı¸skan bölgesi denklemleri için de uygulanmı¸stır. Aksi taktirde tam Newton yöntemini, tam ba˘gla¸sık denklem sisteminde fazladan sıfır olmayan bloklar olu¸smasına neden olmakta ve bu da hafıza gereksinimlerini özellikle 3-boyutlu hesaplamalarda oldukça artırmaktadır. Akı¸skan ve yapı bölgelerine ait denklemlerin çözümü tam ba˘gla¸sık bir yakla¸sıma dayanmaktadır. Bu yakla¸sımda, akı¸skan ve yapı denklemleri tek bir denklem sistemi olu¸sacak ¸sekilde in¸sa edilmektedir ve bu denklem sistemi her zaman adımında tam ba˘gla¸sık ¸sekilde çözülmektedir. Akı¸skan ve yapı bölgeleri arasında ba˘gla¸stırma, akı¸skan-yapı ara yüzü boyunca birbirine uyumlu akı¸skan ve yapı çözüm a˘gları kullanılmasıyla basitle¸stirilmi¸stir. Hali hazırdaki monolitik FSI çözücü, ortaya çıkan tam ba˘gla¸sık denklem sisteminin çözümünde, ön-ko¸sulu bir Krylov alt uzay metodu kullanmaktadır. Hızın diverjansinin sıfır olması ko¸sulu nedeniyle olu¸san sıfır blok diyagonal, bütün sitem için verimli ön-ko¸sullandırıcıların uygulanmasını zorla¸stır-maktadır. Mevcut yöntemde ise, orijinal sistemdeki sıfır blok yerine ölçeklenmi¸s bir ayrık Laplacian olu¸sturan bir üst üçgen sa˘g ön-ko¸sullandırıcı uygulanmaktadır. Bu ön-ko¸sullandırma, matris-matris çarpımları nedeniyle sıfır olmayan eleman sayısında belirgin bir artı¸sa neden oldu˘gu için, ön-ko¸sullandırıcının sıfır olmayan blo˘gu hesaplama açısından daha az pahalı bir matris ise de˘gi¸stirilmektedir. Momentum denkleminde basınç gradyanlarina olan katkı, hız vektörlerinin ayrıkla¸stırıldı˘gı ortak eleman yüzeyini paylasan sa˘g ve soldaki elemanlardan oldu˘gundan, kullanılan bu matris sadece bu katkılardan kaynaklanan terimleri içermektedir. Bu yakla¸sım iteratif çözücünün yakınsama karakterini çok belirgin ¸sekilde etkilemese de, özellikle 3-boyutlu hesaplamalarda, hesaplama zamanı ve hafıza gereksinimlerinde ciddi azamlar sa˘glamaktadır. Mevcut tek seviye iteratif çözüm yakla¸sımı, sistemin simetrik olmayan do˘gası gere˘gi, kısıtlı aditif Schwarz ön-ko¸sullu 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ı¸stır. Algebrik denklemlerin do˘gasından kaynaklanan do˘grusal olmama durumu nedeniyle, yeterince tatmin edici bir yakınsama kriterine ula¸sıncaya kadar her zaman adımında alt-iterasyonlar uygulanmı¸stır. Mevcut ön-ko¸sullu Krylov alt uzay algoritması, matris-matris çarpımları ve kısıtlı aditif Schwarz ön-ko¸sullandırıcı uygulaması için PETSc (Portable, Extensible Toolkit for Scientific Computation) kütüphanesi kullanılmı¸stır. Bu kütüphane, do˘grusal ve do˘grusal olmayan denklem çözücülerinin paralel impilementasyonu için olu¸sturulmu¸s veri yapıları ve rutinleri içermektedir. Bütün akı¸skan-yapı ara yüzünün yapısal olmayan çözüm a˘gı, METIS kütüphanesi
kullanılarak alt parçalara ayrılmı¸stır. Bu kütüphane, yapısal olmayan grafik ve çözüm a˘glarının paralel programlamaya yönelik parçalanması için geli¸stirilmi¸s programlar içeren bir kütüphanedir.
Geli¸stirilen bu FSI çözücünün do˘grulu˘gunu 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ı¸stır. ˙Ilk do˘grulama örne˘gi oldukça popüler 2-boyutlu bir FSI test problemidir. Problem rijit dairesel bir engel arkasına yerle¸stirilmi¸s elastik bir çubuk ile etkile¸sime giren Newtonian bir akı¸stan olu¸smaktadır. Elastik çubuk akis bölgesinin alt ve üst duvarları arasında asimetrik yerle¸stirilmi¸stir ve bu akı¸skan yapı etkile¸simi senaryosunda, engelden kopan ve ilerleyen girdapların elastik çizim üzerinde indükledi˘gi periyodik osilasyonlar olu¸smaktadır. Bu test probleminde hem daimi hem de daimi olmayan durumlar dikkate alinmistir. Geli¸stirilen algoritmanın çözüm a˘gı yakınsama ve ölçeklenme karakterinin ortaya konulması amacıyla, üç fark çözüm a˘gı çözünürlü˘gü dikkate alınmı¸stır. Çözüm a˘gı çözünürlü˘günün, i¸slemci sayisinin, ILU(k) ön-ko¸sullandırıcı seviyesinin ve kısıtlı aditif Schwarz ön-ko¸sullandırıcıdaki overlap miktarının performansa etkileri ortaya konulmu¸stur. Bu test probleminde elde edilen sonuçlar literatürdeki birçok çözümle kar¸sıla¸stırılmı¸s ve mevcut çözücünün do˘grulu˘gu ispat edilmi¸stir. ˙Ikinci FSI test probleminde yine literatürde sıklıkla adres edilen 3-boyutlu bir konfigürasyon dikkate alınmı¸stır. Bu problem konfigürasyonu basitçe elastik arterlerden kan akisini simule etmektedir. Problemde sıkı¸stırılamaz viskoz bir akis esnek dairesel bir tüp ile çevrelenmi¸stir ve empoze edilen sinir ko¸sulları ile, elastik tüpte zamana ba˘glı radyal ve eksenel deplasmanlar olu¸sturacak ¸sekilde ilerleyen bir dalga çözümüne ula¸sılmaktadır. Bu 3-boyutlu test problemi icin de bir ölçeklenme testi yapılmı¸stır. Bu test için iki farklı çözüm a˘gı çözünürlü˘gü dikkate alin mistir. Radyal deplasmanlar için hesaplanan sonuçların literatür ile iyi uyum içinde oldu˘gu gösterilmi¸stir. Takip eden 3-boyutlu FSI test problemlerinden bir di˘geri dikdörtgensel bir kanal içine yerle¸stirilmi¸s elastik bir cisimden olu¸smakta ve daimi bir akis çözümü vermektedir. Cisim üzerinde bir kontrol noktasının deplasman çözümü verilmi¸s ve literatür ile olan uyumu gösterilmi¸stir. Dördüncü test problemi, dikdörtgensel bir engel arkasına yerle¸stirilmi¸s bir bayra˘gın girdap indüklü titre¸simlerini modelleyen 3-boyutlu bir konfigürasyondur. Bu problem için mevcut hesaplamalar oldukça yüksek çözünürlüklü bir çözüm a˘gında yapılmı¸stı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¸stir. Son test problemi, geli¸stirilen FSI algoritmasının, özellikle akı¸skan bölgesi yapı bölgesi ile tamamen çevrelendi˘gi durumlarda, kütle korunum kabiliyetinin ortaya konulması amacıyla tasarlanmı¸stır. Bu amaçla tasarlanan problem, paralel rijit duvarlar arasına simetrik olarak yerle¸stirilmi¸s iki boyutlu dairesel elastik bir yüzükten olu¸smaktadır. Bu dairesel yapı bir akı¸skan ile çevrelenmi¸s olmakla birlikte kendisi de bir akı¸skan bölgesini çevrelemektedir. Bu konfigürasyon basitçe kırmızı kan hücrelerinin dü¸sük Reynolds sayılı akı¸sta deformasyonunu simule etmektedir. Bu son test problemi, geli¸stirilen algoritmanın, kullanılan uyumlu FSI ara yüz ¸sartı ile makine hassasiyetinde kütle korunumunu sa˘gladı˘gını göstermektedir.
FSI çözücüsünün do˘grulanması amacıyla yapılan sayısal deneylerden sonra, mevcut algoritma kardiyovaskular akı¸skan-yapı etkile¸siminde sıklıkla kar¸sıla¸sın gerçekçi bir problemin çözümünde kullanılmı¸stır. Bu akı¸skan-yapı etkile¸simi problemi, damar çatallanma tepesinde anevrizma ihtiva eden bir beyin arterinde impulsif olarak ba¸slatılan bir FSI problemidir. Kan Newtonian bir akı¸skan olarak, damar
duvarı ise Saint Venant Kirchhoff modeline uygun bir elastik malzeme olarak modellenmi¸stir. Ba¸slangıçta kullanılan tamamen hegzahedral konformal çözüm a˘gı octree metodu kullanılarak olu¸sturulmu¸stur. Akı¸skan hız alanı, kan basıncı, duvar kayma gerilmeleri gibi birçok hemodinamik büyüklü˘gün yanında, zamana ba˘glı damar duvar deplasmanları hesaplanmı¸stır. ˙Iyi ölçeklenme karakteri bu problemde de elde edilmi¸stir.
Bu çalı¸smada geli¸stirilen FSI çözücüsünün geli¸stirme ve test a¸samalarında kullanılan yöntemler özetlenmi¸stir. Mevcut yöntemin avantaj ve dezavantajlarına de˘ginilmi¸s ve muhtemel gelecek uygulamalarından bahsedilmi¸stir.
1. INTRODUCTION
This chapter introduces the foundations of this study starting with a qualitative description of fluid-structure interaction (FSI) with underlying physical aspects in a simple manner, as well as its numerous occurrences in real life applications from an engineering point of view. Subsequently, various computational modeling approaches from the literature are discussed, which highlight the basis for the fundamental principles involved in FSI modeling. A brief description of the study is provided in the outline section at the end of this chapter.
1.1 Fluid Structure Interaction and Applications
The FSI problem is mainly characterized by the mutual interaction between a surrounding or internal fluid flow and a solid structure that is movable or deformable. The fluid flow causes deformation of the elastic structure, which, in turn, changes the boundary conditions of the fluid flow, and results in a dynamically coupled problem for most of the cases. As a result, fluid-structure interaction can be interpreted as a multi-field problem composed of computational fluid dynamics (CFD) and computational structural dynamics (CSD) as presented in Figure 1.1, which simply states that FSI is an overlap of these disciplines with a convenient coupling procedure.
Figure 1.2 : Mathematical FSI domain.
The most common engineering applications of fluid-structure interaction involve aeroelastic phenomena such as flutter, buffeting and divergence, design of offshore structures, structural effect of strong wind on bridges and tall buildings, flow induced vibrations in heat exchangers tubes, explosions or high-velocity impacts, large class of acoustic problems, membrane structures and many others. On the other hand, some of the most common applications in biomechanics include blood flow in the veins and arteries, deformations and aggregations of blood cells, dynamics of heart valves and animal locomotion. In Figure 1.3, sample applications with different physicals scales, which involve fluid-structure interaction are presented. An important example of FSI is a dynamic instability called wing flutter, which is not only evident for aircrafts at transonic speeds, but also for low subsonic aircrafts with slender wings of large aspect ratio as depicted by the first sub-figure, which shows a high performance sailplane. Flutter is a result of a positive feedback between force exerted on the wing by the fluid flow and the wing’s deflection leading to large amplitudes and stresses, where accurate prediction of fluid-structure interaction plays a crucial role in order to avoid catastrophic failure of the structure. Among various applications of fluid-structure interaction in biomechanics are the artery aneurysms. The second sub-figure of Figure 1.3 shows the case of cerebral aneurysm, which forms at a weak junction point of a brain artery that bulges out filling with blood. The growing aneurysm poses a risk of leak or rupture causing symptoms from severe headache to stroke-like symptoms, or even death. The initiation, growth and rupture of cerebral aneurysms are directly related to the blood pressure and wall shear stresses on the arterial wall, thus the prediction of these interaction between the blood flow and the aneurysmal wall plays an important role in estimating the risk factors for aneurysm. An other engineering application, that requires a through understanding of fluid-structure interaction is the
design and deployment of airbags (Figure 1.3(c)). Although, for the most of the cases, an airbag is beneficial if not a lifesaver during an accident, there is still a risk of injury for the passengers located within the deployment region of the airbag. Fast deployment and longer period of inflation can be considered to be design criteria for a reliable airbag, and achieving this criteria requires a detailed simulation of the involving fluid-structure interaction.
(a)
(b) (c)
Figure 1.3 : Applications of fluid-structure interaction
1.2 Literature Review
The fluid-structure interaction problem is probably one of the most addressed multi-physics problems in the literature, not only for its numerous applications in engineering and biomechanics, but also for the computational challenges involved in its numerical modeling and simulation. One of the most well-known methods used to model the interaction between structure and fluid is the Arbitrary Lagrangian Eulerian (ALE) method as described in Hirt el al. [1]. 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. This differs from the standard Eulerian formulation in a way that the mesh movement has to fulfill special conditions in order to maintain the
accuracy and the stability of the time integration scheme. This condition is satisfied by the enforcement of the so-called discrete geometric conservation law (DGCL) as coined by Thomas and Lombard [2]. Although the GCL is satisfied easily in the continuous sense, its discrete implementation may not be trivially satisfied. The ALE time integration scheme developed by Koobus and Farhat [3] is based on more continuous time integration of the fluxes and offers second-order accuracy in time obeying the GCL. Geuzaine et al. [4] have showed that the GCL is neither a necessary nor a sufficient condition for an ALE scheme to preserve its order of time accuracy established on fixed meshes. The ALE approach was subsequently adopted within the finite element context to solve free surface problems of incompressible viscous fluid flow [5]. In the case of an FSI problem, the deformable fluid-structure interface is taken into account and the fluid points at the fluid-solid interface are moved in a Lagrangian way [6]. Although the Eulerian approaches have also been used for the simulation of the fluid-structure interaction problems around complex geometrical shapes [7–9], the use of these approaches for large mesh deformations requires relatively high mesh resolution in large percentages of the fluid computational domain, or adaptive mesh refinement, in order to properly capture the viscous effects.
The goal of the present study is to incorporate an unstructured finite volume algorithm based on an Arbitrary Lagrangian Eulerian formulation with the classical Galerkin finite element discretization for the solid domain in order to develop a novel algorithm for the large-scale simulation of fluid-structure interaction problems in a fully coupled form, where the fluid is modeled by the incompressible Navier-Stokes equations and the structure is modeled by the nonlinear Saint Venant-Kirchhoff model. The Arbitrary Lagrangian-Eulerian formulation used for the governing incompressible Navier-Stokes equations of the fluid domain is based on the side-centered unstructured finite volume method where the velocity vector components are defined at the mid-point of each cell face while the pressure is defined at the element centroid. The present arrangement of the primitive variables leads to a stable numerical scheme and it does not require any ad-hoc modifications in order to avoid odd-even pressure decoupling or spurious pressure modes on unstructured meshes [10]. This approach was initially used by Maliska and Raithby [11] to calculate the contravariant velocity components that enter the mass conservation constraint on boundary-fitted meshes and later by Hwang [10]
and Rida et al. [12] for the solution of the incompressible Navier-Stokes equations on unstructured triangular meshes. More recently, Sahin has incorporated the present approach with a parallel non-nested geometric multigrid method in order to simulate the three-dimensional viscoelastic wake instabilities [13]. 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 [14] which is very important for the efficient solution of the large-scale FSI problems in 3D, since the numerical simulation of the FSI problems in general requires large computational resources and it is typically the fluid subproblem that requires the most computational resources in the coupled system [15]. To our best knowledge, the present arrangement of the primitive variables is not considered for the FSI problems. In the present work, a special attention will be given to satisfy the continuity equation exactly within each element and the summation of the continuity equations can be exactly reduced to the domain boundary, which is important for the global mass conservation. The mesh deformation within the fluid domain is achieved by using an algebraic approach based on the minimum distance function at each time level while avoiding re-meshing in order to enhance numerical robustness. The time integration method for the structure domain is based on the Newmark type generalized−α method while the second-order backward difference formula (BDF2) is used in the fluid domain.
There are mainly two solution procedures for the numerical simulation of fluid-structure interaction problems: partitioned (segregated) [16–18] or fully coupled (monolithic) [19–24] methods. In the partitioned approach, separate solvers are utilized for the fluid and structure subproblems. The main advantage of the partitioned approach is the ability to reuse existing solvers which allows the application of different, possibly more efficient, computational methods specifically developed for either the fluid or the structure subproblem. Both explicit or implicit methods can be used in order to couple the fluid and structure solvers in partitioned approaches. In explicit partitioned methods, which are also known as loosely or weakly coupled methods, typically a fixed point (Picard) iteration is employed to obtain a coupled solution. Although the implementation of this approach is relatively easy, it does, however, suffer some serious drawbacks. The fixed point iterations tend to converge
slowly and the iterations may diverge in the presence of strong fluid-structure coupling due to the high fluid/structure density ratio which causes to the so-called artificial added mass effect [25]. For strong coupling in partitioned procedures, on the contrary, several fluid and structure computations are performed in a single time-step until a satisfactory convergence tolerance is reached. This approach, however, requires costly sub-iterations and the sub-iteration convergence may not be guaranteed. In addition, these standard partitioned methods can not satisfy the incompressibility constraint of the fluid [26] during standard alternating FSI iterations when the fluid domain is entirely enclosed by all Dirichlet boundary conditions. In a fully coupled (monolithic) approach, the fluid and structure equations are discretized and solved simultaneously as a single equation system for the entire problem. However, this requires an efficient numerical technique for the solution of a large system of coupled nonlinear algebraic equations, which poses the major challenge of monolithic FSI approaches, especially in large scale problems. Although monolithic solvers are believed to be too expensive for use in large-scale problems, more recent studies [27, 28] demonstrate that monolithic solvers are competitive even in the case of weak fluid-structure interactions problems. In this context, Muddle et al. [24] presented a block preconditioner for the efficient solution of the linear systems by Krylov subspace solvers. Barker and Cai [22] developed a two-level Newton-Krylov-Schwarz algorithm for monolithic coupling in fluid-structure interaction, with special consideration of the application area of simulating blood flow. Gee et al. [23] applied an algebraic multigrid technique to the entire fluid-structure interaction system of equations. Tezduyar et al. [29] described a set of solution techniques for computation of fluid-structure interaction problems and the solution of the fully discretized fluid and structural mechanics equations is done using the block-iterative, quasi-direct, and direct coupling techniques. Badia et al. [20] proposed semi-implicit algorithms based on inexact block-LU factorization of the linear system obtained after the space-time discretization and linearization of the FSI problem. In the present study, the original system of equations is preconditioned with an upper triangular right preconditioner which results in a scaled discrete Laplacian instead of a zero block in the original system due to the divergence-free constraint. Then a one-level restricted additive Schwarz preconditioner with a block-incomplete factorization within each partitioned sub-domains is utilized for the resulting fully coupled system. The implementation of the preconditioned
Krylov subspace algorithm, matrix-matrix multiplication and the restricted additive Schwarz preconditioner are carried out using the PETSc [30] software package developed at the Argonne National Laboratories.
1.3 Outline
In this thesis, in order to keep the integrity of the overall discussion and provide the reader with the basic concepts, some fundamental aspects of continuum mechanics directly related to subsequent analysis are introduced before going into the detailed numerical formulations for the governing equations of the fluid and the structure models. Chapter 2 is devoted to these fundamental concepts of kinematics of deformable continua; Lagrangian, Eulerian and ALE (Arbitrary Lagrangian Eulerian) descriptions; stress and strain measures and constitutive relations in solids and fluids, and finally the conservation laws and governing equations of motion.
In Chapter 3, the numerical methods used the discretize the governing equations for solid and fluid domains are addressed with detailed discussion on the Galerkin finite element based discretization and time integration method for the structure solver, and the unstructured finite volume method used to discretize the ALE based conservation equations for the fluid solver. The concepts of geometric conservation law (GCL) for accurate time integration in moving domain problems are presented with application of GCL to the current finite volume discretization of fluid equations. The mesh deformation method used in the present numerical approach is also introduced.
The fully coupled fluid-structure interaction algorithm used to couple the structure and the fluid solvers is presented in Chapter 4. The fluid-structure interface conditions, the discretization of the interface equations and the assembly of the monolithic system of equations are discussed along with the proposed preconditioned Krylov subspace method.
Chapter 5 is dedicated to the application of the present fully coupled fluid-structure interaction solver with solutions to several FSI benchmark problems addressed in the literature. The first problem is a two dimensional benchmark problem where an elastic bar is placed behind a rigid cylinder and shedding vortices induce oscillations of the elastic bar. Both steady and unsteady FSI is considered in the benchmark.
The second benchmark is three dimensional configuration modeling a pressure pulse propagating in an elastic channel, which is a simplified model of blood flow through elastics arteries. The third reference problem is FSI of a three dimensional elastic solid confined in a channel flow, which results in a steady fluid structure interaction. The fourth benchmark case, which involves three dimensional vortex induced vibrations of an elastic flag, is a rather difficult one requiring a considerably high mesh resolution as a result of a relatively high Reynolds number regime. 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 model is created involving a flexible circular ring placed symmetrically between two parallel rigid plates. This configuration mimics the deformation of a red blood cell inside a capillary wall at a very low Reynolds number. In Chapter 5, the mesh convergence and the scaling properties of the present algorithm is also presented by testing the present monolithic solver on several different mesh resolutions for the benchmark problems, and for different number of processor setting as well. The results are presented in a detailed manner. Chapter 6 is devoted 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. 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.
Final remarks on methods used to develop and test the present FSI solver are summarized in Chapter 7. Advantages and the drawbacks of the present solver are addressed with possible future applications.
2. MATHEMATICAL MODELLING
2.1 Kinematics of a Deformable Body
The general motion of a deformable body is illustrated in Figure 2.1. The initial state of the domain of the deformable body at t = 0 is called the initial configuration and denoted by Ω0. Unless specifically stated, throughout the ongoing discussions,
the initial configuration is taken as the undeformed and the reference configuration, to which the various equations are referred in order to describe the motion and the deformation of the body. The domain occupied by the deformed body at time t is denoted byΩ, which is referred as the current configuration. The position vector X of a material point in the reference configuration with respect to the Cartesian basis EI is
given by
Figure 2.1 : General motion of a deformable body. Deformed (Ω) and undeformed Ω0configurations.
X= XiEi (2.1)
The coordinate X that labels a material point in the reference configuration at t = 0 is called the material coordinate or the Lagrangian coordinate for that point, and for a given material point the material coordinate is time invariant. In general, the positions of the particles in the current configuration are designated by the coordinates x with respect to an alternative Cartesian basis ei. The bases EI and ei are taken to
be coincident in the remainder of this study. The coordinates x are called the spatial or the Eulerian coordinates and give the spatial potion of the particle in the deformed configuration,thus the position vector in the current configuration is given by
x= xiei (2.2)
The motion of the deformable body is described by a mapping functionφ(X,t), which maps the reference configuration into the current configuration. Since coincident Cartesian bases are used for both configurations, the spatial coordinates coincide with the material coordinates for t= 0.
x = φ(X,t) (2.3)
X = x(X, 0) ≡φ(X, 0) (2.4)
The displacement of a material point is defined as the the difference between its current position and its initial(reference) position, and given as a function of material coordinates as
u(X,t) =φ(X,t) −φ(X, 0) =φ(X,t) − X (2.5) Describing the motion of the body as in (2.3), the velocity of a particle can be defined as the time rate of change of its position vector with holding the material coordinate X constant. It should be noted that the velocity of a particle is a spatial vector, while it is expressed in terms of the material coordinates of the particle in (2.6). Thus, a more consistent expression can be defined in terms of the spatial coordinate by inverting (2.3), which is given by (2.7).
v(X,t) = ∂φ(X,t)
∂t =
∂u(X,t)
∂t (2.6)
v(X,t) = v(φ−1(x,t),t) (2.7)
The acceleration of the particle is the time rate of change of its velocity, and can be expressed as a function of material coordinates as
a(X,t) = ∂v(X,t)
∂t =
∂2u(X,t)
∂t2 (2.8)
In Equations (2.6) and (2.8), the time derivatives are calculated keeping the material coordinate X constant. In other words, these expressions measure the time rate of change of the given quantity for a specific point initially at X. A time derivative of this type is called the material time derivative, and frequently denoted by the symbol D()/Dt. When a spatial variable is defined by the spatial coordinate x, as in the case of an Eulerian description of velocity as v(x,t), the material time derivative can be obtained using the chain rule of differentiation as follows
D v D t = ∂v(x,t) ∂t + ∂v(x,t) ∂x ∂φ(X,t) ∂t (2.9)
In Equation (2.9), the fist term is called the spatial time derivative, and calculated at a fixed spatial coordinate, while the second term involving the particle velocity is called the convective derivative. In tensor notation Equation (2.9) can be written as
D v D t =
∂v
∂t + v ·∇v (2.10)
2.2 Lagrangian, Eulerian and ALE Descriptions
In the previous section, two main approaches for describing the motion of a continuum are discussed based on expressing a dependent variable of the domain by either the material coordinates of particles or the spatial coordinates of the domain, where the former approach is called the Lagrangian description, and the latter is called the Eulerian description of the domain. The behavior of the computational grid is directly related to the choice of independent variables for describing the domain.
Since the Lagrangian viewpoint is based on following the material particle during its motion, the computational mesh is selected so that it is permanently attached to the continuum. So, each node of the Lagrangian mesh labels an individual material particle and follows the particle during the deformation of the domain. On the other hand, in an Eulerian description the domain variables are based on the spatial coordinates, and the physical quantities of particles are observed at a fixed region of the domain when the particles are passing through it at an instant of time. Therefore, in an Eulerian description the domain is described by a fixed computational grid with time invariant spatial coordinates. These different behaviours of Lagrangian and Eulerian meshes are simply illustrated in Figure 2.2 for a one-dimensional example. Figure 2.2 shows the space-time description of elements starting from the reference configuration at t = 0 and the location of the element nodes and material particles at an instant of time t. For the Lagrangian description, the computational grid deforms with the material domain and the grid nodes and particles remain coincident during deformation. Whereas for the Eulerian description, the computational grid nodes are coincident with the spatial coordinates which remain fixed, and the material points move relative to the Eulerian mesh.
In most of the cases, the Lagrangian mesh description is more suited to problems encountered in structural mechanics. Since the material particles cannot pass through an element boundary, each finite element in a Lagrangian mesh always contains the same material particles, which is computationally advantageous when treating materials with history dependent constitutive relations, and also simplifies interface tracking between different materials. On the contrary, the problems encountered in fluid mechanics are most widely described by Eulerian algorithms because of the fact that the behavior of the fluid is independent of its history, and there is often no need to express the motion with respect to a reference configuration. Additionally, Eulerian description is more favorable when dealing with complex material motion with excessively large distortions as in the case of fluid flow.
In terms of computational effectiveness, both the Lagrangian and Eulerian descriptions have their own advantages and disadvantages when selecting an appropriate mesh description for a specific problem. For the Lagrangian description, for example, since the mesh deforms with the deforming continuum, the boundary nodes always
Figure 2.2 : Mesh descriptions.
remain on the physical boundary, which permits an easy treatment of boundary conditions. However, if the continuum undergoes excessively large deformations, the Lagrangian mesh can get severely distorted, which causes a decrease in element accuracy and limits the magnitude of deformation that can be resolved by a Lagrangian mesh. In contrast to the Lagrangian description, the Eulerian mesh description can easily deal with material motion with large distortions without any change in element accuracy, but since the spatial node of the Eulerian mesh is fixed, the boundary nodes do not remain on moving boundaries, which poses a significant difficulty in imposing boundary conditions, and requires special tracking methods for treating moving boundaries.
In order to benefit from the domain descriptions of the Lagrangian and the Eulerian approaches, or when it is almost impossible or not practical to define the whole computational domain as either pure Lagrangian or pure Eulerian as in the case of complex fluid flows interacting with moving and/or deforming material boundaries, a combined computing method which utilizes both Lagrangian and Eulerian aspects of
domain description is used, which is called the Arbitrary Lagrangian Eulerian (ALE) method. In an ALE computational domain, the domain vertices can move with the fluid in a Lagrangian way, remain fixed in an Eulerian way, or move arbitrarily prescribed to establish a continuous rezoning at each transient level, which is simply illustrated in the one-dimensional example of 2.2. In a fluid-structure interaction problem, for example, when an ALE based approached is used for domain description, the vertices on the fluid structure interface move with the interface in a Lagrangian way, while the surrounding vertices move arbitrarily using a convenient domain deformation algorithm to obtain a sufficient level of grid quality for accurate representation of the domain. The details of the present ALE based computing method and grid deformation strategies are discussed in detail in section 3.2.
2.3 Structure
2.3.1 Deformation gradient
A key quantity of significant importance in finite deformation theory is the deformation gradient F, which relates the relative spatial position of two neighboring particles after deformation to the relative material position of particles in the reference configuration (see Figure 2.3), and defined as
F= ∂φ
∂X ≡
∂x
∂X (2.11)
An infinitesimal line segment dx in the current configuration can be related to the corresponding line segment dX by the following relation
dx= FdX (2.12)
The deformation gradient can be explicitly given in the following matrix form in a Cartesian coordinate system.
Figure 2.3 : Deformation of continuum. F= ∂x1 ∂X1 ∂x1 ∂X2 ∂x1 ∂X3 ∂x2 ∂X1 ∂x2 ∂X2 ∂x2 ∂X3 ∂x3 ∂X1 ∂x3 ∂X2 ∂x3 ∂X3 (2.13)
2.3.2 Volume change and compresibility
An infinitesimal volume element dV in the current configuration can be related to the corresponding elemental material volume dV0 in the reference configuration by the
following push-forward operations of the undeformed edge vectors
dx1= FdX1= ∂φ ∂X1 dX1 (2.14) dx2= FdX2= ∂∂φ X2 dX2 (2.15) dx3= FdX3= ∂∂φ X3 dX3 (2.16)
The deformed elemental volume in the current configuration can be obtained as the triple product of these elemental vectors as
dV = dx1· (dx2× dx1) = ∂∂φ X1 · ∂φ ∂X2 × ∂φ ∂X3 dX1dX2dX3 (2.17)
In Equation (2.17), the triple product gives the determinant of the deformation gradient, which is also called the Jacobian J , and the volume change can be given in terms of the Jacobian as in Equation (2.19).
J= detF = ∂φ ∂X1 · ∂φ ∂X2 × ∂φ ∂X3 (2.18) dV = JdV0 (2.19)
A powerful property of the Jacobian determinant is that it can be used to relate the integrals in the current configuration to the ones in the reference configuration as
Z
Ω f dV =
Z
Ω0 f JdV0 (2.20)
Since the volume cannot have a negative value, physically the Jacobian determinant should be positive, i.e. J> 0. Additionally, if for a material the volume is preserved for all values of the deformation gradient F, then the material is defined as incompressible, and the value of the Jacobian is given by Equation (2.21).
J= detF = 1 (2.21)
2.3.3 Strain
A general measure for strain can be established by considering the deformation of an infinitesimal line segment dX in the reference configuration to the corresponding line segment dx in the current (deformed) configuration. The lengths of these infinitesimal line segments can be determined by the following dot products of line vectors.
(dS)2= dX · dX (2.22)
(ds)2= dx · dx = dX · FTFdX
In Equation (2.23), the length of the line segment in the current configuration is expressed in terms of the material vectors, where C is the right Cauchy-Green deformation tensor, which is given in terms of the deformation gradient F as
C= FTF (2.24)
As the body deforms from the reference to the current configuration, the length change of the material vector dX can be expressed with respect to the undeformed length as in Equation (2.25), where the material tensor E is called the Green-Lagrange strain tensor, and defined by the relation given in Equation (2.26).
(ds)2− (dS)2= 2 dX · EdX (2.25)
E=1 2(F
TF− I) = 1
2(C − I) (2.26)
The Green-Lagrange strain tensor E is probably the most commonly used strain measure in finite deformatin analysis, and it can be easlly seen from Equation (2.26) that E is symmetric, i.e. E= ET.
2.3.4 Stress measures
In linear small displacement analysis, since there is no distinction made between the deformed and the undeformed configurations, i.e. they are assumed to be geometrically coincident, the state of stress can be defined by a single stress measure, which leads to the well-known Cauchy stress tensor. In finite deformation analysis, on the other hand, alternative stress measures which are based on the reference configuration of the structure can also be defined. In this section, these stress measures are introduced, which will be later used to express the governing equations of the structure with respect to the reference configuration.
Figure 2.4 shows the traction vectors t and T that are acting on area elements ds and dS with normal vectors n and N, which correspond to the current and reference configurations of the structure respectively. The elemental force df in the current
Figure 2.4 : Traction vectors in reference and current configurations.
configuration is defined by Equation (2.27) by either the true traction t or the nominal traction T.
df= tds = TdS (2.27)
t= t(x,t, n), T= T(X,t, N) (2.28)
The traction t, which is the force measured per unit surface are defined in the current configuration, is related to the surface normal n by Cauchy’s theorem given in Equation 2.29, where σσσ is called the Cauchy stress tensor. The tensorσσσ is a true measure of stress in the current configuration.
t=σσσn (2.29)
Using Nanson’s formula (2.30), an other stress measure can be defined by the transformation from the current configuration to the reference configuration. Here,ΠΠΠ is called the first Piola-Kirchhoff stress tensor, which is evident from Equation (2.31) thatΠΠΠis not symmetric, and it is a two field tensor with one basis referred to the current and the other to the reference configuration.
