Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
The projection immersed boundary method for compressible flow and its application in fluid-structure interaction simulations of parachute systems
(USC Thesis Other)
The projection immersed boundary method for compressible flow and its application in fluid-structure interaction simulations of parachute systems
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
The Projection Immersed Boundary Method for Compressible Flow and Its Application in Fluid-Structure Interaction Simulations of Parachute Systems by Hang Yu A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (MECHANICAL ENGINEERING) December 2024 Copyright 2024 Hang Yu Table of Contents List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x Chapter 1: Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Motivation and background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Objectives and thesis outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Related publications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Chapter 2: The projection immersed boundary method for compressible viscous flow . . . . . 5 2.1 Immersed boundary method: a literature survey . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.3 Governing equations of compressible projection IB methods . . . . . . . . . . . . . . . . . 9 2.3.1 Singular source terms and jump conditions . . . . . . . . . . . . . . . . . . . . . . 9 2.3.2 Differential-algebraic system of governing equations . . . . . . . . . . . . . . . . . 12 2.4 Numerical discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.4.1 Spatial discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.4.2 Half-explicit temporal discretization . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.5 Preliminary studies on permeable surface . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.6 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.6.1 Convergence tests of moving ball . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.6.2 Flow past circular cylinder . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.6.3 Acoustic scattering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 2.6.4 Shock reflection from a wedge . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 2.7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 Chapter 3: Tikhonov regularization of the projection immersed boundary method . . . . . . . 42 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 3.2 The projection immersed boundary method . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.3 Challenges of the projection method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.4 A regularized projection immersed boundary method . . . . . . . . . . . . . . . . . . . . . 49 3.4.1 Weighted least square formulation of the projection method . . . . . . . . . . . . . 50 3.4.2 The Tikhonov regularized formulation . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.4.3 Constructing the smoothing seminorm with the graph Laplacian . . . . . . . . . . . 52 ii 3.4.4 The regularized immersed boundary equation . . . . . . . . . . . . . . . . . . . . . 54 3.5 Regularization of compressible IB method . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.6 Preliminary investigation of the oscillatory patterns near boundary . . . . . . . . . . . . . . 59 3.7 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 3.7.1 Flow passing cylinder with regularization . . . . . . . . . . . . . . . . . . . . . . . 63 3.7.2 Impulsively started rotating cylinder . . . . . . . . . . . . . . . . . . . . . . . . . . 68 3.7.3 Fluid-structure interaction simulation of a fixed cylinder and a flexible beam . . . . 71 3.8 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 Chapter 4: Computational framework of fluid-structure interaction simulations . . . . . . . . 77 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 4.2 Thin-shell structure dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 4.3 Immersed boundary method with patch-based structured adaptive mesh refinement (IB-SAMR) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 4.3.1 IB-SAMR grid structures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 4.3.2 IB-SAMR grid hierarchy generating and regridding . . . . . . . . . . . . . . . . . . 82 4.3.3 Coarsening and refining operators . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 4.3.4 Operators between IB Lagrangian grid and SAMR Eulerian grids . . . . . . . . . . 86 4.3.5 IB-SAMR integration cycle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 4.4 Partitioned loosely coupled FSI time integrator . . . . . . . . . . . . . . . . . . . . . . . . 88 4.5 FSI in non-inertial frames of references . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 4.6 Verifications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 4.6.1 Airdrop of sphere . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 4.6.2 Drop of Parachute Test Vehicle (PTV) . . . . . . . . . . . . . . . . . . . . . . . . . 97 4.7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 Chapter 5: Simulation of the drogue parachute cluster of the Orion Capsule Parachute Assembly System (CPAS) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 5.2 Problem setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 5.2.1 Parachute system descriptions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 5.2.2 Flight conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 5.2.3 Mesh information . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 5.2.4 Simulation schedules . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 5.3 Preliminary results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 5.3.1 Reefed stage 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 5.3.2 Reefed stage 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 5.3.3 Fully open stage 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 5.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 Chapter 6: Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 6.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 6.2 Future outlooks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 A Derivative of discontinuous function and the surface singularity . . . . . . . . . . . . . . . . 130 iii List of Tables 2.1 Coefficients of the 3-stage SSPRK [1] formulated as a DAE integrator [2]. . . . . . . . . . . 22 2.2 Drag coefficients for Ma = 0.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 2.3 Comparison of drag coefficients for Re = 20 and 40. . . . . . . . . . . . . . . . . . . . . . 35 2.4 Comparison of drag coefficients for Re = 100 and 200. . . . . . . . . . . . . . . . . . . . . 37 3.1 Drag coefficients, separation angle and wake characteristic lengths for Re = 20 and 40, with comparison to previous literature. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 3.2 Drag/lift coefficients and Strouhal number of vortex shedding for Re = 200, with comparison to previous literature. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 3.3 The amplitude, drag coefficients, and the Strouhal number of the benchmark FSI problem, with a comparison to previous literature . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 4.1 The expected and actual Reynolds number, terminal velocity, and drag coefficients at terminal states. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 5.1 Materials properties of ribbons. The density and Poisson ratio for all ribbons are 1.15 g/cm3 and 0.33. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 5.2 Materials properties of tensile cords. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 5.3 Initial flight conditions of simulations of stages 1, 2, and 3 . . . . . . . . . . . . . . . . . . 105 5.4 Summary of the simulation schedules . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 iv List of Figures 1.1 CPAS systems and deployment sequence [3] . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.1 Two-dimensional schematic layouts of the finite-volume Eulerian grid of the fluid computational domain Ω and the Lagrangian grid of the immersed boundary Γ. . . . . . . . 17 2.2 Velocity field of moving ball problem at a t = 0.5 with no-slip and imposed temperature condition; internal induced flow inside the ball shown. . . . . . . . . . . . . . . . . . . . . 31 2.3 L1 error norm of the moving ball problem non-dimensionalized by the reference fields as a function of ∆t: (a) no-slip and imposed temperature condition and (b) stress and heat flux condition; second-order reference slop is plotted. . . . . . . . . . . . . . . . . . . . . . . . 32 2.4 L1 error norm of the moving ball problem non-dimensionalized by the reference fields as a function of ∆x: (a) no-slip and imposed temperature condition and (b) stress and heat flux condition; second-order reference slop is plotted. . . . . . . . . . . . . . . . . . . . . . . . 32 2.5 L1 error norm non-dimensionalized by the reference fields of the moving ball problem with time-dependent boundary temperature: (a) temporal and (b) spatial error, respectively. . . . . 33 2.6 Streamlines near the cylinder for (a): Re = 20 and (b) : Re = 40. Nb = 161. . . . . . . . . . 35 2.7 Vorticity contour for (top) Re = 100 and (bottom) Re = 200. Nb = 161. . . . . . . . . . . . 36 2.8 Evolution of drag coefficients with time (top) after vortex shedding developed for (left) Re = 100 and (right) Re = 200, and the evolution of lift coefficients with time (bottom) for the corresponding Reynolds number. Solid: Nb = 322; Dash: Nb = 161; Dash-dot: Nb = 81. x-axis is non-dimensionalized and shifted for better comparison. . . . . . . . . . . 37 2.9 Contours of instantaneous non-dimensionalized pressure fluctuations for acoustic scattering probem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 2.10 Pressure fluctuations as a function of time at x = 2 and y = 0: simulations (solid) and analytical solution (dash). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 2.11 Shock reflection from a wedge: (a) pressure and (b) Mach number; domain under the wedge is masked for clarity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.1 Singular values of EBL against the total variation of the right singular vectors. (a): ∆s = 0.0698; (b): ∆s = 0.0980; (c): ∆s = 0.1253;(d): ∆s = 0.1564. . . . . . . . . . . . . . . 48 3.2 Condition number of EBL for various ∆s/∆x. The region where EBL is singular is gray-shaded. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 v 3.3 The boundary force on the cylinder with different ratios of Lagrangian-Eulerian grid spacing. (a): ∆s/∆x ≈ 1.6; (b): ∆s/∆x ≈ 1.3; (c): ∆s/∆x ≈ 1. . . . . . . . . . . . . . . . . . 49 3.4 The eigenvectors of the graph Laplacian for the cylinder problem described in Section 3.3 with ∆s/∆x ≈ 1.3, and its corresponding eigenvalues, non-dimensionalized by ∆x∆y (a) and (b) : the 10th and 40th eigenvectors of the graph Laplacian. . . . . . . . . . . . . . . . . . . 53 3.5 Total variation of eigenvectors and their corresponding eigenvalues for the Graph Laplacian operator constructed for the cylinder grid. . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.6 The oscillatory patterns at boundaries when solving the linear equation for the boundary force, demonstrated with a flat immersed boundary with 0 ≤ ξ ≤ 1. . . . . . . . . . . . . . 61 3.7 The solution of E gL g ˆf = r with operators E g and L g constructed when markers at boundaries are categorized as ghost markers. The right hand side r is the same as in Figure 3.6b. . . . . 62 3.8 The residual norm ∥Eu∥ 2 2 vs. smoothing seminorm ˆf 2 R at different regularization parameters λ for the Tikhonov-type regularization, Re = 40. (a): ∆s/∆x = 1.3; (b): ∆s/∆x = 1.0; (c): ∆s/∆x = 0.7. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 3.9 The residual norm ∥Eu∥ 2 2 vs. smoothing seminorm ˆf 2 R at different regularization parameters λ for the Lavrentiev-type regularization, Re = 40. (a): ∆s/∆x = 1.3; (b): ∆s/∆x = 1.0; (c): ∆s/∆x = 0.7. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 3.10 The x-direction boundary force on the cylinder for ∆s/∆x = 1.3. (a): without regularization; (b): Tikhonov-type regularization with λ = 0.005; (c): Lavrentiev-type regularization with λ = 0.5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.11 The x-direction boundary force on the cylinder for ∆s/∆x = 1.0. (a): without regularization; (b): Tikhonov-type regularization with λ = 0.005; (c): Lavrentiev-type regularization with λ = 0.5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.12 Streamlines of the flow passing a cylinder at Re = 40. The definitions of characteristic lengths a, b, and l of the recirculation wake are annotated. The end of the recirculation zone and the centers of the vortices can trivially be probed from the simulation results because both x and y velocities change sign there. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 3.13 The residual norm ∥Eu∥ 2 2 vs. smoothing seminorm ˆf 2 R at different λ(a): Tikhonov-type regularization; (b): Lavrentiev-type regularization. . . . . . . . . . . . . . . . . . . . . . . 69 3.14 The calculated tangential boundary force at different regularization parameters λ. (a): Tikhonov-type regularization; (b): Lavrentiev-type regularization. . . . . . . . . . . . . . . 70 3.15 Error of the flow field u and the boundary force f compared to the exact solution versus the grid spacing ∆x. (a)-(b): Tikhonov-type regularization; (c)-(d): Lavrentiev-type regularization. 71 3.16 Schematic diagram of the FSI benchmark problem. E1, E2, E3, and E4 are four edges of the beam Ωb. E5 is the boundary arc of the cylinder. . . . . . . . . . . . . . . . . . . . . . . . . 72 3.17 Flow vorticity fields and structural positions at the times when the displacements at beam tip reach maximum. Tikhonov-type regularization with λ = 0.001 is used. . . . . . . . . . . 74 3.18 Time histories of the beam tip displacement for the regularized projection formulation. . . . 75 vi 4.1 Schematic layout of three-level patch-based AMR grids [4]. (a) the composite of mesh; (b) the AMR mesh hierarchy as sequence of levels. Illustrations from SAMRAI documents by LLNL. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 4.2 Re-meshing process in patch-based AMR framework [5]. (a): the cells are tagged for refinement by users; (b): the tagged cells and some untagged cells are grouped into rectangular boxes; (c): rectangular boxes partitioned into patches and dispatched among MPI ranks. Illustrations from SAMRAI documents by LLNL. . . . . . . . . . . . . . . . . 83 4.3 Local two-dimensional AMR hierarchy with refinement ratio 2 for demonstrating the linear refining and conservative coarsening. The finer cells in level l +1 corresponding to (i, j) in level l are (2i,2 j),(2i+1,2 j),(2i,2 j +1),(2i+1,2 j +1). . . . . . . . . . . . . . . . . . . 85 4.4 Loosely coupled partitioned time stepping from n to n+1 for structure data Xs and fluid data U . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 4.5 The inertial, non-inertial (computational), and body-fixed frame of the rigid body. xf and xcm are the origins of the non-inertial and body-fixed frame in the inertial frame, respectively 92 4.6 Partitioned time stepping from n to n + 1 for structure data Xs , fluid data U and the data of rigid body Xb. Data exchange process I and III: passing position/velocity of the body/structure to the fluid as the input of IB solvers; II and V: passing boundary force from IB solvers to the structure and body; IV: passing position/velocity from body to structure; V I: passing force from the structure to the body. . . . . . . . . . . . . . . . . . . . . . . . 93 4.7 The properties of sphere airdrop at Re = 200. (a): the sphere dropping velocity; (b): the drag force on the sphere; (c): fluid velocity near the ball at terminal state. . . . . . . . . . . 95 4.8 The properties of sphere airdrop at Re = 1,000. (a): the sphere dropping velocity; (b): the drag force on the sphere; (c): fluid velocity near the ball at terminal state. . . . . . . . . . . 95 4.9 The vortex structure near the sphere at the terminal stage for Re = 1,000. (a) for vorticity at direction normal to the paper; (b): contour of vorticity magnitude. . . . . . . . . . . . . . . 96 4.10 The geometry of the PTV capsules for CPAS airdrop test . . . . . . . . . . . . . . . . . . . 97 4.11 The orientation or the PTV and the fluid velocity nearby. The snapshot was taken at a time when the wobbling angle is a local maximum. . . . . . . . . . . . . . . . . . . . . . . . . . 98 4.12 The kinematics and dynamics with time in the PTV airdrop simulations. (a): the x-velocity; (b): the aerodynamical drag in x-direction; (c): the pitch angles. . . . . . . . . . . . . . . . 99 5.1 (a): CPAS drogue parachute constructed geometry [6, 7]. (b): the typical construction of a conical parachute, the gore consists of multiple ribbons from the skirt to the vent [8]. . . . . 102 5.2 The initial drogue shape used in this chapter. (a) and (b): the top and side view of the drogue parachute. The 52 ribbons are grouped into four larger bands in the simplified model. (c) and (d): the riser line (yellow), suspension lines (blue), reefing lines (green), skirt/vent bands and radial tapes (purple) of the drogue. . . . . . . . . . . . . . . . . . . . . 103 5.3 The initial setup of two drogue parachutes and the PTV payloads . . . . . . . . . . . . . . . 104 5.4 triangular subdivision surface FE mesh of the CPAS drogue parachutes cluster . . . . . . . . 105 vii 5.5 Simulation schedules of stages 1, 2, and 3. PS[1,2,3] and S[1,2,3] denote pre-stage [1,2,3] and stage [1,2,3] simulations conducted in order, respectively. FC[1, 2, 3] denotes the atmospheric and flight conditions in Table 5.3 for each stage. . . . . . . . . . . . . . . . . . 107 5.6 Initial conditions generated of S1 generated by PS1. (a): the initial flow field velocity. (b): the semi-folded reefed drogue cluster configurations . . . . . . . . . . . . . . . . . . . . . . 109 5.7 The parachute related histories in the Case S1. (a): the average skirt circumference of two drogue parachutes in the cluster. (b): tensions of the riser line at gravitional direction. (c): scatter plot of the riser tension direction at sampling steps. . . . . . . . . . . . . . . . . . . 110 5.8 Snapshots of inflation at initial deployment from semi-folded configurations with von Mise stress of canopy displayed. (a): at the release moment. (b): 0.0202 s after release. (c): 0.0402 s after release. (d): 0.0599 s after release. . . . . . . . . . . . . . . . . . . . . . . . 111 5.9 The histories of the kinematics and dynamics of the PTV payload. (a): the total drag force on the gravitional direction. (b): the velocity along the gravitional direction. (c): the yaw angle α of the PTV. (d): the pitch angle β of the PTV. . . . . . . . . . . . . . . . . . . . . . 112 5.10 The histories of PS2 and S2 simulation before and after transition to reefing stage 2. (a): the average skirt circumference of two drogue parachutes in the cluster. (b): the tension of the riser line at gravitational direction. (c): scatter plot of the riser tension direction at sampling steps. Red and blue points are samples before and after disreefing (Case PS2 and S2). . . . . 113 5.11 Snapshots of inflation at the disreefing from reefed stage 1 to reefed stage 2 with von Mise stress of canopy. (a): at the starting moment. (b): 0.0084 s after disreefing. (c): 0.0169 s after disreefing. (d): 0.0254 s after disreefing. . . . . . . . . . . . . . . . . . . . . . . . . . 114 5.12 The histories of the kinematics and dynamics of the PTV payload. (a): the total drag force on the gravitational direction. (b): the velocity at gravitional direction. (c): the yaw angle α of the PTV. (d): the pitch angle β of the PTV. . . . . . . . . . . . . . . . . . . . . . . . . . 115 5.13 The histories of PS3 and S3 simulation before and after transiting to full open stage 3. (a): the average skirt circumference of two drogue parachutes in the cluster. (b): the tension of the riser line at gravitational direction. (c): scatter plot of the riser tension direction at sampling steps. Red and blue points are samples before and after completely removing the reefing line (Case PS3 and S3), respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . 116 5.14 Snapshots of inflation after reefing line is removed. (a): at the starting moment. (b): 0.0159 s after disreefing. (c): 0.0315 s after disreefing. (d): 0.0629 s after disreefing. . . . . 117 5.15 The histories of the kinematics and dynamics of the PTV payload. (a): the total drag force on the gravitational direction. (b): the velocity at gravitional direction. (c): the yaw angle α of the PTV. (d): the pitch angle β of the PTV. . . . . . . . . . . . . . . . . . . . . . . . . . 118 viii Acknowledgements I would like to express my deepest gratitude to my advisor, Prof. Carlos Pantano-Rubino, for his unwavering guidance, encouragement, and patience throughout these past years. His invaluable insights and active involvement were instrumental in overcoming numerous challenges, making this thesis and related publications possible. I will forever treasure the wisdom and support he has imparted. I am also grateful to my qualifying exam and dissertation committee members, Prof. Ivan BermejoMoreno, Prof. Paul Newton, and Prof. Roger Ghanem, for their constructive feedback and valuable suggestions that significantly enhanced the quality of this thesis and my other publications. I deeply appreciate the collaboration and camaraderie shared with the current and former members of my research group at USC and UIUC: David Petty, Aditya Heru Prathama, Timothy Smith, Ben Shield, Xiaoyi Lu, Sicong Wu, Qingquan Wang, and Junhao Ma. Our discussions have contributed directly or indirectly to this work. In particular, I am profoundly grateful to Qingquan Wang, whose research ideas and direct involvement in coding were essential to the completion of this work. I cherish the friendships with fellow Ph.D. students sharing the office: Anguo Hu, Weiye Wang, Kyu Ho Van, Suhyeon Oh, and Nami Mogharabin. Their support in my work and personal life made my Ph.D. journey all the more rewarding. To Yuan Meng and Yannan Li, I am sincerely thankful for the companionship and countless shared activities that brightened my time here. This research was generously funded by the NASA Early Stage Innovation grant and NASA John Space Center grant. I would like to express my sincere appreciation to Jessica Powell and Jared Daum for their technical expertise and detailed support in this work. Lastly, I extend my deepest thanks to my parents for their unconditional love and unwavering support throughout my academic journey. Their encouragement from afar made the completion of my doctoral degree possible. ix Abstract We present a projection immersed boundary (IB) method for simulating compressible viscous flow around solid walls without body-conforming grids. The immersed boundary is represented with an arbitrary unstructured surface grid, and the fluid dynamic is resolved with a Cartesian grid. The Eulerian and Lagrangian grids are not conformal, and convolutions with discrete delta functions are used to transfer data between the Eulerian and Lagrangian grids. Singular source terms supported on the Lagrangian grid of the immersed boundary, e.g., boundary force vectors, are added to the governing equations of the fluid to enforce the desired boundary conditions. The resultant equation set is a system of differential-algebraic equations (DAE), where the singular source terms are algebraic variables, and the boundary constraints are algebraic constraints. The DAE system is solved with a half-explicit Runge-Kutta method, and the algebraic equation set of each Runge-Kutta stage is solved with a fractional step (projection) method. Various kinematic and thermal boundary conditions will be discussed under the framework of the projection IB method. The projection IB method, however, fails to produce accurate boundary forces that converge to the actual physical values when the grid is refined. This happens because the singular source terms are solutions of discrete Fredholm integral equations, which are not well-posed because the discrete integral operators do not have a bounded inverse as the grid refines. In most cases, the boundary source terms show spurious oscillation. To alleviate this issue, we proposed a generalized Tikhonov regularization approach that reformulates the problem into a minimization problem. The graph Laplacian is used to formulate the seminorm similar to the Sobolev H 1 seminorm that measures the roughness of the solution and serves as the penalty term in regularization. Two- and three-dimensional numerical tests are conducted to demonstrate accuracy. The projection IB method with regularization is implemented under a parallelized patch-based adaptive mesh refinement (AMR) framework for efficient computation of large-scale problems, with flow features and boundary proximity as refinement criteria. The fluid dynamic solver featuring IB and AMR is coupled with a thin-shell finite element Kirchhoff-Love plate solver as a loosely coupled fluid-structure interaction x (FSI) solver where the time integrator is done in a staggered manner. The FSI solver is implemented in noninertial frame of reference for generality. In the end, the FSI solver is used to simulate the drogue parachute cluster of the Orion Capsule Assembly System (CPAS). We focus on the transient dynamic of parachute deployment and dis-reefing at typical flight conditions within the flight envelope. Preliminary results will be reported. xi Chapter 1 Introduction 1.1 Motivation and background Parachutes are the most common decelerators used for atmospheric entry of spacecrafts. Parachutes decelerate, provide specific descent rates, provide stability and proper height/timeline for the payload, and serve as pilot parachutes in more complex parachute systems [9, 10]. The inflation, deceleration, and descent of parachutes involve complex multi-physical phenomena, including the large, rapid deformations of the canopy coupled with a turbulent wake and the unsteady separation flow near the deforming canopy [11]. Understanding these fluid-structure interaction phenomena is beneficial for parachute designs, ensuring that appropriate aerodynamical drags are provided, peak load and stress are below parachute structure tolerance, and fly-over angles do not exceed the design maximums. The fluid-structure interaction of high-speed parachutes has been studied experimentally in the early days of space exploration. In the pioneering work of Maynard [12], various types of decelerators were tested with Mach numbers between 1.6 and 3.0 and altitudes between 15 km and 37 km. A series of airdrop experiments were conducted at the same time and proved that parachutes are viable solutions for decelerators even for supersonic flow; see [9] for the summary. Until today, wind tunnels, mortar shots, and field tests are still required in parachute designs and verifications. Recent Orion Capsule Parachute System (CPAS) experimental studies can be found in [6, 13, 7]. The computational simulations of the fluid-structure interaction of parachutes have gained popularity with increased computational powers. It is not intended to replace the experiments, but it will be helpful in narrowing the parameters in preliminary designs and reducing the number of expensive airdrop tests. Karagiozis [14] studied the dynamic behavior of the disk-gap-band (DGB) parachute, which is a legacy 1 of the Viking mission and typically used in low-density and supersonic environments, at Mach numbers between 1.96 and 2.48. Huang [15] investigated the same DGB parachute at supersonic regimes and focused on its inflation process. Takizawa [16] investigated the drogue parachute and studied the transient dynamics during dis-reefing. This thesis will study the transient dynamics of inflation and disreefing under non-inertial frames of reference using different numerical methods for computations. FSI simulations require high-fidelity numerical schemes for fluid and parachute canopies that are capable of dealing with large and rapid deformations, as well as stable and efficient coupled time integrators for the system. On the fluid side, the immersed boundary (IB) method is used in most parachute FSI simulations. The IB method enables flow simulations around moving solid bodies without dynamically generating body-conformal grids and is preferred for parachute problems where fabrics undergo large deformations. Peskin developed the original IB method [17, 18] for cardiovascular problems, and many variants were later proposed for various applications [19, 20, 21, 22]. See [23, 24, 25, 26, 27] for reviews. Karagiozis [14] used the ghost fluid method, a variant of IB methods, for supersonic parachute FSI simulations. In recent years, the projection IB method [28] has gained popularity after being proposed in 2007, owing to its elegant formulation without ad hoc parameters and easy implementations. The projection IB method has been adopted for FSI problems [29, 30, 31]. The method is designed for incompressible flow with no-slip boundary conditions. To apply the projection IB method for the parachute FSI problems, we extended the projection IB method to compressible viscous flow and introduced a method of enforcing shear stress and heat flux boundary conditions. In addition, a generalized Tikhonov regularization is introduced to accurately compute the boundary sources in the projection method. The IB method is generally used with adaptive mesh refinement (AMR) to control the grid spacing dynamically based on the flow features and geometric proximity without pre-generated grids, especially for flow with high Reynolds numbers where the near-wall grid should be significantly denser to resolve or model the dynamics in the near-wall layer. Both octree-based and patch-based AMR approaches are widely used with IB methods and for parachute FSI simulations, e.g., [32] and [14] for both types of AMR and variants of ghost fluid IB methods. However, limited literature combines the projection IB method with the AMR due to the complexity of data structures associated with the hierarchy Eulerian grid and Lagrangian grid: the thesis will present an efficient algorithm that combines the patch-based structured AMR and the projection IB methods. In addition, a partitioned loosely coupled FSI time integrator is introduced in this thesis as an efficient yet accurate coupling scheme. 2 Figure 1.1: CPAS systems and deployment sequence [3] The Orion Capsule Parachute Assembly System (CPAS) is the first human-rated capsule parachute system since the Apollo program in the U.S. [10, 3]; part of the Artemis Program. The first generation CPAS airdrop test was conducted between 2007 and 2008, and the qualification test for the final CPAS system was conducted in 2018; see [3] for the summary of all airdrop tests. The final CPAS system consists of three mortar-deployed forward bay cover parachutes, two mortar-deployed drogue parachutes, three mortardeployed pilot parachutes, and three pilot-deployed main parachutes. They are deployed in sequence in flight missions as illustrated in Figure 1.1. The two small drogue parachutes are the decelerators and stabilizers before the main parachutes are deployed. The system of two drogue parachutes will be simulated using the IB method and the AMR framework as the final problems, and we will focus on the inflation and disreefing process. 1.2 Objectives and thesis outline This thesis mainly aims to propose numerical methods and computational frameworks for fluid-structure interaction (FSI) problems with reasonable computational cost and acceptable accuracy. Preliminary simulation results for the deployment and descent of drogue parachutes will be included to validate numerical methods and provide direction for future studies. 3 The thesis is organized as follows: Chapter 2 presents the projection immersed boundary method for compressible viscous flow. The immersed boundary method is the foundation element that couples the fluid to the structure in our FSI solver. By removing the necessity of body-conforming mesh, this special numerical method simplifies the process of resolving the fluid dynamic when the structure shows large and rapid deformations. Chapter 3 presents an improvement of the projection IB method that allows accurate computations of boundary sources and relaxes the stringent requirements of the fluid/structure grid spacing relation. Chapter 4 introduces the computational frameworks of the FSI, including the implementation of the fluid solver under the patch-based structure AMR and the partitioned loosely coupled FSI time integrator. In Chapter 5, the drogue parachute cluster in the CPAS will be studied numerically, focusing on the dynamics during inflation and disreefing. A detailed introduction and literature survey will be included in each chapter to discuss the topic more relevantly. Finally in Chapter 6, we offer some concluding remarks on the numerical methods and the parachute FSI simulations in general with some possible future research ideas. 1.3 Related publications The thesis is based on the following published work by the author: • Hang Yu and Carlos Pantano. An immersed boundary method with implicit body force for compressible viscous flow. Journal of Computational Physics, 459:111125, June. 2022. • Hang Yu and Carlos Pantano. A regularized projection immersed boundary method for smooth boundary forces. Journal of Computational Physics, 496:112571, Jan. 2024. 4 Chapter 2 The projection immersed boundary method for compressible viscous flow 2.1 Immersed boundary method: a literature survey The immersed boundary (IB) method is a family of methods in computational fluid dynamics that enables flow simulations around solid bodies without body-conformal grids. The IB method avoids dynamic remeshing for moving boundary problems. The original IB method was developed by Peskin [17] in the 1970s for cardiovascular simulations by introducing boundary forces as source terms of the Navier-Stokes equation and utilizing discrete delta functions to transfer data between the fluid and immersed bodies. The IB method was primarily used by Peskin and his collaborators for cardiovascular problems until the 1990s, when it gained wider recognition and numerous variants began to be published. The IB method is dichotomized into the continuous forcing approach and the discrete forcing approach, according to the classification of Mittal [23, 27]. In the continuous forcing approach, the governing equation in its continuous form was modified to account for the presence of the immersed boundary through added boundary force terms, and no special treatments near the IB are needed in numerical discretization. The continuous forcing methods include Peskin’s original work [17] and some subsequent variants, e.g., [19, 33, 34, 28, 35]. In the discrete forcing approach, the boundary condition is enforced by adjusting the spatial discretizations near the IB, modifying the computational stencils, or altering the solution data. Examples of discrete forcing methods include the ghost cell method [36, 37, 38] and the cut-cell method [39, 40]. Per the classification in [24, 25], the continuous forcing methods belong to the family of diffused interface methods since the boundary interface is smeared to several fluid grid points using continuous kernel functions. The discrete forcing methods generally overlap the family of sharp interface methods, where the boundary is represented without smearing effect and thus is desired for flows with moderate to high Reynolds numbers. 5 In the original method proposed by Peskin [17], the Navier-Stokes equation is modified to add the elastic force that the object exerts on the fluid, i.e., ρ ∂u ∂t +u·∇u = −∇p+ µ∆u+F, (2.1) where u(x,t) is the velocity vector and p(x,t) is the pressure. F(x,t) is the external force function from the body acting on the fluid. Eq.(2.1) is valid in the whole computational domain, including the volume occupied by the solid body, and therefore, numerics can be carried out irrespective of the presence of the body. The body force F(x,t) comes from the interface of the fluid and solid and is written formally as F(x,t) = Z Γ(t) f(s,t)δ(x−xbc(s,t))ds, (2.2) where f(s,t) is the force function on the possibly moving boundary of the immersed body, parametrized by s with coordinate xbc(s,t), and δ is the multi-dimensional Dirac delta. The integral (2.2) should be understood as a surface singularity [41] and interpreted in the sense of distributions, i.e., ⟨F,φ⟩:= R Γ(t) f(s,t)φ(xbc(s,t))ds for any regular enough test function φ(x). The numerical approximation of the operator f 7→ F spreading the boundary force from the solid body to fluid is achieved with discrete delta functions in most IB methods. In the works of Peskin [17, 18, 42, 43], the body force f(s,t) is calculated with the body’s material constitutional relation (Hooke’s law for elastic objects) as a functional of the deformed configurations xbc. The xbc(s,t) is conversely evolved in time by the fluid velocity, i.e., dxbc dt = u(xbc(s,t),t) = Z Ω u(x,t)δ(xbc(s,t)−x)dx, (2.3) where the discrete delta function serves as an interpolation kernel. The interpolation is numerically discretized as the adjoint of the spreading operator in (2.2). Peskin’s original method assumes that the solid body is massless and the interface force equals the body’s internal force. On the fluid computation, no kinematic (velocity) boundary condition is imposed. Therefore, it is limited to a small set of fluid-structure interaction (FSI) problems with small solid-to-fluid density ratios and unsuitable for immersed boundaries with prescribed motion or FSI problems where solidto-fluid density ratios are not small. To impose a kinematic condition on the immersed boundary, Goldstein 6 [19, 33] proposed a feedback PI controller method, where each point of the immersed boundary is treated as a spring-damper with large negative stiffness α and damping ratio β, i.e., f(s,t) = α Z t 0 u(xbc(s,t ′ ),t ′ )−x˙bc(s,t ′ ) dt ′ +β (u(xbc(s,t),t)−x˙bc(s,t)) (2.4) where u(xbc(s,t),t) is the controlled variable and the prescribed velocity x˙bc(s,t) is the target value. This method requires a small time step due to large α and β. Modh-Yosuf [20] proposed the direct forcing method without time step restrictions. The boundary force is calculated only during numerical discretization, and the temporal discretized momentum equation has the form [20, 21] u n+1 −u n ∆t = rhsn+1/2 +f, (2.5) where rhsn+1/2 lumps the convection, pressure, and diffusion terms at some intermediate time level, and the boundary force f is non-zero only at selected grid points and is calculated from f = u (d) −u n ∆t −rhsn+1/2 , (2.6) at these grid points, where u (d) is the desired velocity vector. By adopting the discrete delta functions as transferring kernel, the direct forcing method can be used with arbitrary non-conforming grids, see [22, 44, 45]. Detailed reviews and comparisons of the direct forcing method, the feedback forcing method, and other variants can be found in [21]. Another approach treats the boundary force as the Lagrange multipliers for the no-slip constraint [46, 34]. The method computes the boundary force without ad hoc parameters, e.g., spring damping. Taira [28] noticed that the distributed Lagrange multipliers method leads to a system that is algebraically identical to the incompressible Navier-Stokes equation and presented the IB method as a projection approach. In fact, it is a system of Hessenberg index-2 differential-algebraic equations (DAE); see [47, 48] for more descriptions. The boundary force is the algebraic variable, and the velocity condition is the constraint. The fractional step method is generally used for the time integration of the index-2 DAE system. So far, most IB methods focus on no-slip conditions, and they do not have similar treatments for the shear stress and/or heat flux boundary conditions. Our intention is to extend the projection approach to compressible viscous flow problems for both fixed bodies or bodies with prescribed motions and various kinematic and thermal boundary conditions. 2.2 Introduction Here, we formulate the projection immersed boundary method for compressible viscous flow problems. Similar to the many previous researches based on Lagrange multipliers, singular source terms supported on the immersed boundary are introduced into the governing equations in conjunction with algebraic constraints, and the equation set becomes a system of differential-algebraic equations. Our interest goes beyond the no-slip and isothermal boundary, and we will also focus on the slip wall and adiabatic conditions. In the slip wall and/or adiabatic condition, there are possible discontinuities in tangential velocity and/or temperature fields across the immersed boundary because they are not set to fixed values. Their jump across the immersed boundary is unknown a priori, and it will result in singularities in stresses or heat fluxes on the immersed boundary that need to be canceled. To achieve this, we identify singular source terms in the stress and/or heat flux in addition to the momentum and energy equation. In addition, stress and/or heat flux conditions are incorporated to calculate the singular source terms with a similar projection approach. All singular source terms are identified rigorously with jump conditions, and the temporal discretization is achieved by the half-explicit Runge-Kutta method. This chapter is organized as follows: first, in Section 2.3, we formulate the governing equations in continuous form for the projection immersed boundary method of compressible viscous flow for two typical boundary conditions. The governing equations are systems of differential-algebraic equations, with the singular source terms being the algebraic variable and the constraints being the algebraic equation. In Section 2.4, spatial and temporal discretization of the governing equations are introduced. The spatial discretization part mainly focuses on approximating the singular source terms and the algebraic constraints, and the temporal discretization part focuses on the fractional step (projection) method for DAE systems. In Section 2.5, we extend the IB method to permeable interfaces that follow Darcy’s law and present preliminary results. In Section 2.6, numerical examples and results with our method will be presented in twoand three-dimensional problems with various Mach and Reynolds numbers. We will compare some of these results with other literature. Spatial and temporal convergence tests are included. 8 2.3 Governing equations of compressible projection IB methods This section will formulate the governing equation in continuous form for the projection immersed boundary method for compressible viscous flow. Specifically, we focus on the no-slip isothermal condition and the slip adiabatic condition. Other boundary conditions, e.g., the no-slip adiabatic condition, can be derived following the same path. 2.3.1 Singular source terms and jump conditions The boundary force vector is the only singular source term for incompressible problems with the no-slip wall condition. More boundary source terms are needed for compressible flow, and the jump condition will identify them. We denote the computational domain of the fluid as Ω and the immersed boundary as Γ ⊂ Ω. The original compressible Navier-Stokes equation for d-dimensional problems in its conservative form is valid in Ω\Γ where the immersed boundary is excluded, i.e., ∂U ∂t + d ∑ k=1 ∂F k ∂ xk = 0 in Ω\Γ, (2.7) where the conservative variables U and flux functions F = [F 1 ,··· ,F d ] at each direction are U = ρ ρu ρe , F = F 1 ,··· ,F d = ρu T ρu⊗u T + pI−σ (ρe+ p)u T −u Tσ +q T , (2.8) where ρ is the density, u = [u1,··· ,ud] T is the velocity vector and e is the total specific energy and p is the pressure. The total specific energy and the temperature for the ideal gas are e = p ρ(γ −1) + 1 2 |u| 2 , T = p ρR , (2.9) 9 where γ is the specific heat ratio and R is the specific gas constant. The stress tensor σ = {σi j} and heat flux vector q = [q1,··· ,qk] T are σ = µ ∇u+∇u T − 2 3 (∇·u)I , q = −κ∇T, (2.10) where µ is the dynamic viscosity and κ is the thermal conductivity. Note that Eqs. (2.7)-(2.10) are only valid in Ω\Γ where the immersed boundary Γ has been excluded. Appropriate modifications to these governing equations are needed to extend the applicability to the whole domain Ω so that the immersed boundary method can be naturally applied. To derive the governing equations that are globally valid in Ω, we formally extend the domain of the solution function U(x,t) from Ω\Γ to Ω and denote the extended solution function as Uˆ (x,t). The extended solution Uˆ (x,t) = U(x,t) everywhere for x ∈ Ω \ Γ. In fact, we do not really need to assign the values of Uˆ (x,t) on the zero-measure surface Γ. However, the extended Uˆ (x,t) is not differentiable because it is generally discontinuous across Γ, and the spatial and temporal derivative of Uˆ (x,t) have to be taken in the sense of distributions. Similarly, we extend the fluxes, stresses, and heat fluxes function to Ω and denoted them by Fˆ, σˆ , and qˆ, respectively. Generally, taking the derivative in the sense of distributions for a function with discontinuity across a surface leads to surface singularity; see Appendix.A for the derivations. Applying the temporal derivative to Uˆ and spatial derivative to Fˆ in the sense of distributions gives ∂Uˆ ∂t = ∂U ∂t − Z Γ JUˆ K(s,t)(x˙bc(s,t)·nˆ(s,t))δ(x−xbc(s,t))ds, (2.11) ∂Fˆ k ∂ xk = ∂F k ∂ xk + Z Γ JFˆ k K(s,t)nk(s,t)δ(x−xbc(s,t))ds, (2.12) where J·K(s,t) denotes the jump of a variable across the surface. xbc(s,t) and x˙bc(s,t) denote the position and velocity of the boundary, and nˆ(s,t) = [nˆ1,··· ,nˆd] denotes the unit normal vector of the surface. The convolution with δ(x) over surface Γ should be understood as representing surface singularity, and the added singular source terms cancel the singularity arising from taking derivative across Γ. By comparing the equation against the original N-S equation (2.7), we reach the immersed boundary equation for compressible viscous flow: ∂Uˆ ∂t + d ∑ k=1 ∂Fˆ k ∂ xk = Z Γ φ S (s,t)δ(x−xbc(s,t))ds, (2.13) 10 where the added source term is φ S (s,t) = d ∑ k=1 JFˆ kK(s,t)nˆ k(s,t)−JUˆ K(s,t) (x˙bc(s,t)·nˆ(s,t)) = Jρˆ (uˆ −x˙bc) T Knˆ Jρˆuˆ ⊗(uˆ −x˙bc) T + pˆI−σˆ Knˆ J ρˆ eˆ(uˆ −x˙bc) T + pˆuˆ T −uˆ Tσˆ +qˆ T Knˆ , (2.14) which is a vector-valued function defined on the immersed boundary Γ(t). Similarly, applying the spatial derivative to uˆ and Tˆ in the sense of distributions gives surface singularity related to the jumps Juˆ K and JTˆ K across the boundary. Therefore, the Eq. (2.10) with the immersed boundary method becomes σˆ = µ ∇uˆ +∇uˆ T − 2 3 (∇·uˆ)I − Z Γ φT (s,t)δ(x−xbc(s,t))ds , (2.15) qˆ = −κ ∇Tˆ − Z Γ φQ(s,t)δ(x−xbc(s,t))ds , (2.16) and the added singular source terms are φT = JuK⊗nˆ +nˆ ⊗JuK− 2 3 JuK T n I ˆ , φQ = JTKnˆ, (2.17) which are d ×d dimensional tensor function and d dimensional vector function on Γ, respectively. Eqs.(2.13)-(2.17) with added singular sources are valid throughout Ω. The added singular source terms cancel the singularities that arise when taking the derivative on the discontinuity interface. Herein, all allowed singular source terms for the compressible immersed boundary method have been identified, and they are valid for all types of immersed boundary conditions. The singular source terms can be simplified for specific conditions, knowing that the jumps of specific flow variables are zero. The simplification for the no-slip isothermal condition and the slip adiabatic condition will be discussed in the next section. For clarity of presentation, the hat over the flow variables is subsequently dropped, and all flow variables are interpreted as valid throughout Ω. 11 2.3.2 Differential-algebraic system of governing equations We have derived all possible singular source terms that must be added to the original Navier-Stokes equation, the laws of viscosity, and the laws of heat conduction to account for the presence of the immersed boundary. However, the equations are incomplete because the singular source terms are unknown a priori, and specific kinematic and thermal boundary conditions must be introduced as constraints to achieve closure. Here, we will derive the systems of governing equations for two staple boundary conditions, namely, the no-slip, isothermal condition, and slip, stress, and heat flux condition. Other boundary conditions, such as no-slip and adiabatic, can be trivially derived with the same approach. For clarity, we first project the velocity vector, surface stress vector, and heat flux vector along the normal and tangential directions of the surface according to u ±(s,t)−x˙bc(s,t) = u ± n (s,t)nˆ(s,t) +u ± τ (s,t) (2.18) σ ±(s,t)·nˆ(s,t) = σ ± n (s,t)nˆ(s,t) +σ ± τ (s,t) (2.19) q ±(s,t) = q ± n (s,t)nˆ(s,t) +q ± τ (s,t) (2.20) (2.21) where nˆ(s,t) is the unit normal vector of the surface pointing toward the positive surface direction. The ± denotes the values of the flow variables on each side of the surface. Noting that Ju − x˙bcK = JuK = JunKnˆ +Juτ K, JσK · nˆ = JσnKnˆ +Jστ K and JqK · nˆ = JqnK in terms of the new variables, the function φ S (s,t) in Eq.(2.14) can be rewritten as φ S = JρunK JρunuK+J(p−σn)nˆ −στ K JρuneK+Ju·((p−σn)nˆ −στ )K+JqnK , (2.22) and the tensor function φT (s,t) can be rewritten as φT = JunK(2nˆ ⊗nˆ − 2 3 I) +Juτ K⊗nˆ +nˆ ⊗Juτ K. (2.23) 12 2.3.2.1 No-slip and isothermal conditions Here, we consider the no-slip and isothermal conditions on the surface. The flow velocity vector at the boundary is u(xbc(s,t),t) ≡ x˙bc(s,t), which is the velocity of the immersed boundary itself. The temperature at the boundary is T(xbc(s,t),t) ≡ Tbc(s,t), where Tbc(s,t) is the desired temperature. Noting that JuK = 0 and JTK = 0 since the velocity and temperature are prescribed and are identical across Γ, and therefore, the boundary source terms φT and φQ are identically zero from Eq. (2.17). Also, from the velocity condition that u(xbc(s,t),t) ≡ x˙bc(s,t) the boundary source term φ S (s,t) can be simplified to φ S (s,t) = 0 Jpnˆ −σ ·nˆ K x˙bc · Jpnˆ −σ ·nˆ K+JqnK = 0 f g , (2.24) where f(s,t) = Jpnˆ −σ ·nˆ K is the boundary force vector and g(s,t) = x˙bc · Jpnˆ −σ ·nˆ K+JqnK is the singular boundary energy source, from both the surface friction x˙bc ·Jpnˆ −σ ·nˆ K and heat source JqnK. Using the same distributed Lagrange multiplier method that has been used for incompressible immersed boundary method, e.g., [28, 35], we treat the boundary force f and energy source g as Lagrange multipliers and boundary conditions as the constraints. The constraints are written as Z u(x,t)δ(xbc(s,t)−x)dx = x˙bc(s,t), (2.25) Z T(x,t)δ(xbc(s,t)−x)dx = Tbc(s,t), (2.26) where the convolution with δ(x)interpolates data from the fluid to the boundary, similar to other IB methods. The Navier-Stokes equation with singular source terms Eqs.(2.13)-(2.16) together with the algebraic constraints (2.25)-(2.26) forms a system of differential-algebraic equation after proper spatial discretization, and would be solvable using the projection method. The spatial and temporal discretization will be discussed in Section 2.4. 2.3.2.2 Stress and heat flux conditions Now, we consider situations where the shear stress and the heat flux are specified on the surface instead of the tangential velocity and temperature, and the normal velocity of the fluid at the boundary still equals 13 the normal moving velocity of the boundary, i.e., the fluid does not penetrate the immersed boundary. The no-penetration normal velocity condition u(xbc(s,t),t)· nˆ(s,t) = x˙bc(s,t)· nˆ(s,t) leads to un = 0. The shear stress vectors and heat fluxes at both sides of the surface are specified by σ ± τ,bc(s,t) and q ± n,bc(s,t). The zero-stress slip and adiabatic condition where σ ± τ = 0 and q ± n,bc is the most prominent condition of this kind, but allowing arbitrary wall shear stress and heat flux provides a possible framework for wall modeling for high-Reynolds number flows; this is the main limitation and challenge of the IB method [26], and it is beyond the scope of the chapter. It should be noted that the shear stresses and the heat fluxes are given at both sides of the surface according to the same normal vector nˆ. For example, q + n denotes the heat flux on the positive side that comes out of the surface, while q − n denotes the heat flux on the negative side that goes into the surface. With the normal velocity condition un = 0, the term φ S (s,t) is simplified into φ S = 0 Jp−σnKnˆ −Jστ K Ju·((p−σn)nˆ −στ )K+JqnK = 0 fnnˆ −Jστ,bcK x˙bc ·nˆ fn −u¯ τ · Jστ,bcK−ζ ·σ¯τ,bc +Jqn,bcK , (2.27) where Jστ,bcK = σ + τ,bc − σ − τ,bc, Jqn,bcK = q + n,bc − q − n,bc, and σ¯τ,bc = (σ + τ,bc + σ − τ,bc)/2 are known (specified), while the normal force vector fn = Jp−σnK is unknown and serves as the Lagrange multiplier. ζ (s,t) = Juτ K is the vector of tangential velocity jumps, and u¯ τ = (u + τ +u − τ )/2 is the average tangential velocity. In our application, u¯ τ is approximated by u¯ τ (s,t) ≈ Pτ Z u(x,t)δh(xbc(s,t)−x)dx, (2.28) where δh is the discrete delta function, and the convolution will be spatially discretized in the same manner as the interpolation in Eq.(2.25) or (2.26). The approximation is justified when the boundary is flat enough with small curvature compared to the grid width, which is generally true. Pτ : u(s)7→ u(s)−(u(s)·nˆ(s))nˆ(s) is the pointwise projection operator on the boundary into its tangential direction. Similarly, the singular source φT (s,t) and φQ(s,t) are simplified to φT = Juτ K⊗nˆ +nˆ ⊗Juτ K, φQ = θnˆ, (2.29) 14 where ζ (s,t) = Juτ K(s,t) and θ(s,t) = JTK(s,t). They are unknown and serve as Lagrange multipliers. The unknowns fn(s,t), ζ (s,t), and θ(s,t) serve as Lagrange multipliers to satisfy constraints at the boundary. The normal velocity constraint can be specified as before, according to Z u(x,t)δ(xbc(s,t)−x)dx ·nˆ(s,t) = x˙bc(s,t)·nˆ(s,t), (2.30) which is a scalar constraint equation that allows us to solve for the unknown surface force fn(s,t). If the shear stress and heat flux are identical on both sides, for example, in the slip and adiabatic boundary conditions, the shear stress and heat flux constraints are Pτ Z σ(x,t)δ(xbc(s,t)−x)dx ·nˆ(s,t) = στ,bc(s,t), (2.31) Z q(x,t)δ(xbc(s,t)−x)dx ·nˆ(s,t) = qn,bc(s,t), (2.32) where στ,bc ≡ σ ± τ,bc is the shear stress and qn,bc ≡ q ± n,bc is the heat flux. The right-hand sides of Eqs.(2.31)- (2.32) are zero at a slip (zero shear stress) and adiabatic boundary. However, for the more interesting cases where the desired shear stresses or heat fluxes are different on each side, i.e., σ + τ,bc ̸= σ − τ,bc or q − n,bc ̸= q + n,bc, the constraints above are no longer suitable. We propose to resolve this by constructing the target interface functions q int(x,t) and σ int(x,t) in the region surrounding the immersed boundary that approximates the real pattern of discontinuity of σ and q near the immersed boundary from σ ± τ,bc and q ± n,bc, respectively. Let the signed normal projected distance from point x to the immersed boundary Γ be d Π(x) = dist(x,Γ), and the projection point on Γ be s Π(x), the interface functions are defined by σ int(x,t) = σ − τ,bc +Jστ,bcKHh(d Π(x)) ⊗nˆ(s Π(x)), (2.33) q int(x,t) = q − n,bc +Jqn,bcKHh(d Π(x)) nˆ(s Π(x)), (2.34) where Hh(x) = 0 x < −2h, 1 4 2+ x h + 2sin( πx 2h ) π −2h < x < 2h, 1 x > 2h, (2.35) 15 is a mollified Heaviside function. The constraints are given by requiring the actual stress and heat flux to be similar to the constructed interface function, and we use the discrete delta δh as the test function, i.e., Pτ Z σ(x,t)δh(xbc(s,t)−x)dx ·nˆ(s,t) = Pτ Z σ int(x,t)δh(xbc(s,t)−x)dx ·nˆ(s,t), (2.36) Z q(x,t)δh(xbc(s,t)−x)dx ·nˆ(s,t) = Z q int(x,t)δh(xbc(s,t)−x)dx ·nˆ(s,t) = qn,bc. (2.37) The only difference between Eqs.(2.31)-(2.32) and Eqs.(2.36)-(2.37) after spatial discretization is the righthand side with the same discrete delta function. Note that technically speaking Eqs. (2.36)-(2.37) are not completely well defined in the sense of distributions because they involve products of Dirac delta and Heavyside function. The better treatment is to move the integrands into a single integral, e.g., (σ − σ int), before multiplying by the Dirac delta. This is not completely satisfying but give reasonable results in practice when the Dirac deltas has been approximated numerically. The constraints on the shear stress and heat flux allow us to solve for the unknown ζ (s,t) and θ(s,t), respectively. 2.4 Numerical discretization The flow governing equations will be approximated using a finite-volume discretization. The fractional step (projection) approach is used for the immersed boundary parts. The temporal marching strategy is similar to the projection method applied to incompressible flow [28, 35] but extended to other types of boundary conditions for compressible flow mentioned in the previous section. These elements of temporal and spatial discretizations and solution strategies are described in the following. 2.4.1 Spatial discretization We choose the finite-volume method as the spatial discretization of the compressible Navier-Stokes equation. Figure 2.1 shows a schematic diagram of the Eulerian grid of the fluid computational domain Ω and the Lagrangian grid of the immersed boundary Γ for two-dimensional problems. The domain Ω is uniformly discretized near the immersed boundary by a Cartesian grid where fluid data are associated. The 16 Eulerian grid cell center is positioned at xi, j = (xi , yj), and the faces in the x- and y-direction are positioned at xi+1/2, j = (xi+1/2 , yj) and xi, j+1/2 = (xi , yj+1/2 ), respectively, where numerical fluxes are calculated. The Lagrangian grid of Γ is represented by a set of Lagrangian markers {sα}α, positioned at xbc(sα). The grid layout trivially extends to three-dimensional problems with two-dimensional immersed surfaces. Figure 2.1: Two-dimensional schematic layouts of the finite-volume Eulerian grid of the fluid computational domain Ω and the Lagrangian grid of the immersed boundary Γ. With a finite-volume approach, the conservative law equation with singular source term (2.13) is spatially discretized to dU dt + d ∑ k=1 ∆F k ∆xk = Lφ S , (2.38) where U = Ui, j is the cell-averaged solution vectors and F 1 = n F 1 i+1/2, j o , F 2 = n F 2 i, j+1/2 o are the numerical flux vectors at x- and y-faces. ∆F k/∆xk denotes the difference quotient of the numerical fluxes at k-th direction where ∆F 1/∆x1 i, j = (F 1 i+1/2, j −F 1 i−1/2, j )/∆x1 and ∆F 2/∆x2 i, j = (F 2 i, j+1/2 −F 2 i, j−1/2 )/∆x2. In our numerical tests, the numerical inviscid flux F k is calculated with either skew-symmetric centered schemes [49], upwind schemes with WENO reconstructions and the Roe approximate Riemann solver, or a hybrid scheme that uses skew-symmetric centered schemes in isentropic flow regions without shocks and upwind schemes in regions of shock [50]. The operator L spreads the singular source from the Lagrangian grid to the Eulerian grids and will be introduced later. The viscous part of the numerical flux F k must be calculated with the stress tensor σ and heat flux vector q at Eulerian grid faces; they are evaluated according to Eqs.(2.15) and (2.16), where standard centered 17 finite-difference operators are used to calculate the velocity/temperature gradients with cell averaged solution data U. The singular source terms φT and φQ in Eqs.(2.15) and (2.16) are spread from the Lagrangian grid to the cell faces of k-direction of the Eulerian grid with operator L k . The singular source terms φ S , φT , and φQ on the IB-augmented Navier-Stokes equation (2.13), viscosity law (2.15), and heat flux law (2.16) are spread from the Lagrangian grid to the Eulerian grid by numerical convolutions with discrete delta functions, similar to previous IB work [17]. In the finite volume setting, the singular source term φ S shall be approximated on the grid cells, while singular source terms φT added to the stress function and φQ added to heat flux function shall be approximated on the grid faces. Herein, L denotes the spreading operator from the Lagrangian grid to the grid cells of the Eulerian grid, given by L(t)φ|i, j ≡ ∑α φαδ(xi, j −xbc(sα,t))∆sα, (2.39) for any discrete function φ = {φα} on Lagrangian grid points, and sα is the size of Lagrangian grid element α. The discrete delta δ(x) is defined as the product of one-dimensional discrete delta function [17], δ(x1,··· , xd) = d ∏ k=1 1 ∆xk b xk ∆xk , (2.40) where the width of the discrete delta equals the Eulerian grid spacing ∆xk. We choose the function b(r) given in [17]: b(r) = 1 4 1+cos πr 2 , |r| ≤ 2, 0, otherwise. (2.41) Similarly, L x and L y denote the spreading operators from the Lagrangian grid points to the grid faces of x- and y-direction of the Eulerian grids, achieved similarly by the numerical convolution with discrete delta functions L x (t)φ|i+1/2, j ≡ ∑α φαδ(xi+1/2, j −xbc(sα,t))∆sα, (2.42) and L y (t)φ|i, j+1/2 ≡ ∑α φαδ(xi, j+1/2 −xbc(sα,t))∆sα, (2.43) respectively. The same discrete delta function is used to construct the interpolation operator E from grid cells to the Lagrangian grid points, needed to discretize the algebraic boundary constraints, E(t)Ψ|α ≡ ∑ i, j Ψi, jδ(xbc(sα,t)−xi, j)∆x∆y, (2.44) for any discrete function Ψ = {Ψi, j} on Eulerian grid cells. Interpolation operators from the cell faces are E x and E y , defined similarly as E x (t)Ψ|α,x ≡ ∑ i, j Ψi+1/2, jδ(xbc(sα,t)−xi+1/2, j )∆x∆y, (2.45) or E y (t)Ψ|α,y ≡ ∑ i, j Ψi, j+1/2δ(xbc(sα,t)−xi, j+1/2 )∆x∆y, (2.46) for faces of x- and y-direction, respectively. These interpolation operators from cell faces are used to discretize the algebraic constraints of stress or heat flux defined on the cell walls. For clarity, we may omit the superscripts of L x , L y , E x , and E y later on since they are closely related to the underlying finite-volume scheme instead of the immersed boundary method itself. Note that the interpolation operator E is the transpose of the spreading operator L with scaling, and only one of the matrix E or L needs assembling. 2.4.2 Half-explicit temporal discretization The spatially discretized governing equations are a collection of differential equations in time and algebraic constraints, i.e., a system of differential-algebraic equations. The general structure of the system is of the form dU dt = − ∆F k (U,σ,q) ∆xk +L(t)φ S , (2.47) 0 = G(t,U), (2.48) where the constraints over solution variable U, e.g., the velocity or the temperature constraints, are lumped together and written in a generic form G(t,U). The stress tensor σ and heat flux vector q needed for the flux function F k are algebraic functions of the solution variable U independent of time. However, unlike the original N-S equation, additional singular source terms and constraints may need to be incorporated in 19 addition to the laws of viscosity and heat fluxes. The calculations of stress tensor from U are written in the generic form σ = σ(U)− µL(t)φT , (2.49) 0 = Gσ (t,σ), (2.50) where σ(U) is the original Newton law of viscosity composed by several differential operators, and Gσ is the generic form of the algebraic constraints over the stress tensor. Similarly, the heat flux is calculated from U with q = q(U) +κL(t)φQ, (2.51) 0 = Gq(t,q), (2.52) where q(U) is the original Fourier law of heat conduction, and Gq is the generic form of the algebraic constraints over the heat flux vector. The DAE system (2.47)-(2.48) can be advanced with first order explicit treatment of the fluxes, including the advection and diffusion, and implicit treatment for the algebraic constraint (2.48), this leads to a system that is solvable with the fractional step projection method, similar to the projection IB method for incompressible flow problems. The stress tensor and heat flux vector are calculated first from Eqs.(2.49)-(2.50) and Eqs.(2.51)-(2.52), respectively, for the numerical flux F k . These two algebraic systems are solved with a similar projection approach, although these two systems are not related to temporal marching. A natural extension of the first order method for higher accuracy is the half-explicit Runge-Kutta method [2, 51]. The half-explicit RK methods for DAE treat all terms explicitly except for the constraints in each stage with the same coefficients as the conventional explicit RK methods and require solving similar algebraic systems with the projection approach in each RK stage. For computational efficiency and stability, we use the half-explicit Runge-Kutta method based on the third-order optimal strong stability-preserving RungeKutta (SSPRK) scheme [1]. This method has second-order temporal accuracy and can be implemented in low-storage formats. Half-explicit Runge-Kutta methods with higher-order temporal accuracy can be constructed, but one needs more stages and might not be able to implement that in an efficient low-storage format, see [2]. 20 Applying the half-explicit RK method based on SSPRK3 to Eqs. (2.47)-(2.48) leads to the following time-marching scheme U˜ ν = U˜ ν,⋆ +γν∆tL(t˜ ν−1 )φ˜ ν−1 S , (2.53) 0 = G(t˜ ν ,U˜ ν ), (2.54) with the intermediate (star) state U˜ ν,⋆ = [ρ˜ ν,⋆ ,ρ˜ ν,⋆u˜ ν,⋆ ,ρ˜ ν,⋆e˜ ν,⋆] T computed first from U˜ ν,⋆ = ανU n +βνU˜ ν−1 −γν ∆t ∆xk ∆F k (U˜ ν−1 ,σ˜ ν−1 ,q˜ ν−1 ), (2.55) where t˜ ν = t n +cν+1∆t, U˜ ν and φ˜ ν−1 S denotes the temporary scratch solution variables and algebraic variables. Before calculating the numerical fluxes F k and the intermediate state U˜ ν,⋆, the temporary stress tensor σ˜ ν−1 and heat flux vector ˜q ν−1 need to be calculated first by solving the system arises from Eqs.(2.49)-(2.50) σ˜ ν−1 = σ(U˜ ν−1 )− µL(t˜ ν−1 )φ˜ ν−1 T , (2.56) 0 = Gσ (t˜ ν−1 ,σ˜ ν−1 ), (2.57) and from Eqs.(2.51)-(2.52) q˜ ν−1 = q(U˜ ν−1 ) +κL(t˜ ν−1 )φ˜ ν−1 Q , (2.58) 0 = Gq(t˜ ν−1 ,q˜ ν−1 ), (2.59) with known U˜ ν−1 , respectively. The time integration is accomplished in each step by solving the fully discretized Eqs.(2.53)-(2.59) successively for each stage ν = 1,2,3. For each time step, the iteration starts by setting U˜ 0 = U n and ends with U n+1 = U˜ 3 . The coefficients of the scheme are given in Table 2.1. Note that this scheme does not require the permanent storage of past intermediate values since there is no need to retain U˜ ν−1 once the ν-th stage is updated from U˜ ν−1 to U˜ ν . φ˜ 0 S , φ˜ 0 T and φ˜ 0 Q calculated at the first stage are the boundary source terms at the time step n, the beginning of the step, if needed. The solution strategy for Eqs.(2.53)-(2.54), Eqs.(2.56)-(2.57) and Eqs.(2.58)-(2.59) are achieved with the fractional step (projection) method. Detailed discussion for each type of boundary condition will be discussed below. 21 Table 2.1: Coefficients of the 3-stage SSPRK [1] formulated as a DAE integrator [2]. ν αν βν γν cν 1 0 1 1 0 2 3 4 1 4 1 4 1 3 1 3 2 3 2 3 1 2 2.4.2.1 No-slip and isothermal conditions Here, we discuss the concrete soluting strategy of Eqs.(2.53)-(2.59) for the no-slip and isothermal condition. The surface singular terms φT and φQ are zero for the no-slip and isothermal conditions, and there is no explicit constraint on stress and heat flux. Therefore, Therefore, σ˜ ν−1 and ˜q ν−1 can be calculated trivially from U˜ ν−1 . The only difficulty remains is solving Eqs.(2.53)-(2.54) to correct the intermediate state U˜ ν,⋆ to the final state U˜ ν . In the no-slip and isothermal condition, the expression of the source terms φ S comes from (2.24) and the constraints (2.25) and (2.26) are lumped together as G(t,U). Note that there are algebraic nonlinearities when casting conservative variables to primitive variables to enforce the velocity and temperature constraints. However, these nonlinearities are only apparent since the algebraic equation can be solved progressively from continuity to momentum and energy components. The solution strategy is described as follows: 1. The continuity equation: since there is no surface singularity added to the equation, as seen from the expression φ S in Eq.(2.24), we have: ρ˜ ν = ρ˜ ν,⋆; (2.60) 2. The momentum equation: with the momentum component of the singular source φ S in Eq.(2.24) and the velocity constraint (2.25), we reach the system I −γν∆t diag(1/ρ˜ ν )L(t˜ ν−1 ) E(t˜ ν ) 0 u˜ ν ˜f ν−1 = u˜ ν,⋆ x˙bc(t˜ ν ) , (2.61) where ρ˜ ν−1 at each Eulerian grid cell has already been calculated. The system is solved with a fractional step (projection) approach: first, the boundary force vector ˜f ν−1 is solved from the algebraic 22 equation that arises from applying the Schur complement to the system above, i.e., applying the constraint to the first equation in the equation set (2.61). The equation for ˜f ν−1 is E(t˜ ν )diag(1/ρ˜ ν )L(t˜ ν−1 ) ˜f ν−1 = 1 γν∆t (x˙bc(t˜ ν )−E(t˜ ν )u˜ ν,⋆). (2.62) Instead of solving the linear equation with a matrix-free implementation of a Krylov subspace method. The right hand side operator is constructed explicitly by sparse matrix multiplications. This is because the dimension of the equation is the number of Lagrangian grid points, while the dimension of the range of L and the domain of E is the number of Eulerian grid points, which is generally a larger number. Applying operators L and E sequentially as a composite operator to vectors is far more expensive than pre-multiplying explicit matrices E and L. With the boundary force, the momentum is corrected with a projection step ρ˜ ν u˜ ν = ρ˜ ν,⋆u˜ ν,⋆ +γν∆tL(t˜ ν−1 ) ˜f ν−1 ; (2.63) 3. Finally, the total energy is corrected using the energy component of singular source in Eq.(2.24) and constraint (2.26). The resulting system is: I −γν∆t diag(1/ρ˜ ν )L(t˜ ν−1 ) E(t˜ ν ) 0 e˜ ν g˜ ν−1 = e˜ ν,⋆ RTbc(t˜ ν ) γ−1 + 1 2 E(t˜ ν )(u˜ ν · u˜ ν ) , (2.64) where the density ρ˜ ν and velocity ˜u ν at each Eulerian grid cell have both been calculated. The system is again solved with the projection method. The algebraic equation for ˜g ν−1 is E(t˜ ν )diag(1/ρ˜ ν )L(t˜ ν−1 ) g˜ ν−1 = 1 γν∆t RTbc(t˜ ν ) γ −1 −E(t˜ ν )(e˜ ν,⋆ − 1 2 u˜ ν · u˜ ν ) , (2.65) and the energy is corrected with a final projection step after solving ˜g ν−1 ρ˜ ν e˜ ν = ρ˜ ν,⋆e˜ ν,⋆ +γν∆tL(t˜ ν−1 )g˜ ν−1 . (2.66) 23 2.4.2.2 Stress and heat flux condition Here, we discuss the solution strategy when the stress and heat flux condition is used. The surface singular strengths φT and φQ are generally not zero, and there are explicit constraints on stress and heat flux. They have to be solved by Eqs.(2.56)-(2.57) and Eqs.(2.58)-(2.59) first before Eqs.(2.53)-(2.55) since the stress and heat flux are required in calculating the numerical fluxes F k . The solution strategy is summarized as follows: 1. Solving for the stress tensor σ˜ ν−1 = {σ˜ ν−1 i j } and singular source vector ˜ζ ν−1 = { ˜ζ ν−1 i } according to discretized Eq.(2.29) and (2.31). The system is σ˜ ν−1 i j = σ ⋆ i j −diag(µ)L(Ni ˜ζ ν−1 j +Nj ˜ζ ν−1 i ), ∀1 ≤ i, j ≤ d, Pτ d ∑ j=1 NjEσ˜ ν−1 j = στ,bc(t˜ ν−1 ), (2.67) where Ni = diag(nˆi1,··· ,nˆiα,···) is the diagonal matrix of the i-th component of the normal vector at all Lagrangian grid points, and σ˜ ν−1 j is the stress vector by fixing the second index of the stress tensor to j. On the right hand side, σ ⋆ = {σ ⋆ i j} = σ(U˜ ν−1 ) is the stress calculated from U˜ ν−1 with the laws of viscosity without immersed boundary correction. The system is similar to (2.61) and (2.64), and they are both linear. It is, however, cumbersome to write in matrix form. Note that vector ˜ζ ν−1 is restricted so that it has only tangential components on the boundary, and similarly, the prescribed shear stress vector στ,bc(t˜ ν−1 ) shall not have component normal to the surface. The operator E and L are evaluated at time t˜ ν−1 . The system is solved with a similar projection approach: solving the singular source ˜ζ ν−1 from an algebraic equation first and then calculating the stress tensor σ˜ ν−1 with a projection step. The equation for ˜ζ ν−1 is derived from Eq.(2.67) by applying the second constraint to the first equation, i.e., Pτ d ∑ j=1 NjE diag(µ)L(N ˜ζ ν−1 j +Nj ˜ζ ν−1 ) = Pτ d ∑ j=1 NjEσ ⋆ j −στ,bc,(t˜ ν−1 ), (2.68) 24 which is a vector equation. With some commutation error by taking Nj and N out of operators E diag(µ)L and noting that PτN = 0, the equation can be simplified into d ∑ j=1 NjE diag(µ)LNj ! ˜ζ ν−1 = Pτ d ∑ j=1 NjEσ ⋆ j −στ,bc, (2.69) where on the right hand side, the intermediate stress tensor σ ⋆ is interpolated to the Lagrangian grid points, contracted to the surface traction vector, and the normal component of traction vector is removed by Pτ . With the singular source, the correction step for the stress tensor is σ˜ ν−1 i j = σ ⋆ i j −diag(µ)L(Ni ˜ζ ν−1 j +Nj ˜ζ ν−1 i ) ∀1 ≤ i, j ≤ d,; (2.70) 2. Solving for the heat flux vector ˜q ν−1 = {q˜ ν−1 i } and the singular source θ˜ ν−1 with the projection method according to the discretized version of Eq.(2.29) and constraint (2.32). The system is q˜ ν−1 i = q ⋆ i +diag(κ)LNiθ˜ ν−1 , ∀1 ≤ i ≤ d, d ∑ i=1 NiEq˜ ν−1 i = qn,bc(t˜ ν−1 ), (2.71) where q ⋆ = {q ⋆ i } = q(U˜ ν−1 ) is the heat flux vector calculated from U˜ ν−1 with law of heat conduction. The system is again algebraically similar to (2.61), and the equation for θ˜ ν−1 is d ∑ i=1 NiE diag(κ)LNi ! θ˜ ν−1 = qn,bc(t˜ ν−1 )− d ∑ i=1 NiEq⋆ i , (2.72) where the linear operator on the left hand side is similar to that in Eq.(2.69) For constant κ, we can take the constant diagonal scaling diag(κ) out and solve with the left hand side operator ∑ d i=1NiELNi . With θ˜ ν−1 , the heat flux is finally corrected with a projection step q˜ ν−1 i = q ⋆ i +diag(κ)LNiθ˜ ν−1 ∀1 ≤ i ≤ d; (2.73) 3. Calculating the numerical fluxes F k including the viscous flux according to the calculated σ˜ ν−1 and q˜ ν−1 . The remaining work is to solve Eqs.(2.53)-(2.54) progressively for continuity, momentum, and energy equations. 25 4. The continuity equation: no surface singularity added, and we have ρ˜ ν = ρ˜ ν,⋆; (2.74) 5. The momentum equation: with the source Eq.(2.27) and normal velocity constraint Eq.(2.30), we reach a system u˜ ν i = u˜ ν,⋆⋆ i +γν∆t diag(1/ρ˜ ν )L(t˜ ν−1 )Ni(t˜ ν−1 ) ˜fn ν−1 ∀1 ≤ i ≤ d, d ∑ i=1 Ni(t˜ ν )E(t˜ ν )u˜ ν i = x˙n,bc(t˜ ν ), (2.75) where ˜u ν,⋆⋆ = {u˜ ν,⋆⋆ i } = u˜ ν,⋆−γν∆t 1 ρ˜ ν L(t˜ ν−1 )Jστ,bcK(t˜ ν−1 ) is the intermediate state plus the tangential singular force known from prescribed Jστ,bcK(t˜ ν−1 ), and ˙xn,bc(t˜ ν ) is the prescribed normal velocity. The system is again similar to (2.61), and the equation for ˜fn ν−1 is: d ∑ i=1 Ni(t˜ ν )E(t˜ ν )diag(1/ρ˜ ν )L(t˜ ν−1 )Ni(t˜ ν−1 ) ! ˜fn ν−1 = 1 γν∆t x˙n,bc(t˜ ν )− d ∑ i=1 NiE(t˜ ν )u˜ ν,⋆⋆ i ! . (2.76) With the normal force ˜fn ν−1 , the correction step of momentum is ρ˜ ν u˜ ν i = ρ˜ ν,⋆u˜ ν,⋆⋆ i +γν∆tL(t˜ ν−1 )Ni ˜fn ν−1 ; (2.77) 6. The energy equation is corrected with the source without solving extra source terms since the singular energy source in (2.27) is known or has already been solved. The correction is done by simply adding this singular source after spreading to the Eulerian grid, i.e., ρ˜ ν e˜ ν = ρ˜ ν,⋆e˜ ν,⋆+γν∆tL(t˜ ν−1 ) x˙n,bc(t˜ ν ) ˜fn ν−1 −u˜¯ ν−1 τ Jστ,bcK(t˜ ν−1 )− ˜ζ ν−1σ¯τ,bc(t˜ ν−1 ) +Jqn,bcK(t˜ ν−1 ) . (2.78) 2.5 Preliminary studies on permeable surface Statement of the problem: The case of a permeable surface is considered here as an addition to demonstrate the versatility of the approach. This problem has received substantial attention when the flow is 26 incompressible, e.g., [52, 53, 54, 55, 56, 57]. In the compressible case, both density and velocity may be different on each side of the surface, and the mass flux rate ˙m = ρun is a function of the pressure change across the surface, which is assumed to be well described by the Darcy equation [58] for compressible flow with isothermal porous surface boundary condition, according to p +2 − p −2 2L = p¯ L JpK = −mRT ˙ bc µ k1 , (2.79) where L denotes the thickness of the membrane and ¯p = (p + + p −)/2, while k1 is Darcyan permeability parameter. Note that the velocity at each side of the surface is given by u ± = u ± n nˆ +x˙bc, with u ± n defined before as the penetration speed and u ± n = m˙ /ρ ±, the source φ S is reduced to φ S = 0 m˙ JuK+Jp−σnKnˆ −Jστ K m˙ JeK+Ju·((p−σn)nˆ −στ )K+JqnK = 0 m˙ JunK nˆ + fnnˆ +fτ g , (2.80) where we have used fn = Jp−σnK, fτ = −Jστ K and g = m˙ JeK+Ju·((p−σn)nˆ −στ )K+JqnK. Similarly, the source term φT and φQ are reduced to φT = JunK 2nˆ ⊗nˆ − 2 3 I , φQ = 0. (2.81) Simplification: We assume that the pressure jump is much larger than the normal stresses jump, |JpK| ≫ |JσnK|. Then, Eq.(2.79) reduces to fn = Jp−σnK ≈ JpK = −m˙ µ k1 LRTbc p¯ . (2.82) 27 Furthermore, we express JunK = m˙ J1/ρK = mRT ˙ bcJ1/pK = mRT ˙ bc −JpK p¯ 2− 1 4 JpK 2 ≈ mRT ˙ bc −fn p¯ 2− 1 4 f 2 n . Note that m˙ JunK fn = µ k1 p¯ L 2 1 RTbc f 2 n p¯ 2− 1 4 f 2 n , and for low permeability ˙mJunK ≪ fn. Therefore, we omitted JunK in both Eq.(2.80) and Eq.(2.81). The simplified source terms are φ S = 0 fnnˆ +fτ g ,φT = 0,φQ = 0. (2.83) As previously done, the equations are closed by enforcing velocity constraints, although the normal and tangential directions need to be separated. Since the mass flux is continuous across the surface, the normal velocity condition can be enforced by m˙(s,t) = nˆ(s,t)· Z ρ(x,t)(u(x,t)−x˙bc(s,t))δ(xbc(s,t)−x)dx, (2.84) while the tangential velocity continuity condition is Pτ Z u(x,t)δ(xbc(s,t)−x)dx = Pτx˙bc(s,t). (2.85) The temperature constraint is the same as Eq.(2.26). The average pressure ¯p is estimated with a similar approach in Eq.(2.28), p¯(s,t) = Z p(x,t)δ(xbc(s,t)−x)dx. (2.86) The spatial operators between the Eulerian and Lagrangian grids are the same as described in Section 2.4.1. Temporal marching: Compared to the spatially discretized systems in Section 2.4.2.1 and Section 2.4.2.2, the main difference here is that the normal force fn becomes an index-1 algebraic variable that appears in the constraint Eq.(2.48). Applying the similar half-explicit Runge-Kutta method leads to the temporal marching scheme: U˜ ν = U˜ ν,⋆ +γν∆tL(t˜ ν−1 )φ˜ ν−1 S , (2.87) (2.88) 28 where the intermediate stage U˜ ν,⋆ = ανU n+βνU˜ ν−1−γν ∆t ∆xk F k (U˜ ν−1 ,σ˜ ν−1 ,q˜ ν−1 ) is computed without the considering the immersed boundary. The algebraic variable φ˜ ν−1 S is calculated with a collection of algebraic constraints ˜m˙ ν = d ∑ i=1 NiE(t˜ ν )(ρ˜ ν u˜ ν i )−x˙n,bc(t˜ ν )E(t˜ ν )ρ˜ ν , (2.89) ˜m˙ ν = −ε ν−1 ˜fn ν−1 , (2.90) Pτx˙bc(t˜ ν ) = PτE(t˜ ν )u˜ ν , (2.91) Tbc(t˜ ν ) = E(t˜ ν )T˜ ν , (2.92) where the coefficient of permeability ε ν−1 = k1 LRTbc p¯ ν−1 µ . The temporal marching starts with U˜ 0 = U n and ends with U n+1 = U˜ 3 . The solution strategy for Eqs.(2.87)-(2.92) is as follows 1. No source term exists for the continuity equation, therefore ρ˜ ν = ρ˜ ν,⋆; (2.93) 2. The momentum equation has source terms of normal force ˜f ν−1 n and tangential force vector ˜f ν−1 τ = [ ˜f ν−1 τ,1 ,··· , ˜f ν−1 τ,d ] by examining φ˜ ν−1 S , the system of momentum source and kinematic constraints are summarized as follows: u˜ ν i = u˜ ν,⋆ i +γν∆t diag(1/ρ˜ ν )L(t˜ ν−1 )( ˜f ν−1 n Ni + ˜f ν−1 τ,i ), 1 ≤ i ≤ d ∑ d i=1NiE(t˜ ν )(ρ˜ ν u˜ ν i )−x˙n,bc(t˜ ν )E(t˜ ν )ρ˜ ν = ˜m˙ ν = −ε ν−1 ˜fn ν−1 , Pτx˙bc(t˜ ν ) = PτE(t˜ ν )u˜ ν . (2.94) The equations for ˜f ν−1 n and ˜fτ ν−1 are obtained from Eq.(2.94) E(t˜ ν )L(t˜ ν−1 ) + ε ν−1 γν∆t I ˜fn ν = x˙n,bc(t˜ ν )E(t˜ ν )ρ˜ ν −∑ d i=1NiE(t˜ ν )ρ˜ ν u˜ ν,⋆ i γν∆t , (2.95) E(t˜ ν ) 1 ρ˜ ν L(t˜ ν−1 ) ˜f ν−1 τ = 1 γν∆t Pτ (x˙bc(t˜ ν )−E(t˜ ν )u˜ ν,⋆); (2.96) 29 using the simplification that d ∑ i=1 NiE(t˜ ν )L(t˜ ν−1 ) ˜f ν−1 τ,i ≈ 0, d ∑ i=1 NiE(t˜ ν )L(t˜ ν−1 )Ni ≈ E(t˜ ν )L(t˜ ν−1 ). (2.97) Note that they are not identities because of the commutation error. In our application, however, the ˜fn ν−1 , ˜fτ ν−1 are actually solved by first solving Eq.(2.62) for the force vector, and then solving Eq.(2.95) to correct the normal components. After solving the boundary force, the intermediate velocity field ˜u ν,⋆ i is updated to ˜u ν i . 3. The energy equation is updated with the isothermal constraint. The procedure is identical to that in Section 2.4.2.1. 2.6 Results We numerically validate our compressible immersed boundary method with model problems that have been thoroughly studied, both with stationary and moving boundaries and at different Mach numbers. Temporal and spatial convergences of our method are shown in a two-dimensional problem of the flow around a ball with prescribed harmonic oscillations. The second-order skew-symmetric scheme [59, 49] is used for spatial discretization of the inviscid flux, and the second-order centered scheme is used for the viscous flux in most of the cases, except in the very compressible problems where the hybrid Roe-based Riemann fluxes were used. Other problems are also considered, and the results are compared with those in the literature. 2.6.1 Convergence tests of moving ball We first tested the temporal and spatial convergence of our compressible immersed boundary method with a two-dimensional flow around a moving ball. The domain is Ω = [−2,2]×[−2,2] with periodic boundary conditions in both directions to avoid possible boundary discretization and closure errors. The immersed boundary is a ball with radius R = 0.4 and moving in the y-direction with amplitude A = 1 and angular frequency ω = 40π (period Tp = 0.05), with surface coordinates x(s,t) = Rcoss Rsins+Acos(ωt) , 0 ≤ s < 2π. (2.98) 30 The flow field is initially at rest to be compatible with the velocity of the immersed boundary at t = 0, with ρ0 = 1.225 and p0 = 101,325; standard atmospheric conditions at sea level. The Reynolds and Mach numbers defined by the peak velocity of the ball and the Prandtl number are Re = ρ∞AωD µ = 11, Ma = Aω c = 0.11, Pr = cpµ κ = 0.7. (2.99) The first convergence tests are conducted with the no-slip, isothermal condition, where Tbc = p0/(ρ0R) = 288.144, and the second test uses the slip, adiabatic wall condition. -2 -1 0 1 2 -2 -1 0 1 2 Figure 2.2: Velocity field of moving ball problem at a t = 0.5 with no-slip and imposed temperature condition; internal induced flow inside the ball shown. Figure 2.2 shows the flow at one instant in time for the first case. Temporal convergence is conducted with a fixed Eulerian grid discretized with Nx = Ny = 50 cells in both directions and the immersed grid discretized by 28 uniformly spaced surface markers on the ball. Various time steps ∆t that all satisfy the CFL and diffusion number condition are used. The reference solution is computed with the same spatial grid and a small time step of 1.875×10−6 . The errors with different ∆t are computed by comparing the solution fields against the reference solution at t = 10 Tp. It should be noted that the initial condition of the flow should satisfy the constraint at t = 0 to observe the proper temporal convergence. Figure 2.3 shows the L1 norm of the error of different fields vs. the non-dimensionalized time step ∆t/Tp. We recover second-order temporal error, which is the expected temporal accuracy of the SSPRK3-based half-explicit Runge-Kutta method. 31 10-4 10-3 10-10 10-8 10-6 10-4 density velocity pressure 2 (a) 10-4 10-3 10-10 10-8 10-6 10-4 density velocity pressure 2 (b) Figure 2.3: L1 error norm of the moving ball problem non-dimensionalized by the reference fields as a function of ∆t: (a) no-slip and imposed temperature condition and (b) stress and heat flux condition; secondorder reference slop is plotted. Similarly, spatial convergence tests are conducted with various numbers of grid points from N = 25 to N = 300 in each direction and fixed time step ∆t = 5.0×10−6 (which is sufficiently small to satisfy the CFL and diffusion condition for the finest grid). The immersed boundary grid is refined in lockstep with the flow grid (at the same rate) to keep the spacing of the surface markers comparable to the grid spacing ∆x used to discretize the flow. The solutions are compared against a fine grid of N = 400 cells. Figure 2.4 shows the norms of the error of different flow fields normalized by the reference solution versus the grid width normalized by the diameter of the ball. Second-order spatial error convergence is observed for all fields. 10-2 10-1 10-4 10-3 10-2 10-1 100 2 density velocity pressure (a) 10-2 10-1 10-4 10-3 10-2 10-1 100 2 density velocity pressure (b) Figure 2.4: L1 error norm of the moving ball problem non-dimensionalized by the reference fields as a function of ∆x: (a) no-slip and imposed temperature condition and (b) stress and heat flux condition; secondorder reference slop is plotted. 32 A second set of convergence tests was conducted with no-slip, isothermal condition with time-dependent Tbc. The configuration is the same as the no-slip and imposed temperature test discussed before, except that the temperature Tbc increases linearly from 288.144 at t = 0 to 400 at t = 2Tp = 0.1. The numerical simulation uses the same spatial and temporal discretization as before. The temporal and spatial errors are shown in Figure 2.5 10-4 10-3 10-10 10-8 10-6 10-4 density velocity pressure 2 (a) 10-2 10-1 10-3 10-2 10-1 100 101 2 density velocity pressure (b) Figure 2.5: L1 error norm non-dimensionalized by the reference fields of the moving ball problem with time-dependent boundary temperature: (a) temporal and (b) spatial error, respectively. 2.6.2 Flow past circular cylinder Compressible flow over a stationary circular cylinder is used as another validation for the present method. Here, the adiabatic no-slip boundary condition is used, and the tests are carried out with various Reynolds numbers and resolutions to investigate the properties of the steady state and periodic vortex shedding solutions. A domain size dependence study was also carried out to ensure that the size of the computational domain was sufficiently large such that outflow boundary effects had minimal effect on the solution. Furthermore, the convergence of the solution as a function of the number of surface markers was investigated. The results are compared with experimental and numerical results reported in the literature. We investigate the drag, lift, and frequency of vortex shedding. In the present method, the boundary force is calculated directly from the solution without additional post-processing steps since they are directly related to the integral of f over the surface. The simulations are conducted in a rectangular domain with patch-based adaptive mesh refinements (AMR) for efficiency, implemented in SAMRAI [60]. The AMR implementation of the projection IB method will be discussed in Chapter 4. The mesh is discretized uniformly in the x and y directions with 33 the same grid size in each AMR level. The mesh is refined in the vicinity of the cylinder to a size that matches the spacing of the surface markers of the cylinder. The computational domain has streamwise length Lx = 30D, and the cylinder is placed at xc = 7.5D from the left boundary. The transverse length Ly = 12.5D is sufficient for our further studies as we tested. Table 2.2 shows the drag coefficients from the simulations with different numbers of surface markers Nb and Eulerian mesh resolutions ∆x (measured at the finest AMR level) at Re = 20 and 40. The simulations used 3 (for cases where Nb = 41,81) or 4 (for cases where Nb = 161) AMR levels, concentrated around the surface of the cylinder and the wake. It can be observed that the drag coefficient converges approximately at the rate of first order. The drag coefficients from our test cases with different resolutions are very close to each other, meaning that our method can capture the near-wall flow details with a relatively small number of grid points (for the Reynolds number investigated). The streamlines of the steady state flow for Re = 20 and Re = 40 are shown in Figure 2.6, where the expected recirculation zones can be observed behind the cylinder. Table 2.3 compare the drag coefficients Cd and the angle of the separation point θs with other experimental and numerical studies, and our results are in good agreement with previous literature. Table 2.2: Drag coefficients for Ma = 0.2 Re Nb ∆x/D Cd (% change) 20 41 0.156 2.141 - 20 81 0.078 2.122 -0.92 20 161 0.039 2.112 -0.47 40 41 0.156 1.594 - 40 81 0.078 1.580 -0.90 40 161 0.039 1.572 -0.49 34 (a) (b) Figure 2.6: Streamlines near the cylinder for (a): Re = 20 and (b) : Re = 40. Nb = 161. Table 2.3: Comparison of drag coefficients for Re = 20 and 40. Reference Ma θs Cd Re = 20 Dennis et al. [61] 0 43.7 ◦ 2.05 Taira et al. [28] 0 43.3 ◦ 2.06 Linnick et al. [62] 0 43.5 ◦ 2.06 Hartmann et al. [63] 0.1 - 2.03 Canuto et al. [64] 0.2 43.7 ◦ 2.08 Present 0.2 43◦ ∼ 44◦ 2.11 Re = 40 Dennis et al. [61] 0 53.8 ◦ 1.52 Taira et al. [28] 0 53.7 ◦ 1.54 Linnick et al. [62] 0 53.6 ◦ 1.54 Uddin et al. [65] 0.1 - 1.56 Hartmann et al. [63] 0.1 - 1.52 Canuto et al. [64] 0.2 53.7 ◦ 1.56 Present 0.2 53◦ ∼ 54◦ 1.57 Simulations at Re = 100 and 200 are also carried out with the same AMR grids and spatial discretization to investigate the unsteady vortex shedding. Simulations at both Reynolds numbers are conducted in three 35 different resolutions where Nb = 81, 161, and 322. Figure 2.7 shows the vorticity contour after the vortex street structure has developed and stabilized. Periodic vortex shedding can be observed behind the cylinder. The drag and lift coefficients are calculated and plotted as a function of non-dimensionalized time tU∞/d with phase shift in Figure 2.8. Both Cd and Cl oscillate in a manner that is similar to harmonic functions, which is consistent with the available literature. Table 2.4 reports the coefficients and Strouhal numbers for each resolution and compares the results with available literature. As seen, the results are in close agreement with previously published results. Figure 2.7: Vorticity contour for (top) Re = 100 and (bottom) Re = 200. Nb = 161. 36 0 5 10 15 20 1.39 1.4 1.41 1.42 1.43 0 5 10 15 20 1.36 1.39 1.42 1.45 1.48 0 5 10 15 20 -0.4 -0.2 0 0.2 0.4 0 5 10 15 20 -0.8 -0.4 0 0.4 0.8 Figure 2.8: Evolution of drag coefficients with time (top) after vortex shedding developed for (left) Re = 100 and (right) Re = 200, and the evolution of lift coefficients with time (bottom) for the corresponding Reynolds number. Solid: Nb = 322; Dash: Nb = 161; Dash-dot: Nb = 81. x-axis is non-dimensionalized and shifted for better comparison. Table 2.4: Comparison of drag coefficients for Re = 100 and 200. Reference Ma Cd Cl St Re = 100 Canuto el al. [64] 0.2 1.36±0.0093 ±0.332 0.165 Lai et al. [66] 0 1.447 ±0.330 0.165 Uddin et al. [65] 0.1 1.342 ±0.333 0.164 Karagiozis et al. [67] 0.25 1.336 ±0.319 0.168 Present (Nb = 81) 0.2 1.42±0.0102 ±0.364 0.165 Present (Nb = 161) 0.2 1.41±0.0099 ±0.358 0.166 Present (Nb = 322) 0.2 1.40±0.0101 ±0.354 0.166 Re = 200 Taira et al. [28] 0 1.36±0.043 ±0.69 0.197 Linnick et al. [62] 0 1.34±0.044 ±0.69 0.197 Present (Nb = 81) 0.2 1.43±0.0470 ±0.735 0.197 Present (Nb = 161) 0.2 1.42±0.0469 ±0.728 0.197 Present (Nb = 322) 0.2 1.41±0.0474 ±0.723 0.198 37 2.6.3 Acoustic scattering (a) (b) Figure 2.9: Contours of instantaneous non-dimensionalized pressure fluctuations for acoustic scattering probem. Here, we consider an inviscid problem of acoustic scattering used in [68] to demonstrate the capabilities of the present compressible IBM for inviscid flow. The Euler equation is solved in a rectangular computational domain [−10,10] × [−10,10] with outflow boundary conditions on the four faces. A cylinder of diameter D = 1 is placed at the center of the domain with slip wall boundary conditions. A source of acoustic energy is induced by forcing the energy equation; see [68]. The flow field is initialized with the undisturbed density ρ0 = 1.225, the undisturbed pressure p0 = 101,325, and zero velocity everywhere in the domain. The acoustic energy source is S = 1 4 ρ0c 3 0 D e − ln2 0.04 (x−4) 2+y 2 D2 sin ωc0 D t (2.100) with ω = 4, where c0 = p γ p0/ρ0 with γ = 1.4 the specific heat ratio. The small perturbation solution can be obtained analytically from the linearized Euler equation [68]. In the numerical test, the computational domain is discretized with 801 grid points in each direction, and the immersed cylinder is discretized by 106 uniformly spaced surface markers. The Roe-based approximate Riemann solver is used to compute the inviscid fluxes. The time step is calculated based on CFL = 0.5. 38 Figure 2.9 shows iso-contours and a surface plot of the instantaneous pressure fluctuations at the time tc0/D = 33.6. Note that the pressure fluctuations near the immersed boundary are smooth without developing irregularities despite the flow having no physical viscous diffusion. Figure 2.10 shows the pressure fluctuation as a function of time at the point (x = 2, and y = 0), midway between the center of the cylinder and the acoustic source, compared with the analytical solution. The simulation result shows good agreement with the analytical solution over many periods. 6 8 10 12 14 16 18 20 -5 0 5 10-3 Figure 2.10: Pressure fluctuations as a function of time at x = 2 and y = 0: simulations (solid) and analytical solution (dash). 2.6.4 Shock reflection from a wedge Although the intended target of our method is compressible flow at subsonic Mach numbers, since not much focus has been paid to supersonic conditions and the interaction of immersed boundary methods with shocks, we demonstrate that the method can simulate phenomena with higher levels of compressibility. We consider the problem of the shock reflection from a wedge, which is explained in detail in [69] and used as a validation in other immersed boundary methods, e.g., [65], using the Euler equations. The computational domain is [−4,4] × [−3,3]. The wedge has an angle of 25◦ and is modeled by an inclined line segment across the computational domain with slip wall boundary condition. It should be noted that the domain is essentially split into two parts, where the top outside the wedge is the physical region while the bottom one is the fictitious wedge region. The flow field is initialized uniformly with ρR = 1.225, uR = 0, pR = 101,325 and the left boundary condition is ρL = 2.6922, uL = 3.1527×102 , pL = 3.2475×105 . This setup of initial and boundary conditions initializes a normal shock of Mach number Ms = 1.7 traveling into the domain from the left boundary. The simulation is carried out with 800 and 600 grid points along the x and y directions, respectively, and the time step is calculated based on CFL = 0.5. Figure 2.11 shows the pressure 39 and Mach number contours when the shock has traveled approximately 3/4 of the length of the surface from the wedge’s apex. The characteristic features of the flow are observed, such as a reflected shock, Mach stem, and slip-line. (a) (b) Figure 2.11: Shock reflection from a wedge: (a) pressure and (b) Mach number; domain under the wedge is masked for clarity. 2.7 Conclusions A projection immersed boundary method for compressible viscous flow was developed to accommodate a complete set of kinematic and thermal boundary conditions. The boundary conditions are enforced with added singular source terms identified by the jump conditions. The singular source terms are algebraic variables in the system of differential-algebraic equations, and the boundary conditions are the algebraic constraints. The system of differential-algebraic equations is solved with half-explicit Runge-Kutta methods, and each stage is solved with a fractional-step approach. The method applies to stationary boundaries and boundaries with prescribed motion and can be applied to fluid-structure interaction problems when coupled with structure dynamic solvers. For the two-dimensional test problems, the method achieves second-order spatial accuracy in the L 1 norm when the advection and diffusions are discretized with second-order schemes. Second-order temporal accuracy can be achieved with the half-explicit RK method based on the three-stage SSP-RK method. Theoretically, the temporal order can be made higher with half-explicit Runge-Kutta methods with more stages. This was not investigated further, but the method was developed without a particular choice of discretization in mind. 40 The numerical results from the new immersed boundary method show good agreement with existing experimental, numerical, or analytical results. 41 Chapter 3 Tikhonov regularization of the projection immersed boundary method This chapter introduces a Tikhonov regularization method to accurately compute the boundary force for the projection IB method without ad hoc post-processing. The regularization method will be discussed with incompressible flow problems to avoid unnecessary complexity. However, it can be naturally applied to compressible flow problems, and the extension will be briefly discussed in the end. 3.1 Introduction The governing equations of the projection immersed boundary method for both incompressible [28] and compressible [70] flow are systems of differential-algebraic equations, where the boundary source terms (e.g., boundary force) are algebraic variables, and the prescribed values of the flow variables on the immersed boundary are the algebraic equations (constraints). The system is solvable with a fractional-step (projection) approach, and the boundary force as an algebraic variable has clear physical interpretations. Now, the projection method may not generally produce correct boundary force solutions even though the flow variables and integrals of boundary forces converge as the grid refines in many observations [28, 71]; this issue is well documented [72, 71]. The reason is that the boundary force is the solution of a discretized Fredholm integral equation of first kind [73, 74], where the discretized integral operator does not have a uniformly bounded inverse as the grid is refined. Therefore, the boundary force will not converge to the exact solution as the grid is refined. Solving the discrete integral equation is an inverse problem, and readers may refer to [75, 73] for more details about the difficulties. Moreover, the discretized integral operators generally have small singular values corresponding to highly oscillatory singular vectors, which implies that small high-frequency errors on the right-hand side of the integral equation will be amplified, resulting in 42 spurious oscillations in the boundary force solution. It is observed [72] that the modified Poisson equation is more ill-conditioned when the Lagrangian-to-Eulerian grid spacing ratio gets small. We also observe that the spurious oscillation becomes more violent as the spacing ratio decreases. However, a relatively small Lagrangian-to-Eulerian grid spacing ratio is generally required to prevent transpiration through the immersed boundary since the discrete delta function can only spread the boundary force to a few Eulerian grid points near each Lagrangian grid points, and in many FSI problems, the Lagrangian and Eulerian grid spacing might be determined by different factors that are hard to control. To obtain smoothing and converging boundary forces, Goza et al.[71] proposed a method that explicitly filters out the high-frequency components after the boundary force is solved. The filtering is combined with smeared discrete delta functions so that the actual high-frequency components decay fast, and thus, the filtering does not introduce too much error. The method leads to smooth and converging boundary force in several test cases. The smeared discrete delta function generally has significantly larger support. This method still requires the Lagrangian and Eulerian grids to have comparable spacing, and the discrete integral equation is still ill-conditioned. To address this problem, we propose a regularized projection method based on (generalized) Tikhonov regularization. Compared to the original projection method, the velocity constraint on the immersed boundary is replaced by an optimization problem that minimizes the fluid/immersed boundary velocity mismatches and the irregularity of the boundary forces. The regularization process modifies the discrete integral equation for the boundary force, as opposed to the method of Goza et al. [71], where the boundary force is filtered in a post-processing step. The method produces a smoother boundary force solution for a wide range of Lagrangian-to-Eulerian grid spacing ratios and shows convergence as the grid refines. This chapter will first recapitulate the original projection immersed boundary method in Section 3.2. The limitations of the projection method are discussed in Section 3.3, where the properties of the discretized integral operator for the boundary forces are studied with the singular value decomposition (SVD). The regularized immersed boundary will be introduced in Section 3.4. Some simulation results will be in Section 3.7. 3.2 The projection immersed boundary method We denote the computational domain of the fluid as Ω and the immersed boundary as Γ ⊂ Ω. The immersed boundary Γ is parametrized by s and possibly time-dependent. xbc(s,t) and x˙bc(s,t) denote the position 43 and velocity of the immersed boundary at point s and time t, respectively. The governing equations of the projection immersed boundary method for incompressible flow [28] are ∂u ∂t + (u·∇)u = − 1 ρ ∇p+ν∇ 2u+ 1 ρ Z Γ f(s,t)δ(x−xbc(s,t))ds in Ω, (3.1) ∇·u = 0 in Ω, (3.2) Z Ω u(x,t)δ(xbc(s,t)−x)dx = x˙bc(s,t) on Γ, (3.3) where u(x,t), p(x,t) are the velocity vector and pressure of the fluid, and f(s,t) is the boundary force on the immersed boundary. ρ and ν denote the density and kinematic viscosity, respectively. The convolution with Dirac delta function δ in Eq.(3.1) spreads the boundary force from the immersed boundary Γ to the computational domain Ω, and the convolution in Eq.(3.3) interpolates the fluid velocity from the computational domain Ω to the immersed boundary Γ. The system is spatially discretized into the following system du dt +N(u) = − 1 ρ Gp+ 1 ρ L f, (3.4) Du = 0, (3.5) Eu = x˙bc, (3.6) where u is the discrete velocity vector, p is the discrete pressure vector, and f is the discrete boundary force vector on the Lagrangian grid. ˙xbc discrete boundary velocity vector on the Lagrangian grid. N(u) ≈ (u · ∇)u − ν∇ 2u is the discretized convection-diffusion operator, G and D are the discrete gradient and divergence operators, respectively. The physical boundary terms arising from the differential operators are omitted here. In this chapter, the variables are discretized in staggered grids that are uniform in the vicinity of the immersed boundary, and all differential operators are discretized with a second-order centered difference scheme. L and E are the discretization of the spreading operator f 7→ R Γ f(s)δ(x−xbc(s))ds and interpolation operator u 7→ R Ω u(x)δ(xbc(s) − x)dx, respectively. Like the discussion in Chapter 2, E and L are constructed by the discrete convolution with discrete delta functions. It should be noted that the interpolation operator E is different when applied to velocity components in different directions because of the staggered grid where different velocity components are defined on their respective grid faces. The same 44 goes for the spreading operator L. It has been observed [28] that the operator L is the transpose and right scaling of E by Lagrangian element sizes, i.e., L = E TA, where A = 1 ∏ d k=1 ∆xd diag({∆sα}). Therefore, the implementation only requires the construction of one operator. For simplicity, we consider a temporal discretization that treats the advection and diffusion terms explicitly, although the Crank-Nicolson scheme is more widely used to discretize the diffusion terms to avoid small time steps. The resulting time-marching scheme is I G −L D 0 0 E 0 0 u n+1 pˆ ˆf = u n+1,⋆ 0 x˙ n+1 bc + bc1 bc2 0 , (3.7) where u n+1,⋆ = u n −∆tN(u n ) is the intermediate fluid velocity field. Here, ˆp = ∆t ρ p and ˆf = ∆t ρ f are scaled pressure and boundary force which help simplify the notation. bc1 is the terms from the pressure boundary condition and operator G, while bc2 arises from the velocity boundary condition and operator D. Higherorder methods for the Hessenberg index-2 DAE, e.g., the half-explicit Runga-Kutta method [51, 2], are also possible, which involves solving the system (3.7) at each stage. The system can be factorized by block-LU decomposition as I G −L D 0 0 E 0 0 = I 0 0 D −DG DL E −EG EL I G −L 0 I 0 0 0 I , (3.8) which allows us to solve directly with the fractional step method in three steps [28]: u n+1,⋆ = u n −∆tN(u n ), (solve for intermediate step) (3.9) Mq = b, (solve modified pressure Poisson for ˆp and ˆf) (3.10) u n+1 = u n+1,⋆ −(Gpˆ −L ˆf) +bc1, (projection step) (3.11) where M = −DG DL −EG EL , q = pˆ ˆf , b = bc2 x˙ n+1 bc − D E (u n+1,⋆ +bc1) ≡ b1 b2 . (3.12) 45 3.3 Challenges of the projection method The projection method has been used extensively to simulate the flow around non-deforming immersed boundaries, e.g., cylinders or flat plates [28]. The implicitly calculated boundary force has been used to probe the lift and drag on the bodies. However, the boundary force does not converge when refining the Lagrangian and Eulerian grids, as noted in [71], and in general, the boundary force solutions show spurious oscillation. This could be problematic in the fluid-structure interaction simulations, where accurate boundary forces might be needed. To investigate the reason for the lack of convergence and the source of the spurious oscillation, we rearrange Eq.(3.10) using the Schur complement and isolate the equation for the boundary force ˆf , given by EBL ˆf = b2 −EG(DG) −1 b1 = r, (3.13) where B = I −G(DG) −1D. (3.14) The equation is the discretized version of an integral equation of the first kind, as noted in [71], and the discretized integral operator EBL has small singular values corresponding to high-frequency components. For this reason, the solution of Eq.(3.13) will tend to exhibit oscillatory ˆf behavior even if r is smooth. We investigate the spectrum of EBL to identify the source of the spurious oscillations. Here, we use a test problem where the immersed boundary is a cylinder of diameter 1 centered at the origin, and the computational domain for the fluid is [−1,1] × [−1,1]. The Lagrangian and Eulerian grids are discretized with ∆x = 0.1 and various ∆s’s, and the operators E, B, and L are constructed in the exact manner described in Section 3.2. We calculate the singular value decomposition (SVD) of EBL = UΣV T , where U = [u1,··· ,uNL ], V = [v1,··· , vNL ] are unitary matrices, uk’s, vk’s are the left and right singular vectors, and Σ = diag(σ1,··· ,σNL ) is the diagonal matrix of the singular values σk ≥ 0. Note that zero singular values arise when the rank of E or L is smaller than the number of Lagrangian grid points; in this situation, EBL is singular. This happens when ∆s is finer than ∆x beyond a point so that there are more Lagrangian grid points than the Eulerian grid points affected by the immersed boundary. The SVD gives a solution (or the least-square solution 46 with minimum l 2 norm if EBL is singular) of EBL ˆf = r as ˆf = VΣ †U T r, where Σ † is the pseudoinverse of diagonal matrix Σ, formed by replacing every non-zero diagonal entry by its reciprocal. The solution is ˆf = ∑ σk̸=0 vk 1 σk (u T k r). (3.15) As shown in [71], ˆf does not converge even if the right hand side r converges as the grid is refined because of the small σk’s amplifies the error of r. Here, we plot the singular values σk against the total variation (which measures the degree of oscillation) of the right singular vector vk in Figure 3.1 for different ∆s’s. It can be seen that roughly higher frequency (larger total variation) right singular vectors correspond to smaller singular values, which is what we expect for an integral operator. This leads to oscillatory boundary force solution per Eq.(3.15) because of the large 1/σk associated with vk of large oscillation. Moreover, it can be seen from the plot that for smaller ∆s (i.e., finer Lagrangian grid), the ratio of σmax/σmin is significantly larger, leading to more ill-conditioned equations that are problematic when solving with iterative methods (also noted in [72]), and more oscillatory solutions. Figure 3.2 shows the condition numbers of EBL for various grid spacing ratios. The condition number increases significantly for smaller ∆s/∆x. 47 0 5 10 15 0 1 2 3 4 5 6 (a) 0 5 10 15 0 1 2 3 4 5 6 (b) 0 5 10 15 0 1 2 3 4 5 6 (c) 0 5 10 15 0 1 2 3 4 5 6 (d) Figure 3.1: Singular values of EBL against the total variation of the right singular vectors. (a): ∆s = 0.0698; (b): ∆s = 0.0980; (c): ∆s = 0.1253;(d): ∆s = 0.1564. Figure 3.2: Condition number of EBL for various ∆s/∆x. The region where EBL is singular is gray-shaded. 48 We demonstrate this issue with simulations of the uniform flow of velocity u∞ passing a cylinder of diameter D at Re = u∞D/ν = 40. The domain is sufficiently large such that the exterior boundary does not affect the overall solution around and behind the cylinder. Figure 3.3 shows the force on the cylinder surface. It can be seen that the boundary force is significantly more oscillatory when the Lagrangian markers are slightly more densely distributed, which is in agreement with the previous SVD analysis. In Figure 3.3c, the magnitude of spurious oscillations completely overwhelms the physical surface forces. -2 0 2 -0.2 0 0.2 0.4 0.6 0.8 1 (a) -2 0 2 -0.2 0 0.2 0.4 0.6 0.8 1 (b) -2 0 2 -0.2 0 0.2 0.4 0.6 0.8 1 (c) Figure 3.3: The boundary force on the cylinder with different ratios of Lagrangian-Eulerian grid spacing. (a): ∆s/∆x ≈ 1.6; (b): ∆s/∆x ≈ 1.3; (c): ∆s/∆x ≈ 1. 3.4 A regularized projection immersed boundary method The goal is to modify the projection immersed boundary method as described in Section 3.2, to suppress spurious oscillations in the boundary force and better approximate the actual boundary force distribution. As has been shown before, the issue is due to the nature of the discrete integral equation (3.13), whose small singular values cluster near zero, and a typical strategy to alleviate this issue is the (generalized) Tikhonov regularization [76, 75, 73]. 49 3.4.1 Weighted least square formulation of the projection method We cast the formulation of the projection immersed boundary method, Eq.(3.7), into a minimization problem as the starting point of our regularized projection immersed boundary method. The minimization problem is minimizing Π = D E u n+1 − bc2 x n+1 bc 2 P , (3.16) subject to u n+1 + Gpˆ −L ˆf = u n+1,⋆ +bc1, (3.17) where P is symmetric positive definite and ∥x∥P = √ x TPx is the norm induced by P. The functional Π to be minimized is the square of the residual norm of the divergence-free constraint (i.e., the second row of (3.7)) and the immersed boundary velocity constraint (i.e., the third row of (3.7)). The constraint (3.17) corresponds to the momentum equation, i.e., the first row of (3.7). The functional Π can be re-written in terms of ˆf and ˆp instead of u n+1 by (3.7): Π = ∥Mq−b∥ 2 P = q TMTPMq−b TPMq−q TMTPb+b TPb, (3.18) where M, q, b are defined in Eq.(3.12). Here in Eq.(3.18) the functional Π is a quadratic form with respect to q, and the minimization problem gives a solution of q = [pˆ, ˆf ] by simply letting ∂Π/∂q = 0, i.e., MTPMq = MTPb, (3.19) in case that MTPM is non-singular (or equivalently, the square matrix M is non-singular since P is nonsingular.) The least square formulation Eq.(3.19) is equivalent to Eq.(3.10) with a non-singular M. 5 3.4.2 The Tikhonov regularized formulation We now move to the regularized projection method. The handling of the ill-conditioned discrete integral equation is achieved with the technique of (generalized) Tikhonov regularization, where an additional penalty term is added to the functional Π in Eq.(3.16). This leads to a new minimization problem minimizing Πλ = D E u n+1 − bc2 x n+1 bc 2 P +λ ˆf 2 R , (3.20) subject to u n+1 + Gpˆ −L ˆf = u n+1,⋆ +bc1. (3.21) Equivalently, writing Πλ in terms of ˆp and ˆf to reach an unconstrained minimization problem: minimizing Πλ = ∥Mq−b∥ 2 P +λ ˆf 2 R , (3.22) for arguments p and ˆf . Herein, R is a symmetric positive semi-definite matrix that induces the smoothing seminorm ˆf R = p ˆf TR ˆf that measures the irregularity of the boundary force. It should be noted that in two-dimensional problems, ˆf = [ ˆfx, ˆfy] T and R should be a block diagonal matrix which might be denoted as R = diag(R ′ ,R ′ ), with two identical symmetric positive semi-definite matrices R ′ . The smoothing seminorm ˆf 2 R = ˆf TR ˆf = ˆfxR ′ ˆfx + ˆfyR ′ ˆfy = ˆfx 2 R′ + ˆfy 2 R′ , similarly for three-dimensional problems. For simplicity, we do not distinguish between R and R ′ later. Herein, λ ≥ 0 is the regularization parameter. It can be observed heuristically that the regularized solution is a trade-off between small residual norm ∥Mq−b∥ 2 P of the constraints and small smoothing seminorm ˆf 2 R of the boundary force. The minimization problem (3.22) gives a solution by letting ∂Πλ /∂q = 0, noting that Πλ is a quadratic form with respect to q: MTPM +λR˜ q = MTPb, (3.23) provided that left hand side is non-singular, i.e., N (M)∩ N (R˜) = {0} (see Ref. [73], Chapter 8). Here, R˜ = diag(0,R) and Rq˜ = [0;R ˆf ]. Compared to Eq.(3.19), the extra term λRq˜ accounts for the regularization. 3.4.3 Constructing the smoothing seminorm with the graph Laplacian The smoothing seminorm ∥ f ∥R = p f TR f measures the degree of oscillation of the boundary force vector f = { fα} with respect to the location of markers {xbc(sα)}. The matrix R should have larger eigenvalues λk corresponding to more oscillatory eigenvectors vk, and a zero eigenvalue for the constant eigenvector, since ∥ f ∥R = p ∑k λk⟨f, vk⟩ 2 . The Sobolev seminorms, e.g., ∥∇x f ∥2 , are common choices for measuring the degree of oscillation. Following this idea, R can be chosen from the family of discrete differential operators, e.g., the discrete Laplacian operator. Noting that the immersed surface is simply discretized into the Lagrangian grid, composed by the cluster of Lagrangian markers {xbc(sα)}, without actually parametrizing the surface, we resort to the graph theory and use the graph Laplacian [77] as the operator R. For a Lagrangian grid with NL markers located at xbc(sα) = [xbc(sα), ybc(sα)], α = 1,··· ,NL, we define an un-directed weighted graph over the set of markers with the weight matrix W, where the entry Wαβ depends on the distance between xbc(sα) and xbc(sβ ), given by Wαβ = 1 ∆x∆y w xbc(sα)−xbc(sβ ) ∆x w ybc(sα)−ybc(sβ ) ∆y , (3.24) and w(r) is the Gaussian function with cutoff w(r) = exp −r 2 r < rc, 0 r ≥ rc, (3.25) where rc is the non-dimensional cutoff length. Here, the cut-off length is set to 4, where exp −r 2 c = 1.1×10−7 is already sufficiently small. W is real and symmetric by construction. The Lagrangian markers α and β are considered connected in the graph constructed if Wαβ ̸= 0, i.e., xbc(sα)−xbc(sβ ) < rc∆x and ybc(sα)−ybc(sβ ) < rc∆y. The graph Laplacian R is constructed by subtracting the sum of all entries in each row of W from the diagonal of W, i.e., Rαβ = δαβ NL ∑ γ=1 Wαγ −Wαβ . (3.26) The graph Laplacian R is a symmetric positive semi-definite matrix. The constant [1,1,··· ,1] lies in the nullspace of R since it is easily verifiable that ∑β Wαβ = 0. For three-dimensional immersed boundary problems, the length scaling factor in Eq.(3.24) might be set to 1/(∆x∆y∆z) 2/3 instead of 1/(∆x∆y). To verify some properties of the graph Laplacian, we construct the graph Laplacian for the cylinder described in Section 3.3 with ∆s/∆x ≈ 1.3 and investigate its eigenvectors. Figure 3.4 shows the eigenvectors (normalized by their ℓ ∞ norms and plotted against the angular coordinate with respect to the x axis) of the graph Laplacian corresponding to the 10th and 40th eigenvalues (in ascending order, non-dimensionalized by multiplying with ∆x∆y). The high frequency 40th component has a larger eigenvalue than the low frequency 10th component. For a more rigorous demonstration, we plot the eigenvalues against the total variations of their corresponding eigenvectors. The total variation serves as a measurement of component frequency. It can be seen that the eigenvalues of the graph Laplacian R are proportional to the square of the total variations of the corresponding eigenvectors, which means that ∥·∥ 2 R would be an appropriate seminorm to measure regularity. -2 0 2 -1.5 -1 -0.5 0 0.5 1 1.5 (a) -2 0 2 -1.5 -1 -0.5 0 0.5 1 1.5 (b) Figure 3.4: The eigenvectors of the graph Laplacian for the cylinder problem described in Section 3.3 with ∆s/∆x ≈ 1.3, and its corresponding eigenvalues, non-dimensionalized by ∆x∆y (a) and (b) : the 10th and 40th eigenvectors of the graph Laplacian. 53 101 102 10-3 10-2 10-1 100 Figure 3.5: Total variation of eigenvectors and their corresponding eigenvalues for the Graph Laplacian operator constructed for the cylinder grid. 3.4.4 The regularized immersed boundary equation One can observe that the only difference between the original projection method and the regularized projection method is that the modified Poisson equation (3.10), is replaced by the regularized version, Eq.(3.23). The steps of solving the intermediate velocity, Eq.(3.9), and velocity projection, Eq.(3.11), remain unchanged. The regularization process does not require calculating the SVD, which is impractical owing to the large computational cost. We now focus on the regularized Eq.(3.23). The matrix R that induces the smoothing seminorm has been constructed with the Graph Laplacian, and the remaining problem is to choose the matrix P that induces the residual norm and the regularization parameter λ. Here, we introduce two typical methods of choosing P, which lead to the Tikhonov-type and Lavrentiev-type regularization. 3.4.4.1 Method 1: Tikhonov-type regularization The Tikhonov regularization is the most common way to regularize the discrete integral equation. In the Tikhonov regularization, the residual norm is the ℓ 2 norm (i.e., the residual norm is induced by the identity matrix), and the smoothing seminorm can be the ℓ 2 norm, or more generically, any seminorm. The Tikhonov 54 regularization does not require the matrix M to be symmetric. Following this practice, we set P that induces the residual norm to be a diagonal matrix P = γI 0 0 I , (3.27) where the first diagonal block corresponds to the divergence-free constraint and has the same size, and the second diagonal block corresponds to the immersed boundary velocity constraint, noting Eq.(3.20). γ is a preset constant that will be determined later. The Tikhonov functional under the choice of P in (3.27) is Πλ = γ Dun+1 −bc2 2 2 + Eun+1 −x n+1 bc 2 2 +λ ˆf 2 R . (3.28) The coefficient γ and regularization parameter λ are determined by a simple dimensional analysis. Herein, we denote the characteristic flow velocity and length as U¯ and L¯, and the numbers of Eulerian and Lagrangian grid points as NE and NL. The Eulerian grid spacing ∆x ∼ L¯/N 1/d E and the Lagrangian grid spacing ∆s ∼ L¯/N 1/(d−1) L , and NE/NL ∼ L∆s d−1/∆x d . Noting that Du ∼ U¯ /L¯, Eu ∼ U¯ and ˆf ∼ U¯ L¯, the magnitude of the three terms Dun+1 −bc2 2 2 ∼ U¯ 2/L¯ 2NE, Eun+1 −x n+1 bc 2 2 ∼ U¯ 2NL and ˆf 2 R ∼ U¯ 2NL. Based on the magnitude, we choose γ = CL¯ 2 ∆x d ∆s d−1 1 L¯ , where ∆x and ∆s are chosen as the average Eulerian and Lagrangian grid spacing, and L¯ is chosen to be p LxLy of the Cartesian domain. The non-dimensional C is chosen to be sufficiently large to enforce the divergence-free constraint. The regularization parameter λ is non-dimensional and from the preliminary dimensional analysis λ ∼ 1. The optimal λ might be problem-specific and has to be determined case-by-case. The discrepancy principle [73] is the simplest manner to determine λ. However, it relies on a good estimate of the error norm of the right hand side. In Section 3.7, we will select the optimal λ using the L-curve for the steady problem of the flow passing a cylinder and the impulsively started rotating cylinder problem. The basic principle of the L-curve method is to choose several λ, observe the error and smoothness of the results, and determine the optimal λ accordingly. The results in Section 3.7 indicate that the regularization parameter λ can be selected from a wide range that spans several orders of magnitude without significantly compromising accuracy. Alternatively, we might resort to an iterative regularization method [73] to avoid selecting λ at all; this is beyond the scope of this chapter. 55 3.4.4.2 Method 2: Lavrentiev-type regularization The Tikhonov-type regularization requires calculating matrix-matrix product MTPM. This could be computationally inefficient. Also, this might lead to a larger condition number for small regularization parameter λ. Herein, we introduce another choice of P that does not require the calculation of MTPM. This regularization uses the fact that M could be symmetric with proper right scaling and is similar to the Lavrentiev regularization [78]. Recall that G = −D T , and L = E TA where A is the diagonal scaling matrix of Lagrangian and Eulerian grid spacing ratio, we rewrite M in Eq.(3.10) as M = D E D T E T I 0 0 A , (3.29) and choose P to be P = D E D T E T −1 , (3.30) which is symmetric positive definite. With some manipulation of Eq.(3.23) by substituting P from Eq.(3.30), we obtain a regularized equation −DG DL −EG EL +λ 0 0 0 A −1R pˆ ˆf = b1 b2 , (3.31) where A −1 = ∆x∆ydiag(1/∆s1,···) and ∆sα is the size of the Lagrangian marker α. Compared to the modified Poisson equation (3.10) of the original projection method, A −1R is added to account for the regularization, and the added computational cost is to calculate A −1R ˆf if the modified Poisson equation is solved iteratively. Note that A −1R ˆf can be evaluated on the Lagrangian grid without resorting to the Eulerian grid. For this reason, the Lavrentiev-type regularization is generally more practical in terms of computational cost than the Tikhonov-type regularization and might be preferred in large-scale three-dimensional problems. The preliminary dimensional analysis similar to the Tikhonov regularization in Section 3.4.4.1 also indicates that the regularization parameter λ ∼ 1 in the Lavrentiev-type regularization, and we will demonstrate the process of using L-curve to find the optimal λ for the test problems in Section 3.7. 56 3.5 Regularization of compressible IB method The projection IB method for compressible flow discussed in Chapter 2 requires solving several systems algebraically similar to Eq.(3.7), e.g., Eqs.(2.61), (2.64), (2.67), and (2.71). The dimensions of these systems are significantly smaller compared to Eq.(3.7) of incompressible flow problems, without the pressure and divergent free constraints. The boundary source terms of these systems are solutions of discrete integral equations (e.g., Eqs.(2.62), (2.65), (2.69), and (2.72)) due to the nature of operators E and L, and spurious oscillation will be observed in the boundary source solutions similar to what we observe in Section 3.3. To alleviate these problems, we apply Tikhonov regularization to the projection IB method for compressible flow. With the regularization, systems (2.61), (2.64), (2.67) and (2.71) are replaced by minimization problems with smoothing seminorm as penalty terms. The system (2.61) is replaced by the minimization problem minimizing Π = ∥Eu˜ ν −x˙bc(t˜ ν )∥ 2 P +λ ˆf 2 R′ , (3.32) subject to ˜u ν −diag(1/ρ˜ ν )L ˆf = u˜ ν,⋆ , (3.33) where ˆf = γν∆t ˜f ν−1 . Equivalently, by the expression of ˜u ν , the function Π is written in terms of ˆf only as Π = E diag(1/ρ˜ ν )L ˆf −b 2 P +λ ˆf 2 R′ , (3.34) with b = x˙bc(t˜ ν )−Eu˜ ν,⋆. For easy implementation we choose the matrix P = E 1 ρ˜ ν E T similar to the Lavrentiev type regularization, noting that L = E TA and A = diag(∆s1,···)/∆x∆y is diagonal scaling matrix of grid spacing. For dimensional consistency, the graph Laplacian R has to be scaled by the reciprocal of the density of magnitude comparable to ρ˜ ν near the IB, and the smoothing seminorm is based on the scaled graph Laplacian R. Herein, we choose the R ′ = DρRDρ where the graph Laplacian is left and right scaled by diagonal matrix Dρ = diag({1/ √ Eρ˜ ν }), i.e., interpolating ρ˜ ν to Lagrangian grid points with E and take square root and reciprocal element-wise. This guarantees the consistency of dimension and the positive semi-definiteness of R ′ . With the chosen P and R ′ , the functional becomes Π = ˆf TAE diag(1/ρ˜ ν )L ˆf −2 ˆf TAb+λ ˆf DρRDρ ˆf , (3.35) 57 and the minimization of the quadratic form of ˆf leads to E diag(1/ρ˜ ν )L ˆf +λA −1DρRDρ ˆf = b, (3.36) where the diagonal grid space scaling A and the density scaling Dρ are given above. λA −1DρRDρ is added to the right hand side compared to the original Eq.(2.62) for the boundary force. Regarding implementation, the regularized equation for ˆf replaces Eq.(3.13), and the projection step is not altered. The same modification of regularization goes to the system (2.75) with the same added regularization term. The regularization for system (2.64) is similar, by replacing ˆf with ˆg = γν∆tg˜ ν−1 and the right hand side b = RTbc(t˜ ν ) γ−1 − E(e˜ ν,⋆ − 1 2 u˜ ν · u˜ ν ). The system (2.71) is replaced similarly by the minimization problem minimizing Π = d ∑ i=1 NiEq˜i −qn,bc 2 P +λ θ˜ 2 R′ , (3.37) subject to ˜qi −diag(κ)LNiθ˜ = q ⋆ i , 1 ≤ i ≤ d, (3.38) where the stage superscript ν −1 for tilde variables is omitted for clarity. By substituting the ˜q with θ˜ the function Π is becomes Π = d ∑ i=1 NiE diag(κ)LNi ! θ˜ −b 2 P +λ θ˜ 2 R′ , (3.39) where b = qn,bc − ∑ d i=1NiEq⋆ i . With the same Lavrentiev-type regularization approach, we simply choose P = ∑ d i=1NiE diag(κ)E TNi , noting that L = E TA and the multiplication of diagonal matrices A and Ni commutes. Also, for dimensional consistency, we choose R ′ = DκRDκ where Dκ = diag({ √ Eκ˜}) is the diagonal of the square root of the heat conduction coefficients interpolated to the Lagrangian grid points. Minimizing the quadratic form Π by setting ∂Π/∂θ˜ leads to the linear equation for θ˜ d ∑ i=1 NiE diag(κ)LNi ! θ˜ +λA −1DκRDκ θ˜ = b. (3.40) 58 The regularization of system (2.67) is done similarly by replacing it with the minimization problem minimizing Π = Pτ d ∑ j=1 NjEσ˜ j −στ,bc 2 P +λ ˜ζ 2 R′ , (3.41) subject to σ˜i j +diag(µ)L(Ni ˜ζj +Nj ˜ζi) = σ ⋆ i j, 1 ≤ i, j ≤ d, (3.42) with superscript ν −1 omitted for clarity. With some commutation error by taking Ni and N out, Π is again written in terms of ˜ζ as Π = d ∑ j=1 NjE diag(µ)LNj ! ˜ζ −b 2 P +λ ˜ζ 2 R′ , (3.43) where b = Pτ ∑ d j=1NjEσ ⋆ j −στ,bc. Similarly, by choosing P = ∑ d i=1NiE diag(µ)E TNi , Dµ = diag({ √ Eµ˜ }) and R ′ = DµRDµ and taking ∂Π/∂ ˜ζ = 0, we reach the equation for ˜ζ : d ∑ i=1 NiE diag(µ)LNi ! ˜ζ +λA −1DµRDµ θ˜ = b. (3.44) Herein, λA −1DµRDµ is the extra penalty term for regularization, compared to the original Eq.(2.69). The above equations can be simplified for constant κ and µ we use by taking out the constant diagonal scaling of κ and µ. In all these cases, the regularization parameters λ’s can be chosen similarly by inspecting the L-curve. We need to solve multiple such minimization problems in each step or R-K stage for multiple singular source terms that appear in the compressible version of IB methods, and we choose the same λ since the choice of λ should depend only on the grid layout. 3.6 Preliminary investigation of the oscillatory patterns near boundary In the projection method, the boundary force is the solution of the linear equation K ˆf = r, (3.45) where K = EBL is the discrete integral operator of a continuous version Z Γ K(ξ −ξ ′ ) ˆf(ξ ′ )dξ ′ = r(ξ ). (3.46) 59 It is also observed that the boundary force, calculated from the projection method, shows spurious oscillatory patterns spatially within several Lagrangian grid points touching the edge ∂Γ of the immersed boundary Γ if Γ is open. This is due to the nature of integral operators that amplify high-frequency components, as investigated before, but also the Gibbs phenomenon since by the construction of the discrete convolution operator, ˆf(ξ ′ ) is zero outside Γ 1 . This is demonstrated with a flat boundary with length 1, Eulerian grid spacing ∆x = 0.02, and Lagrangian grid spacing ∆s = 0.0217 when constructing E and L, and K = EL. Figure 3.6a demonstrates the results when a smooth ˆf is prescribed and r is calculated via r = K ˆf , with drop of amplitude near ξ = 0 and ξ = 1 observed in r due to the zero padding of ˆf . Conversely, Figure 3.6b demonstrates the results when a smooth r is prescribed and ˆf is calculated by solving K ˆf = r, the oscillatory patterns near ξ = 0 and ξ = 1 are observed in ˆf . The patterns of ˆf disappear when the prescribed r decays to zero at ξ = 0 and ξ = 1, as demonstrated in Figure 3.6c. 1Similar phenomenon, arising from similar reasons, is observed in digital image processing and the oscillation is generally called ringing artifacts by the community. 60 r = K ˆf (a) ˆf = K −1 r (b) ˆf = K −1 r (c) Figure 3.6: The oscillatory patterns at boundaries when solving the linear equation for the boundary force, demonstrated with a flat immersed boundary with 0 ≤ ξ ≤ 1. To mitigate the issue, we modify the operators E and L such that ˆf will be extrapolated at the boundary ∂Γ instead of treated as zero immediately outside Γ. Herein, the Lagrangian markers are categorized into 61 interior markers set Ii , and ghost markers set Ig. The ghost markers are near ∂Γ, and their boundary force values are interpolated from interior markers with the coefficient matrixC ∈ R |Ig|×|Ii | , where f g m = ∑ |Ii | n=1Cmn fn for force f g m at ghost marker m, where fn is the force at interior marker n. The new spreading operator L g takes values of Ii , interpolates to Ii ∪Ig, and spreading to the Eulerian grids, i.e., L g φ| i, j ≡ |Ii | ∑ n=1 φnδ(xi, j −xbc(sn))∆sn + |Ig| ∑ m=1 |Ii | ∑ n=1 Cmnφn ! δ(xi, j −xbc(s g m))∆sm, (3.47) for any {φn} of interior markers, xbc(sn) and xbc(s g m) are position of interior marker n and ghost marker m, respectively. The new interpolation operator E g only consider the interior markers, i.e., E gΦ| n = ∑ i, j Φi, jδ(xbc(sn)−xi, j)∆x∆y, (3.48) for any {Φi, j} of Eulerian grid. Note that L g is no longer the transpose of E g with scaling. Therefore, only the Tikhonov type regularization in Section 3.4.4.1 can be used alongside this extrapolation treatment, strictly speaking. We revisit the flat plate test and replace E and L with E g and L g . The first two markers at ξ = 0 and the last two markers at ξ = 1 are chosen as the ghost markers, and the values of ghost markers are equal to those of interior markers immediately adjacent to the ghost markers. We choose the same r as used in Figure 3.6b and solve K g ˆf = r where K g = E gL g . Figure 3.7 shows the solution of ˆf , with the boundary oscillation pattern eliminated compared to Figure 3.6b. ˆf = K −1 r Figure 3.7: The solution of E gL g ˆf = r with operators E g and L g constructed when markers at boundaries are categorized as ghost markers. The right hand side r is the same as in Figure 3.6b. 62 3.7 Results Several examples are used to showcase the performance of the regularization and the accuracy of the solutions. 3.7.1 Flow passing cylinder with regularization We revisit the problem of the uniform flow passing a cylinder of diameter D with the regularized projection immersed boundary method. We consider the Reynolds numbers Re = u∞D/ν = 20, 40 for a steady-state solution and 200 for an unsteady periodic vortex shedding solution, where u∞ is the uniform inflow velocity and ν is the kinematic viscosity. The dimensions of the computational domain are [−30D,45D] in the streamwise direction and [−37.5D,37.5D] in the lateral direction. The center of the cylinder is located at (−10D,0). The Dirichlet condition (u, v) = (u∞,0) is prescribed on the rectangular domain’s left, top, and bottom boundaries. The convection boundary condition ∂u ∂t +u∞ ∂u ∂ x = 0 alongside with pressure condition p = 0 is prescribed on the right edge of the computational domain to allow the wake to propagate outward freely. The computational domain is sufficiently large so that the effect of boundaries does not significantly affect the results around the cylinder and the wake. The computational domain is discretized non-uniformly except in the vicinity of the cylinder, where the grid is uniform with spacing ∆x/D = 0.025. The cylinder is discretized with an almost uniformly spaced Lagrangian grid, with various Lagrangian-to-Eulerian grid spacing ratios ∆s/∆x. We now determine the appropriate choice of the regularization parameter λ. The suitable regularization parameter λ is problem-specific. It can even be time-dependent for unsteady problems, depending on the transient local flow structure near the immersed boundary. Here, we first perform the simulations at Re = 40 with three different Lagrangian-to-Eulerian grid spacing ratios, ∆s/∆x = 1.3, 1.0, and 0.7. Both Tikhonov-type and Lavrentiev-type regularizations will be used with various regularization parameters λ, and the suitable λ will be determined with the help of the L-curve, where the smoothing seminorm ˆf 2 R is plotted against the residual norm ∥Eu∥ 2 2 for various λ’s. Figure 3.8 and Figure 3.9 show the L-curves for the Tikhonov-type and Lavrentiev-type regularization, respectively. It can be seen that both types of regularizations achieve similar levels of accuracy and smoothness of the boundary force, although the Lavrentiev-type regularization requires larger λ. The irregularity of the boundary force decreases dramatically with the increase of λ when λ is not too large. For ∆s/∆x = 1.3, 1.0, and 0.7, the L-curves plateau approximately after 63 λ > 0.01 for the Tikhonov-type regularization and λ > 1 for the Lavrentiev-type regularization, meaning that the boundary force could be over-smoothed for larger λ values. 10 -5 10 -4 10 -3 10 -3 10 -2 10 -1 0.0001 0.0005 0.001 0.005 0.01 0.050.1 0.5 (a) 10 -4 10 -3 10 -3 10 -2 10 -1 0.0001 0.0005 0.001 0.005 0.01 0.050.1 0.5 (b) 10 -4 10 -3 10 -3 10 -2 10 -1 0.0001 0.0005 0.001 0.005 0.01 0.05 0.1 (c) Figure 3.8: The residual norm ∥Eu∥ 2 2 vs. smoothing seminorm ˆf 2 R at different regularization parameters λ for the Tikhonov-type regularization, Re = 40. (a): ∆s/∆x = 1.3; (b): ∆s/∆x = 1.0; (c): ∆s/∆x = 0.7. 10 -5 10 -4 10 -3 10 -3 10 -2 10 -1 0.01 0.05 0.1 0.51 5 10 (a) 10 -4 10 -3 10 -3 10 -2 10 -1 0.01 0.05 0.1 0.51 5 10 (b) 10 -4 10 -3 10 -3 10 -2 10 -1 0.01 0.05 0.1 0.5 1 5 (c) Figure 3.9: The residual norm ∥Eu∥ 2 2 vs. smoothing seminorm ˆf 2 R at different regularization parameters λ for the Lavrentiev-type regularization, Re = 40. (a): ∆s/∆x = 1.3; (b): ∆s/∆x = 1.0; (c): ∆s/∆x = 0.7. Figure 3.10 and Figure 3.11 show the boundary force along the x-direction for unregularized, Tikhonovtype regularized and Lavrentiev-type regularized simulations at ∆s/∆x = 1.3 and 1.0. It can be seen visually that the regularized projection method produces significantly smoother results compared to the original projection method for both Lagrangian-to-Eulerian grid spacing ratios. 64 -2 0 2 -0.2 0 0.2 0.4 0.6 0.8 1 (a) -2 0 2 -0.2 0 0.2 0.4 0.6 0.8 1 (b) -2 0 2 -0.2 0 0.2 0.4 0.6 0.8 1 (c) Figure 3.10: The x-direction boundary force on the cylinder for ∆s/∆x = 1.3. (a): without regularization; (b): Tikhonov-type regularization with λ = 0.005; (c): Lavrentiev-type regularization with λ = 0.5 -2 0 2 -0.2 0 0.2 0.4 0.6 0.8 1 (a) -2 0 2 -0.2 0 0.2 0.4 0.6 0.8 1 (b) -2 0 2 -0.2 0 0.2 0.4 0.6 0.8 1 (c) Figure 3.11: The x-direction boundary force on the cylinder for ∆s/∆x = 1.0. (a): without regularization; (b): Tikhonov-type regularization with λ = 0.005; (c): Lavrentiev-type regularization with λ = 0.5 We now investigate flow properties to demonstrate the accuracy of the regularized projection method. Here, the grid ratio is ∆s/∆x = 1.3. We first consider the simulations at Re = 20 and 40 to produce steady results and study the drag coefficients Cd = 2Fx/(ρu 2 ∞D), separation angle θs , locations of wake center a, b and the length of the recirculation zone l; See Figure 3.12 for the definitions. The results are reported in Table 3.1 and compared with previous literature. The drag coefficients agree with previous studies for both types of regularization and various parameters λ in a wide range from the under-smoothing to the oversmoothing region. The characteristic lengths a/D, b/D, and l/D of the recirculation wake also agree with previous results, where the error is less than one grid width ∆x/D = 0.025. Next, we consider the problem at Re = 200 to produce the unsteady vortex shedding and study the lift coefficient CL = 2Fy/(ρu 2 ∞D), drag coefficient CD = 2Fx/(ρu 2 ∞D) and Strouhal number St = f D/u∞. The results are reported in Table 3.2. The lift/drag coefficients and the vortex shedding frequency agree well with previous numerical simulations. 65 Figure 3.12: Streamlines of the flow passing a cylinder at Re = 40. The definitions of characteristic lengths a, b, and l of the recirculation wake are annotated. The end of the recirculation zone and the centers of the vortices can trivially be probed from the simulation results because both x and y velocities change sign there. 66 Table 3.1: Drag coefficients, separation angle and wake characteristic lengths for Re = 20 and 40, with comparison to previous literature. Regularization Type Regularization Parameter θs CD a/D b/D l/D Re = 20 Tikhonov-type 0.001 43◦ 2.08 0.37 0.45 0.95 0.01 43◦ 2.08 0.37 0.45 0.97 0.1 43◦ 2.08 0.37 0.45 0.97 Lavrentiev-type 0.05 43◦ 2.08 0.37 0.45 0.97 0.5 43◦ 2.08 0.37 0.45 0.97 5 43◦ 2.08 0.37 0.45 0.97 Dennis et al. [61] 43.7 ◦ 2.05 - - 0.94 Taira et al. [28] 43.3 ◦ 2.06 0.37 0.43 0.94 Linnick et al. [62] 43.5 ◦ 2.06 0.36 0.43 0.93 Re = 40 Tikhonov-type 0.001 54◦ 1.55 0.76 0.60 2.37 0.01 54◦ 1.55 0.76 0.60 2.37 0.1 54◦ 1.55 0.76 0.60 2.37 Lavrentiev-type 0.05 54◦ 1.55 0.76 0.60 2.37 0.5 54◦ 1.55 0.75 0.60 2.37 5 54◦ 1.55 0.76 0.60 2.37 Dennis et al. [61] 53.8 ◦ 1.52 - - 2.35 Taira et al. [28] 53.7 ◦ 1.54 0.73 0.60 2.30 Linnick et al. [62] 53.6 ◦ 1.54 0.72 0.60 2.28 67 Table 3.2: Drag/lift coefficients and Strouhal number of vortex shedding for Re = 200, with comparison to previous literature. Regularization Type Regularization Parameter St CD CL Re = 200 Tikhonov-type 0.001 0.196 1.36±0.048 ±0.70 0.01 0.196 1.36±0.047 ±0.70 0.1 0.196 1.36±0.048 ±0.70 Lavrentiev-type 0.05 0.196 1.36±0.048 ±0.70 0.5 0.196 1.36±0.047 ±0.70 5 0.196 1.37±0.048 ±0.70 Taira et al. [28] 0.196 1.35±0.048 ±0.68 Linnick et al. [62] 0.197 1.34±0.044 ±0.69 3.7.2 Impulsively started rotating cylinder Here, we consider the two-dimensional problem of a cylinder of radius r = 1 submerged in the fluid. The cylinder rotates at the constant angular velocity ω = 1, and the fluid is quiescent initially. The Reynolds number Re = ωr 2 ν = 10. Goza et al. [71] demonstrate the accuracy of their filtering technique with the problem. It should be noted that the flow will develop both inside and outside the cylinder. The problem has an exact solution provided that the fluid domain is infinite, where the velocity uex = uex(r,t)eˆθ and the azimuthal velocity is uex(r,t) = r +2∑ ∞ n=1 J1( √ λnr) √ λnJ0( √ λn) exp − λnt Re , r ≤ 1 L −1 h K1(r √ sRe) sK1( √ sRe) i , r > 1 (3.49) where Re is the Reynolds number, Jp is the p-th Bessel function of the first kind, √ λn is the n-th root of J1, K1 is the 1-st modified Bessel function of the second kind. L −1 [·] is the inverse Laplace transform with respect to variable s to the time t. The tangential boundary source fex(t) is given by fex(t) = 1 Re ∂ ∂ r uex r r=1− − 1 Re ∂ ∂ r uex r r=1+ (3.50) 68 It should be noted that although the exact tangential boundary force is a constant, the exact fx and fy in the Cartesian coordinate are not constant. Similar to [71], the size of the flow domain is [−20,20] × [−20,20] in the simulations, which is sufficiently large so that the exact solution for the infinite domain could still apply. The flow domain is discretized non-uniformly except in the vicinity of the cylinder, [−2.5,2.5]×[−2.5,2.5], where it is discretized uniformly with the grid spacing ∆x. The cylinder is discretized so that the spacing of the Lagrangian grid ∆s ≈ 1.3∆x. Both Tikhonov-type and Lavrentiev-type regularizations are considered with various regularization parameters λ. To investigate the convergence, simulations with various grid resolutions ∆x are conducted for each λ. The time step is selected to be 1.25×10−4 based on the finest spatial grid of all these simulations. The problem is not steady, and we choose the results at t = 2 for all quantitative investigations. First, we fix ∆x = 0.02 and choose different regularization parameters λ to investigate the regularization effect and try to find the suitable λ for this specific problem. Figure 3.13 shows the L-curve, which could help determine the suitable λ. The boundary forces are projected to the tangential direction of the cylinder and plotted in Figure 3.14. It can be observed that with sufficiently large regularization parameters λ as dictated by the L-curve, the regularized method produces a boundary force that is closer to the exact value. (For the problem, the exact tangential boundary force is a constant of approximately 0.3004.) 10 -6 10 -5 10 -4 10 -7 10 -6 10 -5 10 -4 0.001 0.01 0.05 0.1 1 10 100 (a) 10 -6 10 -4 10 -6 10 -4 0.1 0.5 1 10 1001000 10000 (b) Figure 3.13: The residual norm ∥Eu∥ 2 2 vs. smoothing seminorm ˆf 2 R at different λ(a): Tikhonov-type regularization; (b): Lavrentiev-type regularization. 69 -2 0 2 0.1 0.2 0.3 0.4 0.5 0.6 0.7 (a) -2 0 2 0.1 0.2 0.3 0.4 0.5 0.6 0.7 (b) Figure 3.14: The calculated tangential boundary force at different regularization parameters λ. (a): Tikhonov-type regularization; (b): Lavrentiev-type regularization. We now investigate the convergence of the flow velocity and the tangential boundary force as the grid is refined. Figure 3.15 shows the error of the flow velocity and boundary force compared to the exact solution with both types of regularization and various regularization parameters λ. The regularization parameter λ is chosen from a wide range of several orders of magnitude, from the under-regularized to the over-regularized region as shown in the L-curve Figure 3.13. It is observed that the flow velocity converges at first order without regularization. The presence of the regularization does not affect the convergence of the velocity field. The boundary force, however, does not converge without the regularization, similar to the observation in [71]. With the regularized projection immersed boundary method and suitable choice of λ (as dictated by Figure 3.13), the tangential boundary force also shows first order convergence. 70 10-2 10-1 10-2 10-1 (a) 10-2 10-1 10-2 10-1 100 101 (b) 10-2 10-1 10-2 10-1 (c) 10-2 10-1 10-2 10-1 100 101 (d) Figure 3.15: Error of the flow field u and the boundary force f compared to the exact solution versus the grid spacing ∆x. (a)-(b): Tikhonov-type regularization; (c)-(d): Lavrentiev-type regularization. 3.7.3 Fluid-structure interaction simulation of a fixed cylinder and a flexible beam The benchmark two-dimensional fluid-structure interaction problem FSI3, proposed by Turek and Hron [79], is simulated to showcase the performance of the regularized projection method. The structure consists of a fixed cylinder of diameter D = 0.1 m and a flexible beam of length L = 0.35 m and height H = 0.02 m that is attached behind the cylinder, and the structure is submerged in a two-dimensional incompressible channel flow; see Figure 3.16 for the schematic diagram. The fluid kinematic viscosity ν = 10−3 m2/s and density ρf = 103 kg/m3 . The parabolic inflow profile with mean velocity U = 2 m/s is prescribed on the left boundary, and the top/bottom boundaries are no-slip walls. The fluid domain is discretized uniformly in the vicinity of the beam with ∆x = 0.0025m in both directions. The beam modulus E = 5.6×106 Pa, Poisson’s 71 ratio ν = 0.4 and density ρs = 103 kg/m3 . For these parameters, periodic vortex shedding generated by the fixed cylinder induces a periodic excitation on the beam, observed in many literature [79, 80, 81]. Figure 3.16: Schematic diagram of the FSI benchmark problem. E1, E2, E3, and E4 are four edges of the beam Ωb. E5 is the boundary arc of the cylinder. 3.7.3.1 Dynamic of the beam In our two-dimensional simulation, the beam Ωb is assumed in the state of plane strain where the out-ofplane (z-direction) strain εxz, εyz, εzz = 0, and the governing equation for the displacement ub(X,Y;t) and vb(X,Y;t) are ρs ∂ 2ub ∂t 2 = ∂ σxx ∂X + ∂ σxy ∂Y , ρs ∂ 2 vb ∂t 2 = ∂ σyx ∂X + ∂ σyy ∂Y , (3.51) where X and Y refer to the point in the undeformed reference configuration, and ρs is the material density. In the plane strain state, The stress tensor σi j is related to the displacement ub(X;Y;t) and vb(X;Y;t) by the relation εxx = 1−ν 2 E σxx − ν 1−ν σyy , εyy = 1−ν 2 E σyy − ν 1−ν σxx , εxy = εyx = 1+ν E σxy εxx = ∂ub ∂X , εyy = ∂ vb ∂Y , εxy = εyx = 1 2 ∂ub ∂Y + ∂ vb ∂X . (3.52) The structural dynamic equation is solved alongside the Dirichlet boundary condition ub(X,Y;t) = vb(X,Y;t) = 0 (3.53) 72 on the fixed left edge E1 (see Figure 3.16), and Neumann boundary condition σ(X,Y;t)·nˆ(X,Y;t) = −f(X,Y;t) (3.54) on edges E2 ∪E3 ∪E4, where nˆ is the unit normal of the edge and f(X,Y;t) is the external boundary force, calculated from the fluid with the immersed boundary method. The spatial discretization of Eq.(3.51) is achieved with a FEM with quadratic triangle elements, and the maximum grid size Hmax = 0.003m (grid spacing ratio Hmax/∆x = 1.2). The finite element matrices are assembled with the help of MATLAB’s Partial Differential Equation toolbox. The second-order Newmark scheme is used for the temporal marching of Eq.(3.51). A modal decomposition is performed for the EulerBernoulli beam after the spatial discretization, and only the first three modes are considered for efficiency, similar to [82]. 3.7.3.2 The coupling method The structural boundary that the immersed boundary fluid solver sees is E2 ∪E3 ∪E4 ∪E5. For the moving boundary E2 ∪E3 ∪E4 on the beam, the location and velocity are calculated from the beam displacement [ub, vb] and velocity [u˙b, v˙b], evaluated at the boundaries. The coupling between the fluid and structure is achieved by exchanging the location/velocity and boundary force on E2 ∪ E3 ∪ E4. Herein, we use a partitioned iterative method with relaxation to achieve strong coupling between the fluid and the structure [83]; the procedure of time marching the fluid velocity u n , beam location [u n b , v n b ] and beam velocity [u˙ n b , v˙ n b ] from n to n+1 is summarized as follows: 1. Let [u old b , v old b ] ← [u n b , v n b ], and [u˙ old b , v˙ old b ] ← [u˙ n b , v˙ n b ]; 2. With [u old b , v old b ] and [u˙ old b , v˙ old b ], calculate the immersed boundary location x old bc and velocity ˙x old bc ; 3. With x old bc , ˙x old bc and u n , advance the fluid with the regularized projection IBM for a single step, and obtain new velocity u new and boundary force f new; 4. With f new, [u n b , v n b ] and [u˙ n b , v˙ n b ], advance the beam for a single step, and obtain [u new b , v new b ] and [u˙ new b , v˙ new b ]; 5. Calculate the error u˙ new b −u˙ old b + v˙ new b −v˙ old b ; 73 • If the error is sufficiently small, let u n+1 ← u new, [u n+1 b , v n+1 b ] ← [u new b , v new b ] and [u˙ n+1 b , v˙ n+1 b ] ← [u˙ new b , v˙ new b ] and step out; • Otherwise, let φ old ← αφnew + (1−α)φ old where φ = ub, vb,u˙b and ˙vb, and go to Step 2. Here, 0 < α ≤ 1 is the relaxation factor, which is chosen to be 0.1. Note that for problems with small solid-to-fluid density ratios, relaxation is necessary for the stability of iteration. 3.7.3.3 Results 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 (a) 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 (b) Figure 3.17: Flow vorticity fields and structural positions at the times when the displacements at beam tip reach maximum. Tikhonov-type regularization with λ = 0.001 is used. The FSI simulations are conducted with the Tikhonov-type regularization and parameters λ = 0.0001,0.001, and 0.01. The non-regularized projection method is also used for reference. Figure 3.17 shows the instantaneous vorticity of the flow near the structure and the displacement of the beam when λ = 0.001. The pattern of the vorticity shedding and the beam deformation qualitatively agree with the results by Turek and Hron [79], Bhardwaj and Mittal [80], and Tian et al. [81]. Figure 3.18 show the displacement of the beam tip over time after the beam reaches the state of periodic oscillation. We calculate the amplitude of the tip displacement Am, the non-dimensionalized frequency St = f D/U, and the mean drag coefficients CD = Fx/(0.5ρfU 2D) and report the results in Table 3.3 for a quantitative comparison with previous literature. The results are in reasonable agreement for a wide range of regularization parameters. 74 60 70 80 90 -0.5 0 0.5 Figure 3.18: Time histories of the beam tip displacement for the regularized projection formulation. Table 3.3: The amplitude, drag coefficients, and the Strouhal number of the benchmark FSI problem, with a comparison to previous literature Source Am/D St CD Present work, No regularization 0.38 0.24 2.40 Present work, Tikhonov-type regularization, λ = 0.0001 0.34 0.24 2.33 Present work, Tikhonov-type regularization, λ = 0.001 0.37 0.25 2.34 Present work, Tikhonov-type regularization, λ = 0.01 0.44 0.26 2.40 Turek et al. [79] 0.36 0.26 2.30 Bhardwaj et al. [80] 0.41 0.28 2.20 Tian et al. [81] 0.32 0.29 2.16 3.8 Conclusions A regularized projection immersed boundary method is proposed. The method uses the generalized Tikhonov regularization technique, and the graph Laplacian is used to construct the smoothing seminorm. Compared to the original projection immersed boundary method, the regularized method allows a more flexible choice of Lagrangian-to-Eulerian grid spacing ratios and could lead to converging boundary force solutions. The selection of the regularization parameter is not stringent and can be chosen from a wide range of values, and 75 this chapter presents a preliminary dimensional analysis and the method of L-curve to help select the parameters. The efficacy and accuracy of the regularization are demonstrated with various problems, including a fluid-structure interaction problem. 76 Chapter 4 Computational framework of fluid-structure interaction simulations In this chapter, we introduce the numerics and implementations of the fluid-structure interaction (FSI) solver developed for parachute problems, characterized by rapid deformations and not-too-small solid-to-fluid density ratios. The projection IB method with regularization for compressible viscous flow, present in Chapter 2 and Chapter 3 and articles [70, 84], will be combined with the patch-based structured adaptive mesh refinement (SAMR) for the fluid solver subsystem of the partitioned FSI solver. An efficient, partitioned, loosely coupled time integrator will be used for the coupled system. The FSI solver will be adapted to non-inertial frames of references for the complete parachute inflation, deceleration, and descent simulations. Test cases for validations will be included. 4.1 Introduction FSI simulations require resolving the fluid and structure dynamics simultaneously with kinematic and dynamic constraints between fluid and structure on the interface. On the fluid side, the IB methods have gained popularity recently; see [24, 25, 85] for reviews. The IB method is generally used with adaptive mesh refinements (AMR) so that the fluid grid resolution can be dynamically adjusted based on the resolution requirements, e.g., in turbulent boundary layers. To the best of our knowledge, however, there is limited literature that combines the projection IB method and the AMR due to the inconsistency of grid resolutions across levels. We developed a numerical implementation combining the projection IB method and the AMR. Here, we choose the patch-based structured AMR1 [86, 87] where cells are clustered into rectangular patches. The patch-based AMR grid consists of a collection of patches where data is stored as multidimensional arrays. On the structure side, the Kirchhoff-Love plate theory that considers stretching 1The other type of AMR is cell-based, which has been used with ghost fluid methods [15] 77 and bending is used, and the finite element method based on subdivision surface [88, 89] is used for spatial discretization. An explicit contact search and resolving algorithm that preserves momentum and kinetic energy [90] is integrated into the structure subsystem to handle possible fabric collision. The coupled system can be integrated with monolithic approaches by blending the fluid and structural dynamic into one monolithic solver with unitary temporal treatment [24]. This approach does not introduce additional stability conditions but requires solving nonlinear equations with both fluid and structure contributions, generally with Newton iterations, which is definitely unfavorable for parachute problems. The alternate is the partitioned approaches, where the fluid and structural solvers are treated as black boxes as the building blocks of the FSI solver. Partitioned approaches can be classified into weakly coupled (a.k.a. loosely coupled) and strongly coupled approaches. Strongly coupled approaches enforce all coupling constraints simultaneously via fixed point iteration with relaxation, while loosely coupled approaches advance the structural and fluid variables in a staggered manner without sub-iteration, but kinematic or dynamic coupling constraints are not enforced simultaneously: the structure solver provides kinematic information for the fluid solver to update and the fluid solver provides the boundary force for the structure solver. The apparent unbalance of energy transfer between fluid and structure as a result may bring instabilities, especially for small solid-to-fluid density ratio when the added-mass component is significant [91, 92]. However, it is suitable for many aeronautical applications, e.g., [32], where the solid-to-fluid density ratio is generally not too small to trigger the instabilities that blow up the simulation. The nonlinearity associated with fluid dynamics, structural dynamics with large deformation, and self-contact make the iterative approaches highly unfavorable in parachute problems. We will present an efficient, iterative-free, loosely coupled time integrator with minimal modification to existing fluid and structure solver subsystems. This chapter is organized as follows: in Section 4.2, we will briefly reiterate the plate theory and the structure solver of [88, 89]. In Section 4.3, we will introduce the implementation of the projection IB method, introduced in Chapter 2 and Chapter 3, with patch-based structured adaptive mesh refinement (AMR). The weak coupling scheme for the FSI problem is introduced in Section 4.4. We will briefly introduce the rigid body dynamics and adapt the algorithm to non-inertial frames of coordinates in Section 4.5. Some test problems for verifications and results will be given in Section 4.6. 78 4.2 Thin-shell structure dynamics In the plate theory, the configuration of the structure is described by its mid-surface parametrized by θ1 and θ2 with position xs(t,θ1,θ2) and initial undeformed reference position x¯s(θ1,θ2). The surface basis vector aα = ∂xs ∂ θα and a¯α = ∂x¯s ∂ θα for α = 1,2 and deformed and undeformed configuration, respectively. In the threedimensional body, the deformed position r and the undeformed position r¯ of material points is assumed by extruding the mid-surface along its normal direction a3 = a1 × a2/|a1 × a2| and a¯3 = a¯1 × a¯2/|a¯1 × a¯2| by θ3 ∈ [−h/2,h/2] where h is the constant thickness of the shell, i.e., r(t,θ1,θ2,θ3) = xs(t,θ1,θ2) + θ3a3(t,θ1,θ2) and r¯(θ1,θ2,θ3) = x¯s(θ1,θ2) + θ3a¯3(θ1,θ2). The covariant base vectors in deformed and undeformed configurations are g¯α = ∂ r¯ ∂ θα = a¯α +θ3 ∂a¯3 ∂ θα , g¯3 = ∂ r¯ ∂ θ3 = a¯3, (4.1) gα = ∂ r ∂ θα = aα +θ3 ∂a3 ∂ θα , g3 = ∂ r ∂ θ3 = a3, (4.2) (4.3) for α = 1,2. The strain tensor is computed according to Ei j = 1 2 (gi · gj −g¯i · g¯ j), (4.4) and truncating the strain tensor to the first order of θ3 leads to the linear strain expression Ei j = εi j +θ3κi j, (4.5) where the stretching strain εi j and the bending strain κi j are εi j = 1 2 (ai · aj −a¯i · a¯ j), (4.6) κi j = ai · ∂a3 ∂ θj −a¯i · ∂a¯3 ∂ θj = ∂a¯i ∂ θj · a¯3 − ∂ai ∂ θj · a3, (4.7) which are all related to the mid-surface deformation. The stress tensor is calculated with Young’s modulus and Poisson ratio, assuming isotropic linear elasticity, and so is the volumetric energy density. Integrating the energy density over the θ3 direction from −h/2 to h/2 leads to the energy density W(εi j,κi j) per area 79 for stretching and bending. The total elastic energy is a functional of the displacement u = xs − x¯s of the mid-surface, obtained by surface integration Πint[u] = R SW(εi j[u],κi j[u])dS, and the dynamics equations are obtained by the Euler-Lagrangian equation: Z S ρhu¨δudS+ DuΠ int[u],δu = Z S F ext s · δudS, (4.8) for all δu, where DuΠint[u],δu is the first variation of elastic energy functional along the virtual displacement δu and F ext s is the external force function. By introducing the subdivision surface shape functions for the displacement as trial and test functions, the above weak form is discretized and written in terms of the position vectors xs at nodal points: Msx¨s + f int s [xs ] = Q e→n s f ext s , (4.9) where Ms is the mass matrix, f int s [xs ] is the internal force functional and f ext s is external force given at the elements. The matrix Q e→n s applies the f ext s in the physical space of elements to the proper nodal spaces by quadrature. The explicit Newmark method is used for temporal marching but with added explicit contact detection and resolving technique [90]: x n+1 s = x n s +∆tx˙ n s + 1 2 ∆t 2 x¨ n s , (4.10) ˜x˙ n+1 s = C [x˙ n s + (1−γ)∆tx¨ n s ], (4.11) Msx¨ n+1 s + f int s [x n+1 s ] = Q e→n s f ext s , (4.12) x˙ n+1 s = ˜x˙ n+1 s +γ∆tx¨ n+1 s , (4.13) where Eqs.(4.10)-(4.11) are the prediction step and Eq.(4.13) is the correction step. C [·] denotes the fully elastic contact operator that preserves the kinetic energy and momentum used to modify the predicted velocity. The temporal time step is adaptive, so the prediction of position in Eq.(4.10) does not create a new contact. 80 4.3 Immersed boundary method with patch-based structured adaptive mesh refinement (IB-SAMR) The immersed boundary methods are generally used alongside the adaptive mesh refinement (AMR) technique so that the resolution of the fluid grid can be appropriately adjusted at regions of interest without pre-generated grids. AMR is crucial for immersed boundary simulations of wall-bounded flows with moderate to high Reynolds numbers to avoid exorbitant numbers of grid points. In this section, the patch-based structured AMR [86, 87] is integrated with the regularized projection IB method aforementioned in previous chapters as IB-SAMR. The IB-SAMR employs Structured Adaptive Mesh Refinement Application Infrastructure (SAMRAI) [93, 5, 4] developed by Lawrence Livermore National Laboratory (LLNL). 4.3.1 IB-SAMR grid structures The computational grid in SAMR is a hierarchy of levels 0 ≤ l < lmax with integer refinement ratio rl = ∆xl/∆xl+1 between the coarser level l and finer level l +1. The base level 0 covers the entire computational domain, while each finer level l + 1 covers only a subset of the coarser level l for proper nesting. Each level comprises non-overlapping rectangular grid patches dispatched to unique MPI ranks; see Figure 4.1 for a typical schematic layout. In this section, H r = {L r l } lmax−1 l=0 denotes the grid hierarchy for MPI rank r, where L r l = {Pr p}p denotes the patch level l, and Pr p denotes the patch p. Each rectangular patch Pr p has linear index extents {I 0 min,...,I 0 max −1} × {I 1 min,...,I 1 max −1} × {I 2 min,...,I 2 max −1} and coordinate range [x 0 min, x 0 max] × [x 1 min, x 1 max] × [x 2 min, x 2 max] in the corresponding three directions. The index of cells follows the convention that the finer cells in level l +1 corresponding to the cell (i 0 ,i 1 ,i 2 ) in level l are {rl i 0 ,...,rl i 0 + rl − 1} × {rl i 1 ,...,rl i 1 + rl − 1} × {rl i 2 ,...,rl i 2 + rl − 1}. Solution data is associated with each patch as contiguous multi-dimensional arrays, allowing existing non-AMR codes to be used. Numerical procedures are applied to each patch of the grid hierarchy with possible further thread-level parallelization among the patch’s index range. 81 (a) (b) Figure 4.1: Schematic layout of three-level patch-based AMR grids [4]. (a) the composite of mesh; (b) the AMR mesh hierarchy as sequence of levels. Illustrations from SAMRAI documents by LLNL. Data transferring between patches is needed when filling ghost cells at each patch’s boundary, overwriting data on coarser levels with that on finer levels, or transferring data from the old grid to the new grid after re-gridding. SAMRAI maintains the AMR grid topologies, including the patch connectivities, and performs these operations when requested. Refining or coarsening functions will be invoked when patches of two levels are involved. In addition to the SAMR Eulerian grid hierarchy for the fluid, the immersed boundary is represented with markers α ∈ {0,...,N L−1}, and solution data of the Lagrangian grid is one-dimensional array. The number of Lagrangian grid points is considered small, and for efficiency, the Lagrangian data has an identical local replica instead of a global distribution. 4.3.2 IB-SAMR grid hierarchy generating and regridding The cells tagged for refinements alongside some untagged cells are clustered into boxes, partitioned into non-overlapping patches, and dispatched to MPI ranks in the patch-based AMR; see Figure 4.2. SAMRAI handles the clustering with the Berger-Rigoustos algorithm [94] and the partitioning and dispatching with uniform load balancing and minimizing communications overhead between MPI ranks. The re-meshing process may re-arrange the patch layout and remove some finer cells when their corresponding coarser cells are not tagged, and therefore, no tagging for coarsening is needed. Solution data are transferred from the old to the new grid after re-meshing. Regridding is performed every several steps instead of every step. 82 (a) (b) (c) Figure 4.2: Re-meshing process in patch-based AMR framework [5]. (a): the cells are tagged for refinement by users; (b): the tagged cells and some untagged cells are grouped into rectangular boxes; (c): rectangular boxes partitioned into patches and dispatched among MPI ranks. Illustrations from SAMRAI documents by LLNL. In conventional AMR applications, the refinement taggings are based on the jump, the gradient, and the Hessian of flow variables. This captures the flow features of interest, e.g., the shock or the boundary layers. With the IB method, a proximity-based tagging criterion based on the distance of grid cells and the immersed boundary is adopted to ensure sufficient resolution in the boundary layer and proper grid spacing relation. Gradient-based tagging criteria: Flow variables gradient is calculated for grid cells of all patches except those in the finest level, and the magnitude of the gradient is tested against fixed thresholds el for each level 0 ≤ 0 < lmax −1 for tagging, where e0 ≤ e1 ··· ≤ elmax−2. Theoretically, Higher-order derivatives, e.g., Hessian, can serve as the refinement criteria to uniformly control the upper bound of error on each AMR level since, in finite-volume and finite-difference settings, they are the constant of error terms of spatial discretization. The gradient-based detection in most cases is sufficient. This is done by looping all over all patches and cells and applying the standard two-point centered difference in each direction. Proximity-based tagging criterion: The L ∞ distance between each fluid grid cell and the immersed boundary Γ determines whether the cell will be tagged for refinement. The distance with tessellated Γ is defined by d∞(xi ,Γ) = min 0≤α<NL ∥xi −xα∥∞ , (4.14) for each grid cell multi-indices i. The distance is tested against fixed threshold hl’s for level 0 ≤ l ≤ lmax −2, to find the set of cells {i : d∞(xi ,Γ) < hl} to be refined on each level l. Noting that ∥x∥2 / √ d ≤ ∥x∥∞ ≤ ∥x∥2 for dimension d and therefore {i : d2(xi ,Γ) < hl} ⊆ {i : d∞(xi ,Γ) < hl} ⊆ {i : d2(xi ,Γ) < √ dhl}. The set {i : d∞(xi ,Γ) < hl} = S 0≤α<NL {i : ∥xi −xα∥∞ < hl} for each level l is constructed effectively by Algorithm 1 83 executed concurrently on each MPI rank r. Algorithm 1 loops over the Lagrangian markers α in the uniform Cartesian AMR grid without calculating the distance for each cell. Algorithm 1 Proximity-based cell tagging for refinement for L r l ∈ H r do for Pr p ∈ L r l do ▷ for each patch get {I 0 min,...,I 0 max −1} × {I 1 min,...,I 1 max −1} × {I 2 min,...,I 2 max −1} ▷ index extents get [x 0 min, x 0 max]×[x 1 min, x 1 max]×[x 2 min, x 2 max] ▷ coordinate range for α ∈ [0,N L ) do for d = {0,1,2} do i d α ← I d max−I d min x d max−x d min (x d α −x d min) +I d min − 1 2 ▷ fractional index of the marker α {i d α,min,i d α,max} ← { i d α −hl/∆xl , i d α +hl/∆xl } ▷ refining index range based on distance hl end for tag all cells in {i 0 α,min,...,i 0 α,max} × {i 1 α,min,...,i 1 α,max} × {i 2 α,min,...,i 2 α,max} end for end for end for This distance threshold hl on each level can be determined by the law of the wall of the boundary layers. Similar to [32]. In the IB-SAMR, the distance threshold hlmax−2 on the level hlmax−2 must be greater than the support width of the discrete delta function to ensure that the regions where the IB forcing is applied have the finest grid. 84 4.3.3 Coarsening and refining operators Figure 4.3: Local two-dimensional AMR hierarchy with refinement ratio 2 for demonstrating the linear refining and conservative coarsening. The finer cells in level l + 1 corresponding to (i, j) in level l are (2i,2 j),(2i+1,2 j),(2i,2 j +1),(2i+1,2 j +1). Solution data is interpolated from the coarser level l to the finer level l +1 when filling ghost cells of patches when no adjacent patch on the same level provides these data or when new finer grid cells are generated. Here, we utilize multilinear interpolations with 23 stencils (three-dimensional grid). For the cell (i, j, k) in level l +1, the 23 stencils set in level l are {is ,is+1}× { js , js+1}× {ks , ks+1} where is = ji− rl−1 2 /rl k , js = j j − rl−1 2 /rl k , and ks = jk − rl−1 2 /rl k . The trilinear interpolation follows u l+1 i, j,k =(1−x)(1−y)(1−z)u l is , js ,ks +x(1−y)(1−z)u l is+1, js ,ks +(1−x)y(1−z)u l is , js+1,ks +xy(1−z)u l is+1, js+1,ks +(1−x)(1−y)zul is , js ,ks+1 +x(1−y)zul is+1, js ,ks+1 +(1−x)yzul is , js+1,ks+1 +xyzul is+1, js+1,ks+1 , (4.15) where x = i− rl is + rl−1 2 /rl , y = i− rl js + rl−1 2 /rl , and z = i− rlks + rl−1 2 /rl . The multilinear interpolation is conserved except at the boundary of each patch level. 85 Solution data in the finer level l +1 will be coarsened to level l and overwrite data in corresponding cells in the coarser level l after each step or Runge-Kutta stage for synchronization. The conserved coarsening operator is by averaging u l i, j,k = 1 r 3 l rl−1 ∑ i ′=0 rl−1 ∑ j ′=0 rl−1 ∑ k ′=0 u l+1 rl i+i ′ ,rl j+j ′ ,rlk+k ′. (4.16) To ensure the property of conservation between two levels in the AMR grid hierarchy, the numerical fluxes on the boundary of levels l +1 are coarsened to the coarsen levels once computed before updating the solution variables. On the left boundary of the patch at coarser level l, the flux at cell face I 0 min −1/2, j, k is overwritten by F l I 0 min−1/2, j,k = 1 r 2 l rl−1 ∑ j ′=0 rl−1 ∑ k ′=0 F l+1 rl I 0 min−1/2,rl j+j ′ ,rlk+k ′ , (4.17) for all j, k on the left boundary. The other five patch boundaries follow a similar coarsening approach. 4.3.4 Operators between IB Lagrangian grid and SAMR Eulerian grids The projection IB method requires uniform Eulerian grid spacing in each direction, and comparable structure/fluid grid spacing is preferred, even with regularization. The IB is applied only at the finest level to enforce these constraints. The proximity-based refinement criterion ensures that the AMR grid is the best in regions where IB forcing is applied. At the end of each step or Runge-Kutta stage, the solution data at coarser levels are overwritten with data from finer levels, eliminating the need for IB treatment at coarser AMR levels. Consequently, the IB operators are constructed considering only the Eulerian grid points of the finest AMR level. The interpolation and spreading operators are locally constructed within each MPI rank r as sparse matrices E r and L r . These operators are between the Lagrangian grid, which each rank holds as an identical local replica, and the finest level Eulerian grid owned by rank r. Specifically, E r and L r represent the block column and block row of the globally distributed interpolation and spreading operators, respectively. Each Eulerian grid point on the finest level, L r lmax−1 = Pr p , is assigned a unique linear index by flattening its multi-index subscripts within its patch and adding the total number of grid points from all previous patches. The assembly of E r and L r is performed according to Algorithm 2 on each rank r without requiring MPI calls. 86 Algorithm 2 Assemblying the IB interpolation and spreading operators with SAMR base ← 0: the index base for each patch. for Pr p ∈ L r lmax−1 do ▷ loop over patches on the finest level get {I 0 min,...,I 0 max −1} × {I 1 min,...,I 1 max −1} × {I 2 min,...,I 2 max −1} ▷ index extents get [x 0 min, x 0 max]×[x 1 min, x 1 max]×[x 2 min, x 2 max] ▷ coordinate range for α ∈ [0,N L ) do for d = {0,1,2} do i d α ← I d max−I d min x d max−x d min (x d α −x d min) +I d min − 1 2 ▷ compute fractional index of the marker α {i d α,min,i d α,max} ← { i d α −h/∆xlmax−1 , i d α +h/∆xlmax−1 } ▷ compute forcing ranges of marker α based on width h of discrete delta functions end for for i = {i d} 2 d=0 ∈ {I d min,...,I d max −1} ∩ {i d α,min,...,i d α,max} do xi = {x d i } 2 d=0 : x d i ← x d min + (i d −I d min +1/2)∆xlmax−1 ▷ compute Eulerian point location for the uniformly spaced patch grid ie ← base +∑d i d ∏ d−1 k=0 (I k max −I k min) ▷ compute flattened unique id for point i in patch Pr p E p (α,ie) ← δh(xi −xα) ▷ fill sparse entry for interpolation operator E r end for end for base ← base + ∏ 2 d=0 (I d max −I d min): increment the index base end for L r ← (E r ) T diag {∆sα } (∆xlmax−1) d ▷ obtain spreading operator by transposing and right scaling interpolation operator To enforce the stress and heat flux boundary condition in the finite-volume framework, the spreading and interpolation operators between the Lagrangian grid and cell faces of the Eulerian grids need to be assembled, and the assembling follows a similar algorithm. The equations for boundary source terms similar to the form EL f = Eu⋆ −x˙bc with distributed Eulerian grid and IB operators become ∑rErL r f = ∑r E ru ⋆ r −x˙bc. In each iteration of the Krylov space method, the operator E rL r is applied, followed by the all reduction operation in all MPI ranks. 87 4.3.5 IB-SAMR integration cycle The patch-based SAMR is well-suited for systems of hyperbolic conservation laws with finite volume spatial discretization and explicit time marching. The fluxes are calculated for all uniformly spaced rectangular patches of all levels independently after ghost cells of each patch are filled by transferring solution data from neighboring patches or interpolating solution data from coarser levels. The fluxes in level boundaries are coarsened according to Eq.(4.17) to enforce conservation at the interface of two AMR levels before being used to update the conserved variables. The maximum time step is determined by CFL and diffusion number conditions of all levels and the moving speed of IB ensuring ∆t∥x˙bc∥∞ < ∆xlmax−1, although non-uniform time steps between hierarchy levels are possible due to different grid spacing. Solution data in coarser levels are overwritten with finer levels after being updated with numerical fluxes in each Runge-Kutta stage. The IB-SAMR introduces the Lagrangian grid and enforces immersed boundary conditions after obtaining stress/heat flux and updating the solution variables, as explained in Chapter 2. Here, we denote the discretized conserved variables U n l,p = {U n l,p,i, j,k } for patch p at level l of time step n, the set of all data U n l = {U n l,p } at level l, and the hierarchy U n = {U n l }. After obtaining the intermediate state U ⋆ for the entire hierarchy, the equation for the boundary source is solved, and the projection step is applied directly to U ⋆ lmax−1 . For level l < lmax −1, the IB conditions are automatically enforced with data coarsening. 4.4 Partitioned loosely coupled FSI time integrator The IB-SAMR and the subdivision surface thin-shell finite element solver are coupled for the FSI solver. Kinematic and dynamic constraints have to be satisfied on the interfaces of fluid and structure, which in the thin-plate case is just mid-surface. The kinematic constraints are Q n→e s xs = xbc and Q n→e s x˙s = x˙bc where on the right of the equal signs is the immersed boundary data for fluid, and Q n→e s averages the nodal data to the elements. The dynamic constraints from the IB methods are f ext s = −f , and f ext s is applied to the FE nodal with Q e→n s . Time marching of the FSI system is achieved by a partitioned, loosely coupled time integrator. The time integrations for the structure and fluid are independent, as described in previous chapters and sections (Newmark method for structure dynamics and half-explicit Runge-Kutta for the fluid dynamics with immersed boundary), with the exchange of coupling information per step without iterations: the advantage of the loosely coupled integration. The structure dynamic generally requires smaller ∆t, especially with the 88 contact schemes. For efficiency, several structure dynamics Newmark steps are grouped together for one Runge-Kutta step of fluid. The loosely coupled integrations and subiterations essentially bring the temporal accuracy to first order. In short, the thin-shell mid-surface position x n s and velocity ˙x n s at time n are passed to the fluid side and used to advance the fluid subsystem using the immersed boundary method with halfexplicit Runge-Kutta time integrators for all three stages. The boundary force calculated with the projection immersed boundary method from the first Runge-Kutta stages is passed from the fluid side to the structure side with a negative sign as the external force for N Newmark sub-steps. The loosely coupled discretized systems are summarized as follows with the no-slip condition as described in Chapter 2 and Chapter 3, the thermal boundary conditions parts are not included for clarity: U˜ ν = ανU n + (1−αν ) U˜ ν−1 − ∆t ∆x ∆F˜ ν−1 +L(x n bc)φ˜ S( ˜f ν−1 ) , (4.18) ˜f ν−1 = argmin∥E(x n bc)u˜ ν −x˙ n bc∥ 2 P +λ ∆t ˜f ν−1 2 R , (4.19) (4.20) for ν = 1,2,3 successively, with U˜ 0 = U n , x n bc = Q n→e s x n s , ˙x n bc = Q n→e s x˙ n s and U n+1 = U˜ 3 . The structure Newmarks are advanced with N sub-steps: x n+k/N s = x n+(k−1)/N s + ∆t N x˙ n+(k−1)/N s + 1 2 ∆t 2 N2 x¨ n+(k−1)/N s , (4.21) ˜x˙ n+k/N s = C x˙ n+(k−1)/N s + (1−γ) ∆t N x¨ n+(k−1)/N s , (4.22) Msx¨ n+k/N s + f int s [x n+k/N s ] = Q e→n s f ext,n s , (4.23) x˙ n+k/N s = ˜x˙ n+k/N s +γ ∆t N x¨ n+k/N s , (4.24) for k = 1,··· ,N, and external force f ext,n s = − ˜f 0 . The time integration and the exchange of coupling information are illustrated in Figure 4.4. 89 X n s X n+ 1 N s X n+ k N s X n+1 s N Newmark steps of size ∆t/N U n U˜ 1 U˜ 2 U n+1 3 half-explicit Runge-Kutta stages x n s , x˙ n s f ext s Figure 4.4: Loosely coupled partitioned time stepping from n to n+1 for structure data Xs and fluid data U The weak coupling is known to blow up exponentially regardless of the time step ∆t for low solid-tofluid density ratios. Here, we negate the structure displacement, contact, and variation of fluid advection for a preliminary linear stability study that focuses purely on the loosely coupled integrator and the regularized projection methods. The simplified systems with only the projection IB methods and the loose coupling scheme is ρf u˜ 0 = ρf u n +∆tL ˜f 0 , (4.25) ˜f 0 = argmin Eu˜ 0 −Q n→e s x˙ n s 2 P +λ ∆t ˜f 0 2 R , (4.26) x˙ n+1 s = x˙ n s + 1−γ N ∆tx¨ n s + N −1+γ N ∆tx¨ n+1 s , (4.27) Msx¨ n+1 s = −Q e→n s ˜f 0 . (4.28) (4.29) With the assumption that Eun = Q n→e s x˙ n−1 s from the previous step, ˜f 0 = 1 ∆t ρf [(EL) TP(EL) +λR] −1 (EL) TPQn→e s (x˙ n s −x˙ n−1 s ) = 1 ∆t ρfKQn→e s (x˙ n s −x˙ n−1 s ), (4.30) linearly depends on ˙x n s − x˙ n−1 s with K representing the solution of the regularized projection IB equations and depends only on the mesh layout. Afterward, let ∆x˙ n s = x˙ n s −x˙ n−1 s and the advancement of the structure leads to ∆x˙ n+1 s = −ρf M−1 s Q e→n s KQn→e s 1−γ N ∆x˙ n−1 s + N −1+γ N ∆x˙ n s . (4.31) 90 Stability ensues when µ = ρf M−1 s Q e→n s KQn→e s < 1, where ∆x˙ n+1 s ≤ µ(α∥∆x˙ n s ∥ + (1 − α) ∆x˙ n−1 s ) and α = (1−γ)/N. By telescoping ∥∆x˙ n s ∥ ≲ µ n−2 ( ∆x˙ 1 s + ∆x˙ 2 s ) with decaying ∆x˙ n s and bounded ˙x n s : ∥x˙ n s ∥ ≤ x˙ 0 s + n ∑ k=1 ∆x˙ k s ≲ x˙ 0 s + 1− µ n−2 1− µ ( ∆x˙ 1 s + ∆x˙ 2 s ) < ∞. (4.32) The condition µ < 1 is satisfiable for a small fluid-to-solid density ratio. Note that this is independent of the time step ∆t. 4.5 FSI in non-inertial frames of references For the parachute problems, we model the payload as a rigid body with six degrees of freedom. The external forces acting on it include aerodynamic drag and parachute riser cable tension. The position of the center of mass xcm(t) of the rigid body follows mx¨ cm = F ext, and the unit quaternion q(t) = qw(t) +qx(t)i+qy(t)j+ qz(t)k representing its rotational orientation follows the kinematic relation d dt qw qx qy qz = 0 −ωX −ωY −ωZ ωX 0 ωZ −ωY ωY −ωZ 0 ωX ωZ ωY −ωX 0 qw qx qy qz , (4.33) where the angular velocity ω(t) = ωX (t)eˆX (t)+ωY (t)eˆY (t)+ωZ(t)eˆZ(t) are given by components ωX ,ωY ,ωZ in body-fixed basis eˆX , eˆY , eˆZ. The evolving of the angular velocity follows the Eulerian dynamic equation Iω˙ +ω ×Iω = Mext c , (4.34) where ω˙(t) = ω˙X (t)eˆX (t) +ω˙Y (t)eˆY (t) +ω˙Z(t)eˆZ(t) (note that this is not ∂ω/∂t). I is the inertial operator, a constant matrix under body-fixed basis, and Mext c is the external moment w.r.t. the center of mass. The system of equations for 6DOF dynamics is solved with explicit time marching. 91 The position of point s on rigid body’s surface relative to the center of mass is Xb(s,t) = Xb(s)eˆX (t) + Yb(s)eˆY (t)+Zb(s)eˆZ(t) = xb(s,t)eˆx+yb(s,t)eˆy+zb(s,t)eˆz with Xb(s),Yb(s),Zb(s) being the time-independent coordinate under the body-fixed basis. Its position and velocity in the inertial frame are xb(s,t) = Xb(s,t) +xcm(t), (4.35) x˙b(s,t) = ω(t)×Xb(s,t) +x˙ cm(t), (4.36) Conversely, let the external boundary force on the rigid body surface be f ext b (s), the total force F ext and momentum w.r.t. center of mass Mext c due to the surface force are F ext = Z f ext b (s)ds, Mext c = Z Xb(s)×f ext b (s)ds. (4.37) The surface kinematics and forces will be used to couple with the fluid. Figure 4.5: The inertial, non-inertial (computational), and body-fixed frame of the rigid body. xf and xcm are the origins of the non-inertial and body-fixed frame in the inertial frame, respectively The FSI framework is adapted to non-inertial frames of reference by incorporating inertial forces into the structure and fluid dynamics. In this context, the non-inertial frame undergoes translation but no rotation relative to the inertial frame, with the translation determined to keep the center of mass of the rigid body fixed. See Figure 4.5 for a schematic demonstration. The rigid body’s center of mass in the non-inertial computational frame, x f cm = xcm(t)−xf(t), (4.38) 92 is a pre-set constant, and this determines the position of the non-inertial frame, as well as the velocity and acceleration, noting that x˙ f(t) = x˙ cm(t) and x¨ f(t) = x¨ cm(t). The fluid and structure solvers are implemented in the non-inertial frames, while the rigid body is simulated in the inertial ground frame. X n b ········· X n+1 b Explicit Time marching of rigid body dynamics X n s X n+ 1 N s X n+ k N s X n+1 s N Newmark steps of size ∆t/N U n U˜ 1 U˜ 2 U n+1 3 half-explicit Runge-Kutta stages I II III IV V V I Figure 4.6: Partitioned time stepping from n to n+1 for structure data Xs , fluid data U and the data of rigid body Xb. Data exchange process I and III: passing position/velocity of the body/structure to the fluid as the input of IB solvers; II and V: passing boundary force from IB solvers to the structure and body; IV: passing position/velocity from body to structure; V I: passing force from the structure to the body. The time integrator for the coupled dynamics of rigid body, fluid, and structure is illustrated in Figure 4.6, with an exchange of coupling information one with each other once per step. The time integrator is implemented as follows: 1. The non-inertial frame velocity ˙x n f = x˙ n cm and acceleration ¨x n f = x¨ n cm based on the rigid body’s data at step n 2. Compute the position x n b and velocity ˙x n b on the surface of the rigid body at time n in the non-inertial frames from Eqs.(4.35), (4.36), and (4.38) x f,n b = x n b −x n f = X n b +x f cm, x˙ f,n b = x˙ n b −x˙ n f = ω n ×X n b ; (4.39) where x f cm is the fixed value of the center of mass at the non-inertial frame. 3. Advance the fluid solver with half-explicit RK methods, from Eqs.(4.18) and (4.19). Note that the immersed boundary position x n bc and velocity ˙x N bc here includes both the rigid body surface and the structure. Translational inertial force −ρx¨ n f is added as source terms to the momentum and energy equation. 93 4. Transfer the boundary force acting on the structure f ext,n s = −f n s a rigid body f ext,n b = −f n b , where f n s and f n b are implicitly solved boundary source with the projection method at the first RK stage. 5. The interface between the structure and the rigid body is the controlling points. Transfer the position/velocity of controlling points on the rigid body to the structure (process IV), and the interaction force at these controlling points between is passed from the structure to the body (process V I) for the interaction between. 6. The structure is advanced with f ext,n s from the fluid, and the position and velocity at controlling points are overwritten from that of the rigid body. Meanwhile, the rigid body is advanced with f ext,n b of the fluid and the interaction force at these controlling points. 4.6 Verifications 4.6.1 Airdrop of sphere Airdrop of a sphere of diameter D = 1 m under gravity and aerodynamic drags is simulated as a validation of the IB-SAMR solver under the non-inertial frame of reference. Here, we choose the parameters so that the terminal airdrop velocity u∞ = 20 m/s and Reynolds numbers Re = ρu∞D/µ = 200 and 1,000, in the region of steady and unsteady asymmetric vortex shedding. The density ρ∞ = 1 kg/m3 . The drag coefficients at Re = 200 and 1,000 are approximately 0.78 and 0.47 from previous studies [95, 96, 97], and the masses of the sphere can be inferred from the mg = 1 2 CDρ∞u 2 ∞πD 2/4 for both Reynolds numbers under g = 9.8 m/s, which are 12.50 kg and 7.54 kg. The inertial matrices of both cases are diagonal with Ixx = Iyy = Izz = 2.0833 kg m2 and 1.2567 kg m2 , respectively. The fluid is initially quiescent, and the rigid body has zero velocity and angular velocity. The computational domain for fluid is [−5,15] × [−5,5] × [−5,5] m with the sphere fixed at (0,0,0) m. The AMR hierarchy has five levels with ∆x/D = 0.0125 on the finest level. Fifth-order WENO scheme is used for the advection fluxes, while the five-point fourth-order centered difference scheme is used for the viscous fluxes. The velocity and drag with time for cases with expected terminal Re = 200 and Re = 1,000 are plotted in Figure 4.7 and Figure 4.8, respectively, with the expected terminal velocity u∞ = 20 m/s marked for comparison. The aerodynamic drag force increases monotonically for the case with expected Re = 200, while the drag force for the case with expected Re = 1,000 shows some roughness when the velocity becomes 94 large enough due to the unsteady flow pattern. For the latter case, the unsteady asymmetric vortex shedding structure is plotted in Figure 4.9 at a typical instance. The expectation and actual values at terminal states are tabularized in Table 4.1 for a closer scrutinizing, showing good agreements. 0 2 4 6 8 10 12 -25 -20 -15 -10 -5 0 Terminal Velocity (a) 0 2 4 6 8 10 12 0 20 40 60 80 100 120 140 Gravity (b) (c) Figure 4.7: The properties of sphere airdrop at Re = 200. (a): the sphere dropping velocity; (b): the drag force on the sphere; (c): fluid velocity near the ball at terminal state. 0 2 4 6 8 10 -25 -20 -15 -10 -5 0 Terminal Velocity (a) 0 2 4 6 8 10 0 20 40 60 80 100 120 140 Gravity (b) (c) Figure 4.8: The properties of sphere airdrop at Re = 1,000. (a): the sphere dropping velocity; (b): the drag force on the sphere; (c): fluid velocity near the ball at terminal state. 95 (a) (b) Figure 4.9: The vortex structure near the sphere at the terminal stage for Re = 1,000. (a) for vorticity at direction normal to the paper; (b): contour of vorticity magnitude. Table 4.1: The expected and actual Reynolds number, terminal velocity, and drag coefficients at terminal states. Expected Actual Value Re u∞ (m/s) CD Re u∞ (m/s) CD 200 20 0.78 198.4 19.84 0.79 1000 20 0.47 992 19.84 0.478 96 4.6.2 Drop of Parachute Test Vehicle (PTV) (a) (b) Figure 4.10: The geometry of the PTV capsules for CPAS airdrop test The Parachute Test Vehicle (PTV) is the full-scale capsule-shaped test vehicle designed by the Orion Capsule Parachute System (CPAS) team to mimic the actual shape and weight of Orion spacecraft for parachute test purposes [3]. The airdrop of the PTV under gravity and aerodynamic drags is simulated as the case of another validation for the IB-SAMR framework under the non-inertial frame of reference. Here, the ambient air properties are taken from the US standard atmosphere data at 4000m above sea level, i.e., ρ = 0.8194 m/s, p = 61660 Pa and µ = 1.66×10−5 N m/s 2 . The PTV is a rigid body with mass m = 9979.02 kg and the inertial matrix in the body-fixed frame at the center of mass I = 20777.4 −23.4112 877.919 −23.4112 16680.5 117.056 877.919 117.056 13754.1 kg m2 . (4.40) The geometry is illustrated in Figure 4.10 with D ≈ 5 m. The origin of the coordinate is the center of mass, which is not the centroid of the geometry. The blunt end is designed to face downward in the airdrop, and the position of the center of pressure is behind the center of gravity, which makes the PTV aerodynamically 97 stable. The gravity g = −9.8 m/s 2 in the x direction. The fluid is initially quiescent, and the rigid body initially has zero velocity and angular velocity. The initial Euler angles (in Tait–Bryan z−y ′−x ′′ convention) that describe the orientation of the body are yaw α = 0 ◦ , pitch β = −45◦ , and roll γ = 0 ◦ , i.e., the normal direction of the blunt end (eˆX direction in the body-fixed frame) and the gravity direction (eˆx direction in the inertial frame) has an angle of arccos(cosα cosβ) = 45◦ . The computational domain for the fluid is [−50,10] × [−25,25] × [−25,25] m in three directions, and the center of mass of the rigid body is fixed at (−34.1117,−0.01524,0.14478) m. Four levels of grids with the refinement ratio 2 are used in the AMR hierarchy, and the finest AMR grid ∆x = 0.078125 m. Fifth-order WENO scheme is used for spatial discretization of the compressible N-S equations. Figure 4.11: The orientation or the PTV and the fluid velocity nearby. The snapshot was taken at a time when the wobbling angle is a local maximum. The velocity and drag force at x-direction and pitch angle β of the PTV are plotted in Figure 4.12. The PTV wobbles periodically, as depicted in Figure 4.11 during airdrop because of its initial pitch angle. The yaw angle α remains negligible compared to the pitch angle β in the airdrop and thus has negligible contribution to the total wobbling angles. Because of the wobbling, the projected area of the PTV against the mainstream flow changes with time. This is reflected in the aerodynamic drag, which oscillates periodically in synchronization with the pitch angle, as shown in Figure 4.12c and Figure 4.12b. The PTV reaches the terminal velocity at ∼ 112m/s, corresponding to the drag coefficient CD = 0.97 based on the front area. The pitch angles shows periodic oscillation near the equilibrium angle approximately β = −10◦ , as a result of the position of the center of gravity ahead of the center of pressure. 98 0 5 10 15 20 25 30 35 -120 -100 -80 -60 -40 -20 0 Terminal Velocity (a) 0 5 10 15 20 25 30 35 0 2 4 6 8 10 12 104 Gravity (b) 0 5 10 15 20 25 30 35 -50 -40 -30 -20 -10 0 10 20 (c) Figure 4.12: The kinematics and dynamics with time in the PTV airdrop simulations. (a): the x-velocity; (b): the aerodynamical drag in x-direction; (c): the pitch angles. 4.7 Conclusions The projection IB method with regularization is implemented under a parallelized patch-based adaptive mesh refinement (AMR) framework. The flow features and boundary proximity are refinement criteria of the grid for regions of interest. The fluid dynamic solver featuring IB and AMR is incorporated with a thin-shell finite element solver of Kirchhoff-Love plate as a loosely coupled fluid-structure interaction (FSI) solver where the time integrator is done in a staggered manner. The FSI solver can handle non-inertial frames of references, and a six-degree-of-freedom rigid body dynamic solver is incorporated in the noninertial FSI solver to keep track of moving objects. The solver is verified with problems of a falling sphere to demonstrate the correctness, and the falling of a spacecraft vehicle with wobbling is investigated. 99 Chapter 5 Simulation of the drogue parachute cluster of the Orion Capsule Parachute Assembly System (CPAS) 5.1 Introduction Parachutes are critical devices in space exploration missions that affect payload vehicles’ Entry, Descent, and Landing (EDL) process. They provide aerodynamic drag for deceleration, set specific terminal descent speeds, provide height and timeline for scientific instruments, provide stability that prevents tumbling, and deploy further parachutes in the system [9, 10, 3]. The parachute system for space vehicles can consist of a single stage or multiple stages that are deployed in sequence, and each stage can consist of a single parachute or a cluster of smaller parachutes for increased stability. e.g., the Viking mission to Mars used a single disk-gap-band parachute, while the Venus and Galileo missions used a two-stage parachute system with a pilot parachute that deployed the conical ribbon main parachute due to the thicker atmosphere and larger aerodynamical loads needed. Today, the newly designed human-rated Orion Capsule Parachute Assembly System (CPAS) for Earth re-entry consists of four stages and 11 parachutes. Some of the parachutes have multi-stage reefing that restricts the skirt opening and prevents sudden over-inflation. Summaries of experiments or computational studies of the CPAS system can be found in [10, 13, 6, 7, 16, 3]. Computational studies of parachutes and associated fluid-structure interaction phenomena could provide practical guidance for design and validation. The arbitrary Lagrangian-Eulerian (ALE) method is traditionally used in FSI problems (e.g., [98, 99, 100]). However, in recent literature, the immersed boundary method without dynamic regridding is more frequently used (e.g., [52, 101, 14, 32]). Most computational studies of parachutes aim to capture the dynamic features during deceleration. Few focus on the initial inflation and disreefing despite having more profound dynamic behaviors and larger peak load and stress. In this chapter, 100 we investigate the dynamics of parachutes computationally and focus on the unexplored part based on the FSI solver in Chapter 4. The second stage of the CPAS consists of two mortar-deployed drogue parachutes with conical ribbon shapes designed to decelerate and stabilize the payload. These drogue parachutes undergo two-stage reefing, which helps prevent large peak loads and over-inflation. This chapter will focus on the CPAS drogue parachute cluster. The setup of the problems will be described in Section 5.2, including the fluid/structure properties, flight conditions, and reefing conditions. Preliminary results, including the transient dynamics after inflation and disreefing, will be reported in Section 5.3. 5.2 Problem setup 5.2.1 Parachute system descriptions In the final and qualified design of CPAS, two drogue parachutes are used in the second stage [3]. The drogue parachutes are Variable Porosity Conical Ribbon parachutes with diameter 23 ft and cone angle 25.7 ◦ in its conical constructed shape [6, 7]. See Figure 5.1 for illustrations. The drogue parachute has 24 gores, which are segments of the parachute canopy separated circumferentially by radial reinforcement tapes. Each gore consists of 52 horizontal ribbons from the parachute skirt to the vent, with gaps between two ribbons that provide geometric porosity. Three of these gaps between the ribbons are wider gaps, and the three gaps separate the ribbons into four groups with different materials. The drogue parachute is reinforced with webbing bands on the skirt and vent, and 24 reinforced radial tapes are used to stitch the ribbons between gores. 24 suspension lines are stitched to the radial tapes and combined to a riser line that is attached to the payload. Two sets of reefing lines are placed along the parachute skirt to constrain its diameters. The reefing lines will be cut after some time to disreef the parachute. 101 (a) (b) Figure 5.1: (a): CPAS drogue parachute constructed geometry [6, 7]. (b): the typical construction of a conical parachute, the gore consists of multiple ribbons from the skirt to the vent [8]. The drogue geometry used in this chapter is simplified. Most notably, the small gaps between ribbons are neglected, and only the three wider gaps are considered. The four groups of ribbons separated by wider gaps are combined into larger fabric bands in the model. The radial tapes and the skirt/vent bands are treated as one-dimensional linear elastic tensile cords, the same as the suspension, riser, and reefing lines. See Figure 5.2 for the initial drogue shape used for this chapter that we have in hand, including the fabrics and tensile cords. They should be isometrical to the nominal construction. The ribbon groups from the vent to the skirt are numbered from 1 to 4, and the material properties in our studies are tabularized in Table 5.1. The properties of lines, tapes, and bands used in our studies are tabularized in Table 5.2. The materials data are taken from private communications and published data in [102]. The length of each suspension line is 14.0208 m, exactly twice the nominal diameter of the drogue parachute as used in the CPAS Gen II [7]. The length of the riser line is 16.34 m. 102 (a) (b) (c) (d) Figure 5.2: The initial drogue shape used in this chapter. (a) and (b): the top and side view of the drogue parachute. The 52 ribbons are grouped into four larger bands in the simplified model. (c) and (d): the riser line (yellow), suspension lines (blue), reefing lines (green), skirt/vent bands and radial tapes (purple) of the drogue. Table 5.1: Materials properties of ribbons. The density and Poisson ratio for all ribbons are 1.15 g/cm3 and 0.33. Description Material Thickness (mm) Young’s Modulus (GPa) ribbons group 1 Nylon type 5 Class D 0.06858 0.8310 ribbons group 2,3 Nylon type 5 Class C 0.08128 0.9669 ribbons group 4 Nylon type 2 Class B 0.18288 1.1805 103 Table 5.2: Materials properties of tensile cords. Description Material Density (kg/m) Stiffness (N) Riser lines 24x Kevlar Cord, Type XB 0.42 1.5148×107 Suspension lines Kevlar Cord, Type XB 0.01750 6.31166×105 Reefing lines Kevlar Cord, Type XI 0.02480 7.81444×105 Radial tapes, skirt/vent bands PIA-T-87130 Type VI Class 9 Kevlar Webbing 0.03100 7.11715×105 The two reefing lines at the skirt corresponding to the two reefed stages restrict the drogue parachute’s opening diameter before being cut. In our studies, the reefing line lengths are 329 in and 437 in, respectively, which correspond to the reefing ratios 38% and 50% based on the constructed nominal conical shape with D = 23 ft shown in Figure 5.1, and reefing ratios 67% and 89% based on the initial undeformed shape used in simulations Figure 5.2. Two drogue parachutes are used in the second stage of CPAS. In our studies, the two drogue parachutes with their respective suspension and riser lines are exactly the same as described above, and the two riser lines are anchored to the same point of the payload. In the initial configuration, the two drogue parachutes are rotated at an angle so they do not touch initially; see Figure 5.3 for the setup. The dynamic parameters of the PTV are given in Section 4.6.2. . Figure 5.3: The initial setup of two drogue parachutes and the PTV payloads 104 5.2.2 Flight conditions The flight conditions of the simulations are determined by the altitude-Mach number trajectories of CPAS airdrop experiments [3]. Here, we choose the flight conditions of altitudes 20,000 ft, 16,000 ft, and 13,000 ft and Mach numbers 0.3, 0.25, and 0.2 as the initial conditions of simulations of reefed stage 1, 2, and full open stage 3 without reefing line, respectively. The ambient density and pressure are determined by 1976 US Standard Atmosphere data, and the air is treated as idea gas with specific gas constant 287.058 J/(kg K) and adiabatic ratio 1.4. The dynamic viscosity is chosen as constant µ = 1.661×10−5 N s/m2 and Prandtl number cpµ/κ = 0.7. The data are tabularized in Table 5.3 as FC1, FC2, and FC3. Simulations at each stage are conducted independently since a unified simulation that includes all three stages and the entire descent process is unnecessarily expensive. For simplicity, the ambient density and pressure do not change during the descent; we only track the velocity changes. Table 5.3: Initial flight conditions of simulations of stages 1, 2, and 3 Flight condition Altitude (ft) Mach Density (kg/m3 ) Pressure (Pa) Velocity (m/s) FC1 20,000 0.3 0.6527 46563.3 94.81 FC2 16,000 0.25 0.7460 54915.2 80.26 FC3 13,000 0.2 0.8224 61942.9 64.95 5.2.3 Mesh information Figure 5.4: triangular subdivision surface FE mesh of the CPAS drogue parachutes cluster 105 The subdivision surface thin-shell finite element mesh of the CPAS drogue parachutes cluster used in the FSI simulations are shown in Figure 5.4. The mesh consists of 14,112 vertices and 25,728 triangular elements in total, with equal numbers of vertices and elements for each drogue parachute. The riser line, suspension lines, skirt/vent bands, and radial tapes are treated as tensile cords and discretized as collections of linear elastic links, each connecting two vertices on the FE mesh or independent vertices of the cords, and the mass of cords are transferred to the vertices. The suspension lines consist of 1,008 link elements and the riser lines consist of 2 link elements. The computational domain of the fluid is [−60,40] × [−40,40] × [−40,40] m with the PTV center of mass fixed at (−34.1117,−0.01524,0.1447) m, noting that the geometric center of PTV (which does not colocate with the center of mass) is positioned at y = 0 and z = 0 initially. The PTV anchor point is initially (−32.6769,0,0) m on the geometric central axis. The AMR grid hierarchy consists of four patch levels with a refinement ratio of 2 in between. On the finest level ∆x = 0.078125 m. 5.2.4 Simulation schedules Stages 1, 2, and 3 simulations are conducted independently with the same mesh and varying atmospheric conditions and initial velocities as described in previous sections. The high computational costs do not warrant a single, unified simulation of the entire descent process and all three stages: the entire descent timeframe is too long and uninteresting, with only repeated features most of the time. Here, we are interested in flow features and parachute behaviors during the initial inflation and disreefing and sometime shortly after. The drogue parachutes cluster FSI simulation schedule is illustrated in Figure 5.5. 106 Time PS1 S1 FC1 PS2 S2 FC2 PS3 S3 FC3 Stage 1 Stage 2 Stage 3 Deployment Disreefing Disreefing Figure 5.5: Simulation schedules of stages 1, 2, and 3. PS[1,2,3] and S[1,2,3] denote pre-stage [1,2,3] and stage [1,2,3] simulations conducted in order, respectively. FC[1, 2, 3] denotes the atmospheric and flight conditions in Table 5.3 for each stage. Pre-stage simulations PS1, PS2, and PS3 are conducted before non-inertial airdrop simulations S1, S2, and S3 of stages 1, 2, and 3, respectively, to provide proper initial conditions of fluid and structure for the abrupt deployment or disreefing to bootstrap S1, S2, and S3. The atmospheric condition of each simulation and their corresponding pre-stage simulations are the same. Stage 1 of airdrop starts from the moment the drogue cluster is mortar deployed, and the initial inflation process immediately follows. The precise mortar deployment process is beyond the current capability of FSI simulations, and the study in this article uses a strategy where the fluid is started with a standalone CFD simulation with a non-moving, semi-folded, stage 1 reefed drogue cluster geometry, which is the PS1 simulation. Once the flow fields of PS1 are fully developed, the non-inertial FSI simulation S1 starts, with the flow fields from PS1 and the semi-folded drogue configuration as the initial conditions. The semi-folded, reefed configuration is generated through a standalone structural mechanics simulation. In this simulation, uniform loads—approximately equal to the dynamic load at stage 1—are applied to the parachute, with a gradual reduction of the reefing ratio. After this process, the uniform loads are removed, allowing for the release of internal stress. It has to be emphasized that the semi-folded configuration generated in such a manner very likely does not represent the actual mortar deployed geometry, and a large discrepancy between the simulated and actual deployment data is possible. The PS2 and PS3 are FSI simulations under inertial frames. In these simulations, the atmospheric parameters and velocity boundary conditions are set to the constants of FC2 and FC3, respectively. PS2 and 107 PS3 use the reefing ratio corresponding to the stage 1 and 2, respectively. The non-inertial FSI simulations, S2 and S3, begin after PS2 and PS3 are fully developed. These simulations use the fluid and structural variables from PS2 and PS3 as their initial conditions. During the transition from PS2 to S2, the reefing ratio abruptly changes from stage 1 to stage 2. Similarly, when transitioning from PS3 to S3, the reefing line is completely removed. The schedules are summarized in Table 5.4. Table 5.4: Summary of the simulation schedules Case Type Flight Condition Stage Initial Condition PS1 Standalone CFD FC1 1 - S1 Non-inertial FSI FC1 1 PS1 PS2 Inertial FSI FC2 1 - S2 Non-inertial FSI FC2 2 PS2 PS3 Inertial FSI FC3 2 - S3 Non-inertial FSI FC3 3 PS3 108 5.3 Preliminary results 5.3.1 Reefed stage 1 (a) (b) Figure 5.6: Initial conditions generated of S1 generated by PS1. (a): the initial flow field velocity. (b): the semi-folded reefed drogue cluster configurations The initial configuration of S1 is shown in Figure 5.6, where the drogue parachutes are semi-folded with the reefing ratio of stage 1, intended to mimic the mortar deployment. Figure 5.7 shows the histories of drogue parachute skirt circumference (average of the two drogue parachutes) and the riser line tensions (combining the two riser lines of the two drogue parachutes). It can be observed from Figure 5.7a that the parachute cluster inflates rapidly within 0.08 s from the semifolded shape we generated under FC1. Snapshots of drogue cluster configurations during the inflation are 109 illustrated in Figure 5.8. Again, the transient inflation dynamics might not agree with the actual data due to the discrepancy between the simplified simulation schedule and the actual mortar deployment process. The parachute riser tension histories at the gravitational direction (negative x-direction in our setting) and the scatter plot of the riser tension direction at sampling time steps are plotted in Figure 5.7b and Figure 5.7c, respectively. It can be observed that in the reefed stage 1 under FC1, the riser tension fluctuates around a constant, and the riser tensions align with the gravitational direction with q |F2 y +F2 z |/|F2 x +F2 y +F2 z | less than 0.15 at all times. 0 0.5 1 1.5 2 2.5 6.4 6.6 6.8 7 7.2 7.4 7.6 7.8 8 8.2 8.4 (a) 0 0.5 1 1.5 2 2.5 0 1 2 3 4 5 6 7 8 104 (b) 0° 30° 60° 90° 120° 150° 180° 210° 240° 270° 300° 330° 0 0.05 0.1 0.15 (c) Figure 5.7: The parachute related histories in the Case S1. (a): the average skirt circumference of two drogue parachutes in the cluster. (b): tensions of the riser line at gravitional direction. (c): scatter plot of the riser tension direction at sampling steps. 110 (a) (b) (c) (d) Figure 5.8: Snapshots of inflation at initial deployment from semi-folded configurations with von Mise stress of canopy displayed. (a): at the release moment. (b): 0.0202 s after release. (c): 0.0402 s after release. (d): 0.0599 s after release. 111 The PTV-related histories are plotted in Figure 5.9. The total drag in Figure 5.9a includes the riser tension and the aerodynamic drags directly exerted on the surface of PTV. Sufficient drag is generated for small deceleration most of the time under FC1 and reefed stage 1, as shown. 0 0.5 1 1.5 2 2.5 5 6 7 8 9 10 11 12 13 104 Gravity (a) 0 0.5 1 1.5 2 2.5 93 93.2 93.4 93.6 93.8 94 94.2 94.4 94.6 94.8 95 (b) 0 0.5 1 1.5 2 2.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 (c) 0 0.5 1 1.5 2 2.5 -16 -14 -12 -10 -8 -6 -4 -2 0 2 (d) Figure 5.9: The histories of the kinematics and dynamics of the PTV payload. (a): the total drag force on the gravitional direction. (b): the velocity along the gravitional direction. (c): the yaw angle α of the PTV. (d): the pitch angle β of the PTV. 5.3.2 Reefed stage 2 PS2 and S2 simulate the transition from reefed stage 1 to stage 2 under FC2. The PS2 set the initial condition of fluid and structure under FC2 and reefing ratio of stage 1. The reefing line length is modified to stage 2 instantaneously at the beginning of S2. The drogue parachute rapidly inflates after disreefing, as depicted 112 in Figure 5.10a, reaching the circumference set by the reefing ratio of stage 2 within 0.025 s. Snapshots of parachute configurations during inflations are plotted in Figure 5.11. During the inflation, a sudden surge of the riser tension is observed as shown in Figure 5.10b, with the peak riser tension approximately 170% that of the average tension during the stable deceleration under FC2. The angle between the riser tension and gravitational direction remains small, as shown in Figure 5.10c. -1 -0.5 0 0.5 1 1.5 2 8 8.5 9 9.5 10 10.5 11 11.5 Disreefing (a) -1 -0.5 0 0.5 1 1.5 2 3 4 5 6 7 8 9 10 11 104 Disreefing (b) 0° 30° 60° 90° 120° 150° 180° 210° 240° 270° 300° 330° 0 0.02 0.04 0.06 0.08 0.1 (c) Figure 5.10: The histories of PS2 and S2 simulation before and after transition to reefing stage 2. (a): the average skirt circumference of two drogue parachutes in the cluster. (b): the tension of the riser line at gravitational direction. (c): scatter plot of the riser tension direction at sampling steps. Red and blue points are samples before and after disreefing (Case PS2 and S2). 113 (a) (b) (c) (d) Figure 5.11: Snapshots of inflation at the disreefing from reefed stage 1 to reefed stage 2 with von Mise stress of canopy. (a): at the starting moment. (b): 0.0084 s after disreefing. (c): 0.0169 s after disreefing. (d): 0.0254 s after disreefing. 114 The PTV-related histories are displayed in Figure 5.12. The riser tension and the aerodynamic drag on PTV provide sufficient drag for deceleration under FC2, as shown in Figure 5.12a, with a small deceleration rate. 0 0.5 1 1.5 2 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 105 Gravity (a) 0 0.5 1 1.5 2 79 79.2 79.4 79.6 79.8 80 80.2 80.4 (b) 0 0.5 1 1.5 2 -2.5 -2 -1.5 -1 -0.5 0 0.5 (c) 0 0.5 1 1.5 2 -15 -10 -5 0 (d) Figure 5.12: The histories of the kinematics and dynamics of the PTV payload. (a): the total drag force on the gravitational direction. (b): the velocity at gravitional direction. (c): the yaw angle α of the PTV. (d): the pitch angle β of the PTV. 5.3.3 Fully open stage 3 PS3 and S3 simulate the disreefing process from the reefed stage 2 to full open stage 3 under FC3. The PS3 set the initial condition of fluid and structure under FC3 and reefing ratio of stage 2. The reefing line is removed instantaneously at the beginning of S3. The drogue parachute rapidly inflates after disreefing, 115 as depicted in Figure 5.13a and snapshots in Figure 5.14, with significantly increased and fluctuating skirt circumference without the reefing line restricting its maximum opening. A sudden surge of the riser tension is observed as expected, shown in Figure 5.13b. The peak riser tension reaches approximately 300% that of the average tension during the stable deceleration under FC3 and disreefed drogue clusters. The riser tension fluctuates in the same trend as the skirt circumference. The fly-out angles of riser tensions increase significantly after removing the reefing line but are still small. -1 -0.5 0 0.5 1 1.5 2 11 11.5 12 12.5 13 13.5 14 14.5 15 15.5 16 Disreefing (a) -1 -0.5 0 0.5 1 1.5 2 2 4 6 8 10 12 14 104 Disreefing (b) 0° 30° 60° 90° 120° 150° 180° 210° 240° 270° 300° 330° 0 0.1 0.2 (c) Figure 5.13: The histories of PS3 and S3 simulation before and after transiting to full open stage 3. (a): the average skirt circumference of two drogue parachutes in the cluster. (b): the tension of the riser line at gravitational direction. (c): scatter plot of the riser tension direction at sampling steps. Red and blue points are samples before and after completely removing the reefing line (Case PS3 and S3), respectively. 116 (a) (b) (c) (d) Figure 5.14: Snapshots of inflation after reefing line is removed. (a): at the starting moment. (b): 0.0159 s after disreefing. (c): 0.0315 s after disreefing. (d): 0.0629 s after disreefing. The PTV-related histories are displayed in Figure 5.15. The riser tension and the aerodynamic drags on PTV provide sufficient drag for deceleration under FC3 most of the time, as shown in Figure 5.15a. 117 0 0.5 1 1.5 2 0.6 0.8 1 1.2 1.4 1.6 1.8 105 Gravity (a) 0 0.5 1 1.5 2 62 62.5 63 63.5 64 64.5 65 65.5 (b) 0 0.5 1 1.5 2 -14 -12 -10 -8 -6 -4 -2 0 2 (c) 0 0.5 1 1.5 2 -10 -8 -6 -4 -2 0 2 4 6 8 10 (d) Figure 5.15: The histories of the kinematics and dynamics of the PTV payload. (a): the total drag force on the gravitational direction. (b): the velocity at gravitional direction. (c): the yaw angle α of the PTV. (d): the pitch angle β of the PTV. 5.4 Conclusions The drogue parachute cluster of the Orion Capsule Parachute Assembly System (CPAS) is simulated with a fluid-structure interaction solver under the non-inertial frame of references, with typical flight conditions within its flight envelope and designed reefed stages. The transient dynamics of parachute deployment and disreefing involve complex multi-physics interactions, and capturing these phenomena is computationally 118 challenging. We designed the computational schedules to capture these dynamics, and preliminary simulation results are presented. Under typical flight conditions of the reefed stage 1, 2, and full open stage 3, the drogue parachutes cluster provides sufficient drags for slow decelerations of the PTV, with controlled fly-out angles and no tumbling. The peak riser tensions during the initial deployment and disreefing are significantly higher than during normal flight conditions, and it remains to be investigated whether they are within the tolerance of the materials and the designs. 119 Chapter 6 Concluding Remarks 6.1 Conclusions A projection immersed boundary (IB) method for compressible viscous flow is presented for the fluid dynamic simulation without using body-conforming grids, where the immersed boundary is represented with arbitrary unstructured surface grid, and the fluid dynamic is resolved with a Cartesian grid, with a discrete delta function transferring data in between. The IB method developed accommodates a wide set of kinematic and thermal boundary conditions, and the boundary condition is enforced with added singular source terms identified by the jump conditions. The governing equation of the projection IB method is a system of differential-algebraic equations after proper spatial discretization, where the singular source terms are algebraic variables and the boundary conditions are the algebraic constraints. The system of differential-algebraic equations is solved with half-explicit Runge-Kutta methods, and each stage is solved with a fractional-step approach. In the projection method, the singular sources on the Lagrangian grid are solved from discrete Fredholm integral equations, which are not well-posed in the sense that the discrete integral operators do not have a bounded inverse as the grid refines. Therefore, the boundary sources calculated are not accurate and show no spatial convergence, usually with spurious oscillations, especially with small Lagrangian-to-Eulerian grid spacing ratios. This fact has been noted by some researchers before. To alleviate this problem, we proposed a generalized Tikhonov regularization approach that reformulates the problem into a minimization problem. The graph Laplacian is used to formulate the seminorm similar to the Sobolev H 1 seminorm, which is used to quantify the smoothness of the solution and serves as the penalty term in regularization. Compared to 120 the original projection immersed boundary method, the regularized method allows a more flexible choice of Lagrangian-to-Eulerian grid spacing ratios and leads to converging boundary source solutions. The projection IB method with regularization is implemented under a parallelized patch-based adaptive mesh refinement (AMR) framework for efficient computation of large-scale problems, with flow features and boundary proximity as refinement criteria. The fluid dynamic solver featuring IB and AMR is incorporated with a thin-shell finite element solver of Kirchhoff-Love plate as a loosely coupled fluid-structure interaction (FSI) solver where the time integrator is done in a staggered manner. The FSI solver can handle non-inertial frames of coordinates, and a six-degree-of-freedom rigid body dynamic solver is incorporated in the noninertial FSI solver to keep track of moving objects. The IB method with regularization is validated with typical test cases, including a strongly coupled FSI problem, and the implementation of the FSI solver with the projection IB method and patch-based AMR is verified with a falling sphere problem to demonstrate the correctness. With the FSI tool readily available, the drogue parachute cluster of the Orion Capsule Parachute Assembly System (CPAS) is simulated under the non-inertial frame of references, with typical flight conditions within its flight envelope and designed reefed stages. The transient dynamics of parachute deployment and dis-reefing involve complex multi-physics interactions, and capturing these phenomena is computationally challenging. We designed the computational schedules to capture these dynamics, and preliminary simulation results are presented. Under typical flight conditions of the reefed stage 1, 2, and full open stage 3, the drogue parachutes cluster provides sufficient drag for slow decelerations of the PTV, with controlled fly-out angles and no tumbling. The peak riser tensions during the initial deployment and dis-reefing are significantly higher than during normal flight conditions, and it remains to be investigated whether they are within the tolerance of the materials and the designs. 6.2 Future outlooks Here, some possible future investigations based on the current work are enumerated. 1. Incorporating the wall modeling for large-eddy simulations. In the future, we would like to use the compressible immersed boundary method at high Reynolds numbers in large-eddy simulations. This requires wall modeling near the immersed wall. The compressible immersed boundary method can enforce both kinematic and thermal boundary conditions, 121 and we demonstrated the accuracy with test cases in low Reynolds numbers. The method naturally allows imposing the stress conditions instead of the velocity conditions, and the τw obtained from the wall model can be inserted trivially. Probing the near-wall velocity for the input of wall models, however, can be challenging, as the ”wall” in the diffused IBMs is smeared over several grid widths. 2. Parameter-free iterative Tikhonov regularization. The generalized Tikhonov regularization for projection IB method here is achieved by explicitly forming the minimization problem and constructing the left hand side operators for solving the singular source terms. A more general approach would be using iterative Tikhonov regularization when solving the linear equation with a Krylov space method. 3. Parametric studies of the drogue parachute cluster The drogue cluster is currently simulated under a single typical flight condition for each given reefed/full open stage. More investigations would be needed to verify and optimize the parachute designs with a wider range of flight conditions and reefing schedules. 122 Bibliography [1] Sigal Gottlieb, Chi-Wang Shu, and Eitan Tadmor. Strong Stability-Preserving High-Order Time Discretization Methods. SIAM Review, 43(1):89–112, jan 2001. [2] V. Brasey and E. Hairer. Half-Explicit Runge-Kutta Methods for Differential-Algebraic Systems of Index 2. SIAM Journal on Numerical Analysis, 30(2):538–552, 4 1993. [3] Carol T. Evans. Orion capsule assembly system (CPAS) test program summary. In AIAA Aviation 2019 Forum. American Institute of Aeronautics and Astronautics, 2019. [4] Brian T.N. Gunney and Robert W. Anderson. Advances in patch-based adaptive mesh refinement scalability. Journal of Parallel and Distributed Computing, 89:65–84, March 2016. [5] R.W. Anderson, W.J. Arrighi, N.S. Elliott, B.T. Gunney, and R.D. Hornung. SAMRAI Concepts and Software Design. Technical report, Lawrence Livermore National Laboratory, 2013. [6] Randy Olmstead, Aaron Morris, Kristin Bledsoe, and Megan Englert. Overview of the Crew Exploration Vehicle Parachute Assembly System (CPAS) Generation I Drogue and Pilot Development Test Results. In 20th AIAA Aerodynamic Decelerator Systems Technology Conference and Seminar, pages 1–14, Seattle, Washington, may 2009. American Institute of Aeronautics and Astronautics. [7] Aaron Morris, Kristin Bledsoe, Usbaldo Fraire, Jr., James Moore, Eric Ray, and Leah Olson. Summary of CPAS Gen II Testing Analysis Results. In 21st AIAA Aerodynamic Decelerator Systems Technology Conference and Seminar, pages 1–14, Dublin, Ireland, may 2011. American Institute of Aeronautics and Astronautics. [8] E. G. Ewing, H. W. Bixby, and T. W. Knacke. Recovery Systems Design Guide. Technical Report AFFDL-TR-78-151, IRVIN INDUSTRIES INC, GARDENA, CA, December 1978. [9] Juan R. Cruz and J. Stephen Lingard. Aerodynamic decelerators for planetary exploration: Past, present, and future. Collection of Technical Papers - AIAA Guidance, Navigation, and Control Conference 2006, 8(August):5342–5361, 2006. [10] Anthony P (Tony) Taylor, Ricardo (Koki) Machin, Paul Royall, and Robert Sinclair. Developing the Parachute System for NASA’s Orion: An Overview at Inception. In 19th AIAA Aerodynamic Decelerator Systems Technology Conference and Seminar, pages 772–786, Williamsburg, Virigina, may 2007. American Institute of Aeronautics and Astronautics. [11] C W Peterson, J H Strickland, and H Higuchi. The fluid dynamics of parachute inflation. Annual Review of Fluid Mechanics, 28(1):361–387, jan 1996. [12] Julian D. Maynard. Aerodynamic Characteristics of Parachutes at Mach Numbers from 1.6 to 3. Technical report, National Aeronautics and Space Administration, 1961. 123 [13] Kristin Bledsoe, Megan Englert, Aaron Morris, and Randy Olmstead. Overview of the Crew Exploration Vehicle Parachute Assembly System (CPAS) Generation I Main and Cluster Development Test Results. In 20th AIAA Aerodynamic Decelerator Systems Technology Conference and Seminar, pages 1–13, Seattle, Washington, may 2009. American Institute of Aeronautics and Astronautics. [14] K. Karagiozis, R. Kamakoti, F. Cirak, and C. Pantano. A computational study of supersonic disk-gapband parachutes using Large-Eddy Simulation coupled to a structural membrane. Journal of Fluids and Structures, 27(2):175–192, 2 2011. [15] Wei-Xi Huang, Soo Jai Shin, and Hyung Jin Sung. Simulation of flexible filaments in a uniform flow by the immersed boundary method. Journal of Computational Physics, 226(2):2206–2228, oct 2007. [16] Kenji Takizawa, Tayfun E. Tezduyar, and Ryan Kolesar. FSI modeling of the orion spacecraft drogue parachutes. Computational Mechanics, 55(6):1167–1179, 2015. [17] Charles S. Peskin. Flow patterns around heart valves: A numerical method. Journal of Computational Physics, 10(2):252–271, oct 1972. [18] Charles S. Peskin. Numerical analysis of blood flow in the heart. Journal of Computational Physics, 25(3):220–252, 1977. [19] D. Goldstein, R. Handler, and L. Sirovich. Modeling a No-Slip Flow Boundary with an External Force Field. Journal of Computational Physics, 105(2):354–366, 4 1993. [20] J. Mohd-Yusof. Combined immersed boundary/B-spline methods for simulation of flow in complex geometries. Technical report, Center for Turbulence Research, NASA Ames/Stanford Univ, 1997. [21] E. A. Fadlun, R. Verzicco, P. Orlandi, and J. Mohd-Yusof. Combined Immersed-Boundary FiniteDifference Methods for Three-Dimensional Complex Flow Simulations. Journal of Computational Physics, 161(1):35–60, 2000. [22] Markus Uhlmann. An immersed boundary method with direct forcing for the simulation of particulate flows. Journal of Computational Physics, 209(2):448–476, 2005. [23] R. Mittal and G. Iaccarino. Immersed boundary methods. Annual Review of Fluid Mechanics, 37:239–261, 2005. [24] Fotis Sotiropoulos and Xiaolei Yang. Immersed boundary methods for simulating fluid–structure interaction. Progress in Aerospace Sciences, 65:1–21, February 2014. [25] Wei-Xi Huang and Fang-Bao Tian. Recent trends and progress in the immersed boundary method. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, 233(23):7617–7636, 2019. [26] Roberto Verzicco. Immersed Boundary Methods: Historical Perspective and Future Outlook. Annual Review of Fluid Mechanics, 55(1):129–155, January 2023. [27] Rajat Mittal and Jung Hee Seo. Origin and evolution of immersed boundary methods in computational fluid dynamics. Physical Review Fluids, 8(10):100501, October 2023. [28] Kunihiko Taira and Tim Colonius. The immersed boundary method: A projection approach. Journal of Computational Physics, 225(2):2118–2137, 8 2007. 124 [29] Ugis L ˇ acis, Kunihiko Taira, and Shervin Bagheri. A stable fluid–structure-interaction solver for low- ¯ density rigid bodies using the immersed boundary projection method. Journal of Computational Physics, 305:300–318, January 2016. [30] Andres Goza and Tim Colonius. A strongly-coupled immersed-boundary formulation for thin elastic structures. Journal of Computational Physics, 336:401–411, 2017. [31] Andres Goza, Tim Colonius, and John E. Sader. Global modes and nonlinear analysis of inverted-flag flapping. Journal of Fluid Mechanics, 857:312–344, December 2018. [32] Daniel Z. Huang, Philip Avery, Charbel Farhat, Jason Rabinovitch, Armen Derkevorkian, and Lee D. Peterson. Modeling, Simulation and Validation of Supersonic Parachute Inflation Dynamics during Mars Landing. In AIAA Scitech 2020 Forum, volume 1 Part F, Orlando, Florida, jan 2020. American Institute of Aeronautics and Astronautics. [33] D. Goldstein, R. Handler, and L. Sirovich. Direct numerical simulation of turbulent flow over a modeled riblet covered surface. Journal of Fluid Mechanics, 302:333–376, 11 1995. [34] R. Glowinski, T. W. Pan, T. I. Hesla, D. D. Joseph, and J. Periaux. A Fictitious Domain Approach to ´ the Direct Numerical Simulation of Incompressible Viscous Flow past Moving Rigid Bodies: Application to Particulate Flow. Journal of Computational Physics, 169(2):363–426, 2001. [35] D. V. Le, B. C. Khoo, and K. M. Lim. An implicit-forcing immersed boundary method for simulating viscous flows in irregular domains. Computer Methods in Applied Mechanics and Engineering, 197(25-28):2119–2130, 2008. [36] Sekhar Majumdar, Gianluca Iaccarino, and Paul Durbin. RANS solvers with adaptive structured boundary non-conforming grids. Technical report, Center for Turbulence Research, 2001. [37] Yu-Heng Tseng and Joel H. Ferziger. A ghost-cell immersed boundary method for flow in complex geometry. Journal of Computational Physics, 192(2):593–623, December 2003. [38] R. Mittal, H. Dong, M. Bozkurttas, F.M. Najjar, A. Vargas, and A. Von Loebbecke. A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries. Journal of Computational Physics, 227(10):4825–4852, May 2008. [39] D. Keith Clarke, M. D. Salas, and H. A. Hassan. Euler calculations for multielement airfoils using Cartesian grids. AIAA Journal, 24(3):353–358, March 1986. [40] H.S. Udaykumar, R. Mittal, P. Rampunggoon, and A. Khanna. A Sharp Interface Cartesian Grid Method for Simulating Flows with Complex Moving Boundaries. Journal of Computational Physics, 174(1):345–380, November 2001. [41] Anna Karin Tornberg and Bjorn Engquist. Numerical approximations of singular source terms in ¨ differential equations. Journal of Computational Physics, 200(2):462–488, 2004. [42] Charles S. Peskin and David M. McQueen. A three-dimensional computational method for blood flow in the heart I. Immersed elastic fibers in a viscous incompressible fluid. Journal of Computational Physics, 81(2):372–405, 1989. [43] Charles S. Peskin. The immersed boundary method. Acta Numerica, 11:479–517, 2002. [44] Antonio Posa, Marcos Vanella, and Elias Balaras. An adaptive reconstruction for Lagrangian, directforcing, immersed-boundary methods. Journal of Computational Physics, 351:422–436, 2017. 125 [45] H. Riahi, M. Meldi, J. Favier, E. Serre, and E. Goncalves. A pressure-corrected Immersed Boundary Method for the numerical simulation of compressible flows. Journal of Computational Physics, 374:361–383, 2018. [46] R. Glowinski, T. W. Pan, and J. Periaux. Distributed Lagrange multiplier methods for incompressible ´ viscous flow around moving rigid bodies. Computer Methods in Applied Mechanics and Engineering, 151(1-2):181–194, 1998. [47] Uri M. Ascher and Linda R. Petzold. Projected Implicit Runge-Kutta Methods for DifferentialAlgebraic Equations. SIAM Journal on Numerical Analysis, 28(4):1097–1120, 8 1991. [48] Uri M. Ascher and Linda R. Petzold. Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations. SIAM, 1998. [49] Albert E. Honein and Parviz Moin. Higher entropy conservation and numerical stability of compressible turbulence simulations. Journal of Computational Physics, 201(2):531–545, 12 2004. [50] C. Pantano, R. Deiterding, D.J. Hill, and D.I. Pullin. A low numerical dissipation patch-based adaptive mesh refinement method for large-eddy simulation of compressible flows. Journal of Computational Physics, 221(1):63–87, 1 2007. [51] Ernst Hairer, Michel Roche, and Christian Lubich. The Numerical Solution of Differential-Algebraic Systems by Runge-Kutta Methods, volume 1409 of Lecture Notes in Mathematics. Springer Berlin Heidelberg, Berlin, Heidelberg, 1989. [52] Y. Kim and C.S. Peskin. 2-D parachute simulation by the immersed boundary method. SIAM J. Sci. Comput., 28:2294–2313, 2006. [53] J.M. Stockie. Modelling and simulation of porous immersed boundaries. Comput. Struct., 87(11- 12):701–709, 2009. [54] X. Gong, Z. Gong, and H. Huang. An immersed boundary method for mass transfer across permeable moving interfaces. Journal of Computational Physics, 278:148–168, 2014. [55] S. Miyauchi, S. Takeuchi, and T. Kajishima. A numerical method for mass transfer by a thin moving membrane with selective permeabilities. Journal of Computational Physics, 284:490–504, 2015. [56] S. Miyauchi, S. Takeuchi, and T. Kajishima. A numerical method for interaction problems between fluid and membranes with arbitrary permeability for fluid. Journal of Computational Physics, 345:33– 57, 2017. [57] X. Wang, X. Gong, K. Sugiyama, S. Takagi, and H. Huang. An immersed boundary method for mass transfer through porous biomembranes under large deformations. Journal of Computational Physics, 413(109444), 2020. [58] B E Schmidt. Compressible Flow Through Porous Media with Application to Injection. Technical report, California Institute of Technology, Pasadena, CA, 2015. [59] F. Ducros, F. Laporte, T. Souleres, V. Guinot, P. Moinat, and B. Caruelle. High-Order Fluxes for Con- ` servative Skew-Symmetric-like Schemes in Structured Meshes: Application to Compressible Flows. Journal of Computational Physics, 161(1):114–139, 6 2000. 126 [60] Richard D. Hornung and Scott R. Kohn. The use of object-oriented design patterns in the SAMRAI structured amr framework. In SIAM Workshop on Object-Oriented Methods for Inter-Operable Scientific and Engineering Computing, pages 235–244, Yorktown Heights, New York, 1998. SIAM. [61] S. C. R. Dennis and Gau-Zu Chang. Numerical solutions for steady flow past a circular cylinder at Reynolds numbers up to 100. Journal of Fluid Mechanics, 42(3):471–489, 7 1970. [62] Mark N. Linnick and Hermann F. Fasel. A high-order immersed interface method for simulating unsteady incompressible flows on irregular domains. Journal of Computational Physics, 204(1):157– 192, 2005. [63] Daniel Hartmann, Matthias Meinke, and Wolfgang Schroder. A strictly conservative Cartesian cut- ¨ cell method for compressible viscous flows on adaptive grids. Computer Methods in Applied Mechanics and Engineering, 200(9-12):1038–1052, 2011. [64] Daniel Canuto and Kunihiko Taira. Two-dimensional compressible viscous flow around a circular cylinder. Journal of Fluid Mechanics, 785:349–371, 2015. [65] H. Uddin, R. M.J. Kramer, and C. Pantano. A Cartesian-based embedded geometry technique with adaptive high-order finite differences for compressible flow around complex geometries. Journal of Computational Physics, 262:379–407, 2014. [66] Ming-Chih Lai and Charles S. Peskin. An Immersed Boundary Method with Formal Second-Order Accuracy and Reduced Numerical Viscosity. Journal of Computational Physics, 160(2):705–719, 5 2000. [67] K. Karagiozis, R. Kamakoti, and C. Pantano. A low numerical dissipation immersed interface method for the compressible Navier-Stokes equations. Journal of Computational Physics, 229(3):701–727, 2010. [68] Nek Sharan, Carlos Pantano, and Daniel J. Bodony. Time-stable overset grid method for hyperbolic problems using summation-by-parts operators. Journal of Computational Physics, 361:199–230, 2018. [69] Eleuterio F. Toro. Riemann Solvers and Numerical Methods for Fluid Dynamics. Springer Berlin Heidelberg, Berlin, Heidelberg, 2009. [70] Hang Yu and Carlos Pantano. An immersed boundary method with implicit body force for compressible viscous flow. Journal of Computational Physics, 459:111125, jun 2022. [71] Andres Goza, Sebastian Liska, Benjamin Morley, and Tim Colonius. Accurate computation of surface stresses and forces with immersed boundary methods. Journal of Computational Physics, 321:860– 873, September 2016. [72] Bakytzhan Kallemov, Amneet Bhalla, Boyce Griffith, and Aleksandar Donev. An immersed boundary method for rigid bodies. Communications in Applied Mathematics and Computational Science, 11(1):79–141, 2016. [73] Per Christian Hansen. Discrete Inverse Problems: Insight and Algorithms. Society for Industrial and Applied Mathematics, 2010. [74] J. E. Marsden and L. Sirovich, editors. Linear Integral Equations, volume 82 of Applied Mathematical Sciences. Springer New York, New York, NY, 1999. 127 [75] Heinz W. Engl, Martin Hanke, and Andreas Neubauer. Regularization of Inverse Problems. Springer Netherlands, Dordrecht, 1996. [76] Rainer Kress. Numerical Analysis, volume 181 of Graduate Texts in Mathematics. Springer New York, New York, NY, 1998. [77] Turker Biyiko ¨ glu, Josef Leydold, and Peter F. Stadler. ˇ Laplacian Eigenvectors of Graphs, volume 1915. Springer Berlin Heidelberg, Berlin, Heidelberg, 2007. [78] Silvia Noschese and Lothar Reichel. Lavrentiev-type regularization methods for Hermitian problems. Calcolo, 52(2):187–205, 2015. [79] Stefan Turek and Jaroslav Hron. Proposal for numerical benchmarking of fluid-structure interaction between an elastic object and laminar incompressible flow. In Hans-Joachim Bungartz and Michael Schafer, editors, ¨ Fluid-Structure Interaction, pages 371–385, Berlin, Heidelberg, 2006. Springer Berlin Heidelberg. [80] Rajneesh Bhardwaj and Rajat Mittal. Benchmarking a Coupled Immersed-Boundary-Finite-Element Solver for Large-Scale Flow-Induced Deformation. AIAA Journal, 50(7):1638–1642, July 2012. [81] Fang-Bao Tian, Hu Dai, Haoxiang Luo, James F. Doyle, and Bernard Rousseau. Fluid–structure interaction involving large deformations: 3D simulations and applications to biological systems. Journal of Computational Physics, 258:451–469, February 2014. [82] Sebastian Geller, Jonas Tolke, and Manfred Krafczyk. Lattice-boltzmann method on quadtree-type ¨ grids for fluid-structure interaction. In Hans-Joachim Bungartz and Michael Schafer, editors, ¨ FluidStructure Interaction, pages 270–293, Berlin, Heidelberg, 2006. Springer Berlin Heidelberg. [83] Michael Schafer, Marcus Heck, and Saim Yigit. An implicit partitioned method for the numerical ¨ simulation of fluid-structure interaction. In Hans-Joachim Bungartz and Michael Schafer, editors, ¨ Fluid-Structure Interaction, pages 171–194, Berlin, Heidelberg, 2006. Springer Berlin Heidelberg. [84] Hang Yu and Carlos Pantano. A regularized projection immersed boundary method for smooth boundary forces. Journal of Computational Physics, 496:112571, 2024. [85] Boyce E. Griffith and Neelesh A. Patankar. Immersed Methods for Fluid–Structure Interaction. Annual Review of Fluid Mechanics, 52(1):421–448, January 2020. [86] Marsha J Berger and Joseph Oliger. Adaptive mesh refinement for hyperbolic partial differential equations. Journal of Computational Physics, 53(3):484–512, 1984. [87] M.J. Berger and P. Colella. Local adaptive mesh refinement for shock hydrodynamics. Journal of Computational Physics, 82(1):64–84, 1989. [88] Fehmi Cirak, Michael Ortiz, and Peter Schroder. Subdivision surfaces: a new paradigm for thin-shell ¨ finite-element analysis. International Journal for Numerical Methods in Engineering, 47(12):2039– 2072, 4 2000. [89] Fehmi Cirak and Michael Ortiz. Fully C1-conforming subdivision elements for finite deformation thin-shell analysis. International Journal for Numerical Methods in Engineering, 51(7):813–833, 2001. [90] Qingquan Wang and Carlos Pantano. A parallel geometric contact algorithm for thin shell finite elements in explicit time integration. Computers and Structures, under review. 128 [91] P. Causin, J.F. Gerbeau, and F. Nobile. Added-mass effect in the design of partitioned algorithms for fluid–structure problems. Computer Methods in Applied Mechanics and Engineering, 194(42):4506– 4527, 2005. [92] Iman Borazjani, Liang Ge, and Fotis Sotiropoulos. Curvilinear immersed boundary method for simulating fluid structure interaction with complex 3D rigid bodies. Journal of Computational Physics, 227(16):7587–7620, August 2008. [93] Andrew M. Wissink, Richard D. Hornung, Scott R. Kohn, Steve S. Smith, and Noah Elliott. Large scale parallel structured AMR calculations using the SAMRAI framework. In Proceedings of the 2001 ACM/IEEE conference on Supercomputing, pages 6–6. ACM, 2001. [94] M. Berger and I. Rigoutsos. An algorithm for point clustering and grid generation. IEEE Transactions on Systems, Man, and Cybernetics, 21(5):1278–1286, October 1991. [95] M. TABATA and K. ITAKURA. A precise computation of drag coefficients of a sphere. International Journal of Computational Fluid Dynamics, 9(3):303–311, 1998. [96] Susumu Shirayama. Flow past a sphere - topological transitions of the vorticity field. AIAA Journal, 30(2):349–358, 1992. [97] Rubens Campregher, Julio Militzer, Sergio Said Mansur, and Aristeu Da Silveira Neto. Computations ´ of the flow past a still sphere at moderate reynolds numbers using an immersed boundary method. Journal of the Brazilian Society of Mechanical Sciences and Engineering, 31(4), 2024. [98] Keith Stein, Richard Benney, Vinay Kalro, Tayfun E. Tezduyar, John Leonard, and Michael Accorsi. Parachute fluid-structure interactions: 3-D computation. Computer Methods in Applied Mechanics and Engineering, 190(3-4):373–386, 2000. [99] Michael Barnhardt, Travis Drayna, Ioannis Nompelis, Graham Candler, and William Garrard. Detached Eddy Simulations of the MSL Parachute at Supersonic Conditions. In 19th AIAA Aerodynamic Decelerator Systems Technology Conference and Seminar, pages 263–273, Williamsburg, Virigina, may 2007. American Institute of Aeronautics and Astronautics. [100] John Lingard, Matt Darley, John Underwood, and Glen Brown. Simulation of the Mars Science Laboratory Parachute Performance and Dynamics. In 19th AIAA Aerodynamic Decelerator Systems Technology Conference and Seminar, pages 1–12, Williamsburg, Virigina, may 2007. American Institute of Aeronautics and Astronautics. [101] Yongsam Kim and Charles S. Peskin. 3-D Parachute simulation by the immersed boundary method. Computers and Fluids, 38(6):1080–1090, 2009. [102] Jason Rabinovitch, Faisal As’ad, Philip Avery, Charbel Farhat, Navid Ataei, and Marcus Lobbia. Update: Modeling Supersonic Parachute Inflations for Mars Spacecraft. In 26th AIAA Aerodynamic Decelerator Systems Technology Conference, Toulouse, FRANCE, May 2022. American Institute of Aeronautics and Astronautics. 129 Appendices A Derivative of discontinuous function and the surface singularity A Dirac delta appears when differentiating a function with a jump in the sense of distributions. When differentiating a function on higher dimension space R d with discontinuity alongside a surface Γ(t), surface singularity on Γ(t) appears. Herein, let Γ(t) be a differentiable surface parametrized by s with function xbc(s,t). Let function f(x,t) be the function with a discontinuity on Γ(t) and continuously differential elsewhere in R d \ Γ(t). Taking the derivative ∂ f /∂ xk as distribution and applying test function ϕ(x), we have Z Rd ∂ f ∂ xk ϕ(x)dx ≡ −Z Rd f(x,t) ∂ ϕ(x) ∂ xk dx = − Z Rd \Γε (t) f(x,t) ∂ ϕ(x) ∂ xk dx− Z Γε (t) f(x,t) ∂ ϕ(x) ∂ xk dx = Z Rd \Γε (t) ∂ f ∂ xk ϕ(x)dx | {z } Iε + I ∂Γε (t) f(x,t)ϕ(x)nˆε,kds | {z } IIε − Z Γε (t) f(x,t) ∂ ϕ(x) ∂ xk dx | {z } IIIε , (A.1) where Γε (t) = n y ∈ R d : dist(y,Γ) < ε o , (A.2) and nˆ ε = [nˆε,1,··· ,nˆε,d] is the unit outward normal on ∂Γε (t). We take the limit ε → 0 for the first term: lim ε→0 Z Rd \Γε (t) ∂ f ∂ xk ϕ(x)dx = Z Rd \Γ(t) ∂ f ∂ xk ϕ(x)dx, (A.3) 130 provided that ∂ f /∂ xk is bounded. For the second term, we have lim ε→0 I ∂Γε (t) f(x,t)ϕ(x)nˆε,kds = Z Γ(t) [ f +(xbc(s,t),t)− f −(xbc(s,t),t)]ϕ(xbc(s,t))nˆkds = Z Γ(t) J f K(s,t)ϕ(xbc(s,t),t)nˆkds. (A.4) The third term is simply zero when taking ε → 0 since the integrand is bounded. Therefore, combining the results we have Z Rd ∂ f ∂ xk ϕ(x)dx = Z Rd \Γ(t) ∂ f ∂ xk ϕ(x)dx+ Z Γ(t) J f K(s,t)ϕ(xbc(s,t),t)nˆkds, (A.5) for all test function ϕ(x) of R d , implying that the surface singularity that arises from ∂ f /∂ xk is Z Γ(t) J f K(s,t)nˆk(s,t)δ(x−xbc(s,t))ds. (A.6) This surface singularity is also denoted as δ(Γ(t), J f Knk,x) in [41]. Similarly, we differentiate f in time and ∂ f /∂t has surface singularity when the discontinuous surface Γ(t) is moving. Applying ∂ f /∂t with test function ϕ(x) gives Z Rd ∂ f ∂t ϕ(x)dx = Z Rd \Γε (t) ∂ f ∂t ϕ(x)dx+ Z Γε (t) ∂ f ∂t ϕ(x)dx = Z Rd \Γε (t) ∂ f ∂t ϕ(x)dx | {z } Iε + d dt Z Γε (t) f(x,t)ϕ(x)dx | {z } IIε − I ∂Γε (t) f(x,t)ϕ(x) (x˙ ε ·nˆ ε )ds | {z } IIIε , (A.7) where x˙ ε is the velocity of boundary ∂Γε (t). Now we take the limit ε → 0 for the first term: lim ε→0 Z Rd \Γε (t) ∂ f ∂t ϕ(x)dx = Z Rd \Γ(t) ∂ f ∂t ϕ(x)dx. (A.8) 131 While for the second term we have lim ε→0 d dt Z Γε (t) f(x,t)ϕ(x)dx = d dt lim ε→0 Z Γε (t) f(x,t)ϕ(x)dx = 0, (A.9) where the swap of order is possible if | f(x,t)| is bounded by a function of x for all t, and the integral goes zero when ε → 0 because of the boundedness of the integrand. Finally, for the third term, we obtain lim ε→0 I ∂Γε (t) f(x,t)ϕ(x) (x˙ ε ·nˆ ε )ds = Z Γ(t) J f K(s,t)ϕ(xbc(s,t)) (x˙bc ·nˆ)ds. (A.10) where x˙bc(s,t) = ∂xbc ∂t is the velocity of the surface. Combining the three terms, we end up with Z Rd ∂ f ∂t ϕ(x)dx = Z Rd \Γ(t) ∂ f ∂t ϕ(x)dx− Z Γ(t) J f K(s,t)ϕ(xbc(s,t)) (x˙bc ·nˆ)ds, (A.11) for all test functions. Therefore, the surface singular term arising from the temporal derivative is − Z Γ(t) J f K(s,t)(x˙bc(s,t)·nˆ(s,t))δ(x−xbc(s,t))ds. (A.12) 132
Abstract (if available)
Abstract
We present a projection immersed boundary (IB) method for simulating compressible viscous flow. The immersed boundary is represented with arbitrary surface grids, and the fluid dynamic is resolved with Cartesian grids. Convolutions with discrete delta functions are used to transfer data between the two grids. Singular source terms are added to the governing equations of the fluid to enforce the desired boundary conditions. The resultant equations are a system of differential-algebraic equations (DAE). The DAE system is solved with a half-explicit Runge-Kutta method, and each stage is solved with a fractional step (projection) method.The projection IB method, however, fails to produce accurate boundary forces that converge to the actual physical values when the grid is refined. This happens because the singular source terms are solutions of discrete Fredholm integral equations, which are not well-posed because the discrete integral operators do not have a bounded inverse as the grid refines. To alleviate this issue, we proposed a generalized Tikhonov regularization approach that reformulates the problem into a minimization problem. The graph Laplacian is used to formulate the seminorm measuring the roughness of the solution.The method developed is implemented with a parallelized patch-based adaptive mesh refinement (AMR) framework. The fluid dynamic solver with IB and AMR is coupled with a thin-shell finite element Kirchhoff-Love plate solver as a loosely coupled fluid-structure interaction (FSI) solver. In the end, the FSI solver is used to simulate the drogue parachute cluster of the Orion Capsule Assembly System (CPAS). Preliminary results will be reported.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Numerical study of shock-wave/turbulent boundary layer interactions on flexible and rigid panels with wall-modeled large-eddy simulations
PDF
Supersonic flow study of a capsule geometry using large-eddy simulation
PDF
Hydrodynamic stability of two fluid columns of different densities and viscosities subject to gravity
PDF
Tracking and evolution of compressible turbulent flow structures
PDF
Dynamic modeling and simulation of flapping-wing micro air vehicles
PDF
Direct numerical simulation of mixing and combustion under canonical shock turbulence interaction
PDF
A high-fidelity geometric contact method for thin fabrics using shell finite elements in explicit and implicit time scheme
PDF
Passive rolling and flapping dynamics of a heaving Λ flyer
PDF
Mach limits and free boundary problems in fluid dynamics
PDF
Multiscale modeling of cilia mechanics and function
PDF
Pulsatile and steady flow in central veins and its impact on right heart-brain hemodynamic coupling in health and disease
PDF
Development and applications of a body-force propulsor model for high-fidelity CFD
PDF
Modeling and analysis of parallel and spatially-evolving wall-bounded shear flows
PDF
Computational fluid dynamic analysis of highway bridge superstructures exposed to hurricane waves
PDF
On regularity and stability in fluid dynamics
PDF
Oscillatory and streaming flow due to small-amplitude vibrations in spherical geometry
PDF
Large eddy simulations of turbulent flows without use of the eddy viscosity concept
PDF
Development and evaluation of high‐order methods in computational fluid dynamics
PDF
Scattering of a plane harmonic SH-wave and dynamic stress concentration for multiple multilayered inclusions embedded in an elastic half-space
PDF
Grid-based Vlasov method for kinetic plasma simulations
Asset Metadata
Creator
Yu, Hang
(author)
Core Title
The projection immersed boundary method for compressible flow and its application in fluid-structure interaction simulations of parachute systems
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Mechanical Engineering
Degree Conferral Date
2024-12
Publication Date
09/26/2024
Defense Date
08/23/2024
Publisher
Los Angeles, California
(original),
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
computational fluid dynamics,fluid-structure interaction,immersed boundary methods,OAI-PMH Harvest,Parachutes
Format
theses
(aat)
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Pantano-Rubino, Carlos (
committee chair
), Bermejo-Moreno, Ivan (
committee member
), Ghanem, Roger (
committee member
), Newton, Paul (
committee member
)
Creator Email
yuhang@usc.edu,yuhang94@outlook.com
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-oUC11399B9ZU
Unique identifier
UC11399B9ZU
Identifier
etd-YuHang-13547.pdf (filename)
Legacy Identifier
etd-YuHang-13547
Document Type
Dissertation
Format
theses (aat)
Rights
Yu, Hang
Internet Media Type
application/pdf
Type
texts
Source
20240926-usctheses-batch-1214
(batch),
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the author, as the original true and official version of the work, but does not grant the reader permission to use the work if the desired use is covered by copyright. It is the author, as rights holder, who must provide use permission if such use is covered by copyright.
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA
Repository Email
cisadmin@lib.usc.edu
Tags
computational fluid dynamics
fluid-structure interaction
immersed boundary methods