Eigensolution analysis of immersed boundary method based on volume penalization: applications to high-order schemes

This paper presents eigensolution and non-modal analyses for immersed boundary methods (IBMs) based on volume penalization for the linear advection equation. This approach is used to analyze the behavior of flux reconstruction (FR) discretization, including the influence of polynomial order and penalization parameter on numerical errors and stability. Through a semi-discrete analysis, we find that the inclusion of IBM adds additional dissipation without changing significantly the dispersion of the original numerical discretization. This agrees with the physical intuition that in this type of approach, the solid wall is modelled as a porous medium with vanishing viscosity. From a stability point view, the selection of penalty parameter can be analyzed based on a fully-discrete analysis, which leads to practical guidelines on the selection of penalization parameter. Numerical experiments indicate that the penalization term needs to be increased to damp oscillations inside the solid (i.e. porous region), which leads to undesirable time step restrictions. As an alternative, we propose to include a second-order term in the solid for the no-slip wall boundary condition. Results show that by adding a second-order term we improve the overall accuracy with relaxed time step restriction. This indicates that the optimal value of the penalization parameter and the second-order damping can be carefully chosen to obtain a more accurate scheme. Finally, the approximated relationship between these two parameters is obtained and used as a guideline to select the optimum penalty terms in a Navier-Stokes solver, to simulate flow past a cylinder.


Introduction
Dispersion and dissipation behavior is of paramount importance for all numerical schemes that solve Partial Differential Equations (PDEs). This behavior can be characterised using eigensolution analysis, which quantifies how the amplitude and frequency of a wave-like solution evolve with time, providing useful information on the numerical errors. Eigensolution analysis, also known as von Neumann or Fourier analysis, has been widely applied to different numerical schemes for space discretization, including finite difference [1,2], finite volume [3], finite element [4], Lattice Boltzmann [5,6], as well as high-order methods [7,8,9,10,11,12,13]. Typically, eigendecomposition is applied to the global discretization matrix, where dispersion and dissipation characteristics reflect in the corresponding eigenvalues. The dispersion error, which is the error in wave advection, is represented by the modified frequencies of the solution. The dissipation error, which corresponds to the nonphysical damping or amplification effect, is represented by the modified amplitudes of the solution. This dispersion-dissipation behavior is crucial to understand and quantity the numerical error of a numerical scheme, and can be used to evaluate its stability and robustness. In addition, insights from the analysis can be extended to design better numerical schemes for turbulence simulation [14,15,16]. Eigensolution analysis in multi-dimensions can also be used to evaluate the effect of mesh quality for high-order schemes [17]. Most analyses belong to temporal eigensolution analysis, since the focus of such analysis is on the temporal evolution of the solution by considering periodic boundary conditions. To investigate the evolution in space, spatial eigensolution analysis [18,19] has been proposed recently for high-order schemes with inflow-outflow boundary conditions. Both spatial and temporal eigensolution analyses have been extended in a non-intrusive and data-driven manner to evaluate dispersion-dissipation behavior only from simulation data [20]. Another alternative perspective for the dissipation error is to look at the non-modal behavior characterized by the short-term dissipation of a numerical scheme [21], which has been shown to obtain useful insights for the hybridized Discontinuous Galerkin (DG) scheme. When the time integration scheme is incorporated into the analysis, a fully-discrete analysis [22,23] is performed and takes into account the numerical error both in space and time.
So far, eigensolution analyses are still limited to body-fitted grids. The immersed boundary method (IBM) [24,25,26] provides an alternative to body-fitted meshes. This technique has gained attention during recent years and has the potential to handle complex geometry and moving bodies on simple Cartesian grids, thus reducing the cost in mesh generation. IBM originates from the work of Peskin [27] to introduce a singular source term to the background flow grid in the vicinity of the solid body. So far, IBM been extensively studied in various applications including benchmark geometries [28,29], turbulent flows [30,31], fluid-structure interaction (FSI) [32,33,34], acoustics [35,36], multiphase flows [37,38], etc. To understand the numerical behaviors and better implement IBM treatment in practical numerical simulation, dispersion-dissipation analysis can be adapted to IBM, which is a subject that has been deeply explored in the present study.
In general, IBM can be achieved in multiple ways, which are broadly divided into two strategies, including cutting the cell by the solid boundary [39,40] or introducing additional forcing terms [41,42,43,44] to mimic the effect of solid objects. As a typical approach of the second kind, the volume penalization method [42,45,46,47] has been widely adopted for boundary treatment due to its good robustness, simplicity and rigorous theoretical background. It is based on the physical intuition that the solid wall can be modelled as a porous medium with vanishing diffusivity [48]. A source or penalization term is introduced inside the solid domain, whose value depends on the diffusivity coefficient (the penalization parameter) that needs to be decided by the user. Applications of volume penalization include flapping wings [49], two-phase flow [50], FSI [51] and thermal flows [52]. Recently, volume penalization has been applied to high-order flux reconstruction to evaluate the accuracy of IBM in high-order numerical schemes [53,54]. This type of boundary treatment is the focus of this study.
There have been several attempts to investigate the the numerical properties of IBM. For classical IBM where the singular force is described by delta functions, linear stability analysis is performed to study the stability and stiffness of IBM coupled with time-stepping methods [55,56,57]. It was found that the explicit time stepping scheme has the worst stability and can be unstable even though the underlying physical system is stable due to the stiffness of the IBM source term. Spectral analysis for Laplace and Stokes operators based on volume penalization and FD discretization has been studied in [58,59] for either Dirichlet or Neumann boundary conditions, where it is suggested that if one wants to reach higher accuracy, the underlying discretization should behave well in the presence of discontinuous solutions. Linear stability analysis for the second-order Adams-Bashforth time integration of the penalization term in the incompressible Navier-Stokes equations was performed in [49], where the stability condition ∆t < η was proposed. The eigenvalue problem for elliptic problems based on immersed finite element method is studied in [60]. Recently, eigenvalue and error analysis of the immersed boundary method based on direct forcing is studied by Zhou and Balachandar [61]. All of these efforts provide valuable insights to deepen the understanding of IBM, while some remaining issues still exist. For example, the IBM treatment based on high-order discretization has not been discussed, as well as some guidelines to improve the implementation in these settings.
In addition to commonly used low-order methods based on FD and FV, there is an emerging interest in developing high-order methods for computational fluid dynamics (CFD) due to their potential in providing higher accuracy with relatively low cost compared to low-order methods [62]. Currently, different high-order methods have been developed, including DG [63,64], spectral difference (SD) [65,66], flux reconstruction (FR) [67], and correction procedure via reconstruction (CPR) [68]. Despite vast efforts, applying high-order methods on unstructured grids remains a challenge due to the difficulty of mesh generation for complex geometries. This makes the development of IBM under high-order frameworks a very appealing alternative. As a representative high-order method, FR is chosen in the present study, which unifies nodal DG and SD schemes in certain conditions [69] (e.g. pure advection, as considered in this work). Therefore our analysis covers both FR and DG formulations. This scheme will be combined with volume penalization to simulate the advection equation with a solid wall in the middle of computational domain, thus providing insights into the numerical behavior of IBM on high-order schemes.
The remainder of this paper is organized as follows. In Section 2, the FR method for space discretization and volume penalization method are introduced. Next, Section 3 introduces the eigensolution and non-modal analyses, as well as the novel analysis approach considering the IBM treatment. Additionally, both semi-discrete and fullydiscrete analysis are considered to study the behaviors of various factors, including the effects of polynomial order and selection of penalization parameter. Furthermore, numerical experiments are performed in Section 4 to validate and better understand the results from the analysis. A novel approach to include artificial viscosity (second-order term) inside the solid is also proposed and evaluated. Finally, conclusions are drawn in Section 5.

The governing equation
We consider a one-dimensional advection equation defined in space x and time t: where u(x, t) is the transported unknown solution and f = cu is the flux, with advection speed c. To apply eigensolution analysis to this equation, we first obtain the analytical solution for a harmonic wave under the periodic boundary conditions and the initial condition u(x, 0) = exp(ikx): where k is the wavenumber, ω = ck is the (angular) frequency and i = √ −1 is the imaginary unit. In the present work, this equation is discretized in space based on the high-order flux reconstruction method [67,69] (which is equivalent to DG). Details in space and time discretization of the advection equation are given in Appendix Appendix A. In classic eigensolution analyses, one only needs to look at one element rather than the global domain since it represents the global behavior of the system under periodic or inflow-outflow boundary conditions. However, when IBM is considered, the global matrices for all elements will be taken into account, since some of the solution points should be penalized by additional source terms to impose the IBM conditions. This will be detailed in the following sections.

Immersed boundary method based on volume penalization
Volume penalization imposes boundary conditions by introducing source terms inside the solid region, which is assumed to be a porous medium whose permeability tends to zero. As shown in previous studies [42,49,48], this method is easy to implement with rigorous theoretical foundations. A mask function χ(x, t) that distinguishes between the fluid region Ω f and solid region Ω s is defined first: For the one-dimensional advection equation with penalized Dirichlet boundary conditions, this problem is formulated as follows: where u s refers to the penalized solution in the solid, which is fixed to 0 in this study to approximate the no-slip wall boundary condition. The additional source term is usually a pointwise operation imposed to each solution point. The extension to Neumann boundary conditions is discussed in [48,70]. The penalization parameter is defined as η, which drives the solution to u s as η → 0. For solution inside the solid, this results in the following governing equation similar with convection-diffusion-reaction equations [71] ∂u ∂t + c ∂u ∂x The analytical solution of this equation is indicating that the penalization term introduces additional damping to the solution, controlled by the penalization parameter η. This indicates that volume penalization mimics the porous media with a very low permeability which dissipate the wave-like solution entering the solid and drive it to the expected solution u s . From the following eigensolution analysis, this penalization term will have an impact on both the primary and secondary modes by proving additional dissipation.
There are a few works that have discussed the convergence and error estimation of volume penalization, where rigorous proofs have been given [42,72]. From Angot et al. [42] and Carbou and Fabrie [72], it was proven that, as the penalization parameter η approaches 0, the solution of the penalized Navier-Stokes equations will converge to the solution of the Navier-Stokes equations with no-slip boundary conditions. This is one of the advantages of volume penalization over other IBM approaches since the numerical error introduced from the penalization term can be controlled a-priori [45]. The total error of the solution of penalized equation compared with the body-fitted simulation (non-penalized equation) includes two parts [51] where the first part is the penalization error and the second part is the discretization error. u exact , u η and u N η are the exact analytical solution of the governing equations, the exact and numerical solution of the penalized equations, respectively. The error is quantified by the norm · between solutions. The penalization error depends on the penalization parameter [47]: This indicates that as the penalization parameter approaches zero, the error between exact solutions of penalized and original equation will converge to zero, i.e., lim η→0 u exact − u exact η → 0. It has been proved that the volume penalization gives α = 1 2 for the Dirichlet boundary conditions, indicating the penalization error has a decay rate of O( √ η). For Neumann boundary condition, O(η) can be obtained [48,59]. The discretization error refers to the error between the exact solution and the numerical solution of the penalized equations. With consistent discretization and a stable numerical scheme, the discretization error usually follows (β > 0): However, as pointed out by Schneider et al. [48,47], the convergence for the discretization error is not only determined by the numerical scheme, but also limited by the regularity of the solution, which refers to the continuity (smoothness) of the exact solution u η at the boundary of the penalized equation. This regularity can be improved by designing proper u s inside the solid [73,74]. However, in this study we will discuss the classical volume penalization method with u s = 0 for a no-slip wall. This analysis also indicates the difficulty of IBM to achieve high order convergence near the boundary, with any type of space discretization. From these theories, it is suggested to use a penalization parameter η which is small enough to ensure low penalization error (high penalization since η is in the denominator). In the meantime, η should also be limited since the penalization source term can become very stiff, leading to numerical stabilities. When explicit time integration is considered, it is suggested to use ∆t < η [49] or ∆t ≈ η [75] to ensure stability. An interpretation for this choice is that the penalty term acts as a strong damping term with order η on the velocity, which has to be resolved by the time discretization scheme [51]. We will provide further insights into the selection of the penalty parameter in what follows.

Eigensolution analysis
In the present temporal eigensolution analysis, we assume periodic boundary condition and consider solution u(x, t) = exp[i(kx − ωt)] with real wavenumber k. Therefore, for each element with index n, we have u n−1 = e −ikh u n and u n+1 = e +ikh u n . The space discretization results in du n dt = Le −ikh + C + Re +ikh u n = Mu n .
This semi-discrete matrix formulation is the basis for the eigensolution analysis and the non-modal analysis. It can be transformed into the compact form of an eigenvalue problem where ω * = ck * becomes a complex value due to the dispersion and dissipation errors of the space discretization. This wavenumber ω * relates to the eigendecomposition of the coefficient matrix M (with P + 1 solutions) where λ m and v m are the mth eigenvalues and eigenvectors of matrix M, respectively. By defining the modified wavenumber k * m for each ω * m and λ m , we have: For the advection equation, the difference between Real(k * ) and k is due to the dispersion error, indicating the change in wavenumber of the solution. The difference between Imag(k * ) and 0 corresponds to the dissipation (diffusion) error, where Imag(k * ) ≤ 0 holds for a stable scheme. Note that for the analytical solution we have Real(k * ) = k and Imag(k * ) = 0. For high-order methods, the wavenumber is normalized by the h/(P + 1) which represents the smallest length scale that can be captured by the scheme. In addition, M will have P + 1 eigenmodes, where the one that recovers k * = k as k → 0 will be identified as the physical mode, while others are secondary modes which are indeed translated replicas of the primary mode [11].
When the IBM is included in the analysis, one has to consider the whole computational domain rather than one element, while penalizing the solution points inside the solid domain. To account for multiple (equispaced) elements, we consider the direct sum of N-element contributions as shown in [13], with solution vector defined as U = [u 1 , u 2 , ..., u N ] T . We focus on the computational domain [−T, T ] with periodic boundary condition for the advection equation. In this case, the effect of the physical boundary conditions need to be introduced through two ghost elements for external states, u 0 and u N+1 , derived as u 0 = exp[−2ikT ]u N and u N+1 = exp[2ikT ]u 1 . For cases without IBM, this represents the continuation of a hypothetical infinite mesh. Since volume penalization imposes pointwise source terms, the penalization term is only added to the diagonal of the global matrix. Following [13], we can reach the global matrix form of semi-discrete formulation considering IBM as follows: where I denotes the identity matrix with dimension P + 1. The volume penalization term is imposed when the corresponding solution point lies inside the solid (χ = 1). For example, assuming we have one solid element in the middle, the global matrix becomes d dt This is equivalent to an advection problem with a solid wall in the middle of a periodic computational domain. By investigating this global matrix, the dispersion-dissipation behavior of IBM for high-order schemes can be analyzed. Note that this matrix will contain N(P + 1) eigenmodes, but as pointed out in [13], only one of the modes tracks the physical propagation speed and damping, which is referred to as the primary mode or physical mode (again, the one that recovers k * = k as k → 0 [11]). All remaining modes are referred to as secondary modes.

Non-modal analysis
Non-modal analysis [21] is an approach to analyze the numerical diffusion, which is more relevant in underresolved turbulence simulations. A short-term diffusion is defined that describes the change of magnitude of the numerical solution evolving right after time t = 0. It differs from classic eigensolution analysis since the contribution of all eigenmodes is taken into account, which leads to bigger difference in the high wavenumber range. The shortterm diffusion parameterω * is defined as [21] where || · || is the L 2 norm and τ * = τ(P + 1) = tc(P + 1)/h is a non-dimensional time based on the non-dimensional convection time per degree of freedom. Forω * ≤ 0 the scheme is stable. Considering the orthogonality of solution polynomials and the wave-like behavior of the numerical solution, the above equation can be rewritten as the form of a Rayleigh quotient: where u n,0 is the initial condition, and † denotes the complex conjugate. M is the discretization operator in Equation  12 and Equation 15. The short-term diffusion parameter can be understood as the combined contribution of all eigenmodes of the discretization operator identified in a von Neumann analysis [21]. Therefore, non-modal analysis is consistent with the von Neumann analysis whenever only one eigenmode exists (at P = 0) or when the initial solution is an eigenvector of the discretization operator M.

Semi-discrete analysis
In this section, a semi-discrete eigensolution analysis of the FR scheme with volume penalization for the advection equation is performed. Since semi-discrete analysis only looks at the space discretization, it performs eigensolution analysis and non-modal analysis on the space discretization operator matrix. The analysis of the standard advection equation focuses on Equation 12, while the analysis of IBM focuses on Equation 15. In the following analysis, we consider the advection equation with c = 1 and a computational domain defined in x ∈ [−1, 1] discretized by equispaced elements. All the cases considered here are based on the fully upwind flux with λ = 1. The solid region lies in the middle, starting from x = 0, whose width is defined as ∆: In the present study, we will consider ∆ as integer multiples of element size h (∆ = Zh, where Z is an integer), thus the solid boundaries will appear exactly at the interface between elements. This allows us to define the solid ratio r = Z/N as the ratio between the solid region and the computational domain. For a wavelike initial condition, the solution is expected to approach 0 when x > 0 and x < t (consider a short period). As time evolves, the global solution will eventually become 0 in the whole computation domain, due to the solid and the periodic boundary condition. Moreover, when the solid wall exists, it is expected to have c = 0 for all solution points inside the wall. This indicates that the actual fluid domain is shorter than the whole computational domain, therefore the exact wavenumber in the fluid region is higher than the initial wavenumber. To obtain the real wavelength in the fluid region, the initial wavenumber k should be re-scaled by the solid ratio r This definition is consistent with the standard advection equation where we have r = 0 thusk = k. Note that this re-scaling only applies to the initial wavenumber k to match the modified wavenumber obtained from the eigensolution analysis k * . It should also be noted that the inclusion of IBM sometimes requires shifting the wavenumber by π/N to accurately trace the physical mode in eigensolution analysis. In summary, a schematic illustration of the considered problem is shown in Figure 1, where we consider an initial wave for a short period convection. The IBM wall is introduced in the middle, therefore the initial wave in the middle starts to be damped and will move towards right as the solution evolves. This problem is different from the assumption of classical eigensolution analysis where the solution is decomposed into the combination of harmonics. Therefore to better understand the present problem settings, we can think of eigensolution analysis as transforming the present problem into an equivalent problem in the right of Figure 1, where the solution can be represented as a harmonic wave with modified wavelengthk and equivalent amplitude (damping). The wavelength is re-scaled due to the existence of the solid region, as defined in Equation 19. Since the analysed equation is linear, the equivalent damping is an average dissipation effect coming of both the local volume penalization term (averaged over the whole computational domain) and the FR scheme itself. By looking as this equivalent problem, we have the following questions: 1) how does the penalization parameter influences the dispersion-dissipation behavior? 2) will IBM change the spectral behaviors of the original FR scheme? Both of them will be answered in the following analysis and numerical tests.
As the first example, we consider a coarse discretization, where the computational domain is discretized into N = 10 elements, with the solid width ∆ = h = 0.2 and r = 1/10. The polynomial order is set to P = 3, and the penalization parameter is set to η = 1 × 10 −4 . Results from eigensolution and non-modal analysis are shown in Figure  2, where all eigenmodes are shown. In general, the dispersion-dissipation behavior agrees well with those of the standard advection equation without IBM, as reported in [10,11,13], where one physical / primary mode is identified while other modes are mostly replications of the primary mode. One distinction of the present case is that the volume penalization term will introduce several solid modes (depending on the number of degrees of freedom inside the solid) with constant dispersion and dissipation across all wavenumbers, as shown in Figure 2a and Figure 2b. They are new secondary modes when IBM is considered, which also agrees with the analysis of Equation 6 that constant dissipation is introduced inside the wall. To understand the spectral behavior, the physical mode is more interesting. It can be seen that IBM retains reasonable dispersion behavior at low and medium wavenumbers, and introduces additional dissipation for the primary mode across the wavenumber range, as shown in Figure 2c. This is reasonable because when any solution (including the constant solution k = 0) passes through the IBM wall, it will be damped in order to satisfy the penalized solution u s . Therefore, the dissipation at k = 0 refers to the damping only induced by the IBM condition, which is the shift of a constant solution. When we look at the effect of IBM on the FR scheme, this constant dissipation needs to be subtracted from the total dissipation, as illustrated in Figure 1. Finally, the total dissipation effect for all fluid and solid modes can be seen in the short-term dissipation curve in Figure 2d, reflecting the combined damping effect in the whole computational domain. It can be concluded that the short-term dissipation is dominated by the IBM source term, resulting in very large negative dissipation depending on η. This indicates that most of the damping is due to the effect of IBM, resulting very large short-term dissipation at k = 0. Again, the effect of IBM on the FR scheme should be obtained by subtracting this initial dissipation. The magnitude of this short-term damping will be smaller than that of the single solid mode in the dissipation curve in Figure 2b, since the short-term damping is an averaged damping in the computational domain, but the solid mode is only a characteristic of the solid body. It should be mentioned that the magnitude of these damping factors rely heavily on η which determines the decay of the solution inside the solid. Following this analysis, we can define the IBM-induced dissipation γ IBM for the primary mode as where k * (0) refers to the modified wavenumber k * at k = 0. Therefore, the total dissipation γ of primary mode is a combination of IBM (constant dissipation across all wave numbers) and the numerical scheme (modified behavior due to the existence of IBM) which also applies to short-term dissipation based on the non-modal analysis. By comparing γ * FR with the dissipation of FR for standard advection equation, the influence of IBM on the dissipation of the FR scheme can be studied. In the following investigation, for eigensolution analysis, we will concentrate on the primary mode in the semi-discrete analysis, and the primary and IBM modes in the fully-discrete analysis, since the latter is the main cause of instability. Moreover, to compare the dissipation behavior more reasonably, the IBM-induced dissipation γ IBM of the physical mode and the short-term dissipation should be subtracted, in order to only look at the dissipation effect of IBM on FR itself (γ * FR ). The influence of the penalization parameter is studied first. We consider a finer discretization with N = 40 elements in the computational domain (h = 0.05). We set ∆ = h with r = 1/40 and P = 3. Three penalization parameters 1 × 10 −3 , 1 × 10 −4 and 1 × 10 −5 are selected, which are relatively small because we need η to be small enough to guarantee good accuracy [42]. The results are compared against the standard analysis obtained for the standard advection equation, which are shown in Figure 3. As explained, the initial dissipation (γ IBM ) is subtracted for comparison purpose, in order to only look at the dissipation effect of the numerical scheme, as shown in Figure  3b. It can be seen that in well resolved regions (small to medium wavenumber), the dispersion-dissipation behaviors similar with the standard advection scheme are recovered, while a slight difference is seen in the high wavenumber region. This observation applies for both dispersion-dissipation and non-modal analyses, as shown in Figure 3a, 3b, and 3d. This is expected because it indicates that IBM does not affect too much the underlying FR scheme in the resolved wavenumber range. Another information is about total dissipation (both the IBM-induced dissipation and the dissipation of numerical scheme) for the physical mode, as shown in Figure 3c. It is clear that as η → 0, more damping is added to the primary mode, indicating that the solution will be further damped and the boundary condition will be imposed more strongly. This analysis can justify the use of volume penalization, since this method affects the spectral behavior of the underlying numerical scheme marginally, but only adds constant dissipation inside the wall to mimic the no-slip wall boundary condition.
We then look at the dispersion-dissipation and non-modal behaviors changing with polynomial order P, in order to show the consistency with the standard FR scheme. We fix the penalization parameter to η = 1 × 10 −4 , and set N = 40 with r = 1/40. The results are shown in Figure 4, where the IBM-induced dissipation is subtracted to show γ * FR . Good agreement with the standard advection equation in the resolved wavenumber range is shown. This range is also increasing as we go to higher P. These results indicate that volume penalization guarantees a good spectral behavior for high-order methods in the resolved wavenumber range, thus can be a good candidate for high-order schemes involving IBM treatment.

Fully-discrete analysis
Fully-discrete analysis considers both time and space discretization and is able to evaluate the stability of the space-time system. To march the solution in time, we choose a standard explicit three-stage, third-order Runge-Kutta scheme for time integration, which leads to the following fully-discrete operator matrix [23,76] A = I + ∆tM + 1 2 (∆tM) 2 + 1 6 (∆tM) 3 .
Note here we take the standard advection equation M as an example, while it can be easily generalized to penalized advection equation by changing the discretization operator matrix to the one in Equation 15. In addition, since the short-term dissipation is only defined when t → 0, non-modal analysis will not be considered. For a given k the for P = 1, 2, 3. c) Short-term dissipation subtracted by the value at k = 0 for P = 1, 2, 3. d) Dispersion for P = 4, 5, 6. e) Dissipation subtracted by the IBM-induced dissipation γ IBM (dissipation at k = 0) for P = 4, 5, 6. f) Short-term dissipation subtracted by the value at k = 0 for P = 4, 5, 6. eigenvalues related to the modified wavenumber are given by λ m = e −iω * m ∆t , therefore we have the numerical dispersion and dissipation defined as The fully-discrete dispersion and dissipation curves are again normalized by the smallest length scale h/(P+1) and the initial wavenumber is again re-scaled based on r. Since the dissipation curve obtained from fully-discrete analysis indicates the stability of the space-time scheme, it can help to investigate the stability criterion for the selection of penalization parameter. The scheme becomes unstable when any one of the eigenmodes shows positive dissipation. Again, we consider the computational domain discretized into N = 40 elements with one solid element r = 1/40. The polynomial order P = 3 is used which is a typical value used for high-order schemes. To set the penalization parameter we firstly consider the guideline η = ∆t [49,51]. This case is compared with results from [23], where fully-discrete analysis of the advection equation based FR discretization is considered. When time discretization is considered, the effect of time step ∆t should be taken into account, which is determined by the Courant-Friedrichs-Lewy (CFL) condition. The CFL number is defined as [77,78] Following [23], the maximum time-step for stability is firstly determined and given by CFL(max). Then we select three typical CFL numbers through multiplying this CFL number by 0.5, 0.7 and 0.9. The physical modes are compared in Figure 5, where the dispersion-dissipation behavior between IBM and standard advection equation shows good agreement. Very slight difference is seen when k > 2, which again belongs to the under-resolved wavenumber region. These results indicate that IBM also has very slight impact on the spectral behavior when both space and time discretizations are considered, and the penalization parameter is chosen within the stability limit.
Furthermore it has been shown that even though the time step is chosen below the critical CFL number, instability can still occur when η is too small since the problem becomes very stiff in the limit η → 0 [47]. Therefore, the critical value for penalization parameter, with respect to the time step, needs to be carefully studied to guide the selection of penalization parameter. From previous studies, we introduced some variations to η around η = ∆t and look at the instability behavior based on the fully-discrete dissipation curve. As shown in Figure 2, there are several solid modes among secondary modes induced by volume penalization, with constant dispersion and dissipation across all wavenumbers. When the penalization parameter varies, these solid modes will move and once η becomes smaller than the critical value η critical , these modes will show positive dissipation, indicating the numerical instability occurs. In order to investigate the range of η critical and obtain more accurate estimation of this parameter, the variation of both penalization parameter and the CFL number is studied. The dispersion behavior of the physical mode, and the dissipation behavior of the physical and the most unstable solid mode at two CFL numbers are shown in Figure 6 and Figure 7. This solid mode is identified as the most unstable eigenmode with maximal dissipation among all eigenmodes that have a constant value across the wavenumber range. It should be noted that this mode can be merged into other secondary modes when a relatively large η is selected with very large negative dissipation, and it will appear when η is approaching ∆t. As shown in the figures, when η decreases from 0.5∆t to 0.4∆t, the dissipation of this mode becomes positive, indicating that instability occurs. Therefore we can conclude that the instability of the numerical scheme is mainly driven by the instability of the solid mode introduced by IBM. This is consistent with the basic understanding of IBM, where the instability is induced by the stiffness of the IBM source term. Moreover, it also agrees with previous observations [13] that the secondary modes are responsible of the instability, rather than the primary mode. Finally, we compute the dissipation of the solid mode under different CFL numbers and compare the dissipation in Figure 8. It is found that for different CFL numbers, the critical η for stability based on the present time integration is always about 0.4∆t < η critical < 0.5∆t. As the CFL number increases, this value will slightly increase within this range. Compared with previous guidance of η = ∆t, this gives a more relaxed guidance to select the penalization parameter η.

Advection equation with volume penalization
In this subsection, the numerical simulation of the one-dimensional advection equation is presented, in order to validate the conclusions drawn in previous sections. We select the same case as before, where N = 40, ∆ = h (r = 1/40) and P = 3. The final simulation time is set to 1.1, to guarantee that the solution in the computational domain x ∈ [∆, 1] is sufficiently penalized. Since we consider the no-slip wall boundary condition u s = 0, the expected solution in this domain should be 0. The time integration is based on the third-order Runge-Kutta scheme used for fully-discrete analysis. To reduce the temporal error, a sufficiently small time step is set as ∆t = 1 × 10 −5 . Three initial conditions ranging from low to medium wavenumbers, including kh/(P + 1) = 0.3927, kh/(P + 1) = 0.7854 and kh/(P + 1) = 1.9635, are considered as test cases. The simulation results with different penalization parameters η are shown in Figure 9. It can be seen that as η decreases, more damping inside the solid is imposed, therefore the solution shows a smaller amplitude. This agrees with the observation in Figure 3c where more damping is seen on the primary mode when η is decreasing. In addition, the influence on frequency is not evident among these cases, indicating the fact that volume penalization only provides additional damping but does not influence significantly the dispersion behavior. It should also be noted that the resulting amplitude among these cases is due to the combined dissipation effects of IBM and the FR scheme, where the pure damping effect from IBM can be extracted by considering a constant initial condition. As η → 0, the solution will converge to zero, indicating that the incoming wave is sufficiently damped by the wall. However, due to the stiffness of the source term, as indicated from the fully-discrete analysis, it is not feasible to make η = 0 but to use a small η instead that balances the accuracy and stability. Moreover, as we go from low to high wavenumber, the amplitude of the solution for a given η decreases. This is due to the dissipation effect of the FR scheme, since IBM only has a constant damping across all wavenumbers.
The second set of test cases helps to validate the conclusions obtained from the fully-discrete analysis. We take the same test case as in Figure 9a, where the initial wavenumber is set to kh/(P + 1) = 0.3927. Three CFL numbers are chosen for investigation, including 0.1, 0.5, and 0.7, which correspond to time step ∆t = 6.5 × 10 −4 , 0.00325 and 0.00455, respectively. The penalization parameter is gradually reduced from 2∆t to 0.4∆t, with some typical solutions shown in Figure 10. From the figure it can be seen that the numerical method remains stable until η = 0.5∆ for all CFL numbers, which agree with the fully-discrete analysis. Instabilities are observed when η approaches 0.4∆t, whose critical value increases with increasing CFL number. Another interesting observation is that for large CFL numbers, see Figure 10b and Figure 10c, reducing η does not lead to a reduction of the error, which means for the present problem that there is no need to set η too small when the time step is large.

Improvement of the volume penalization through adding second-order terms inside the solid
From the eigensolution analysis and the numerical experiments, it can be concluded that the IBM works as a porous medium to provide damping inside the solid, where the damping effect can be evaluated through the eigensolution analysis. Therefore, alternative methods that provide additional damping in the solid region can be considered to develop better IBM schemes, like adding artificial viscosity or other wave absorption approaches. Note that all of these strategies must be imposed inside the solid to avoid modifying the real physics. In the present study we use a second-order term to provide additional damping to better satisfy Dirichlet boundary conditions. This strategy has also been mentioned in [45], where the constitutive equations are removed from the solid region and the nonphysical diffusion term is introduced to ensure the smoothness and continuity of the solution. This is different from the present work since the constitutive equation still exists for the solid body. This results in the following equation where η v is the diffusivity coefficient. We focus on two types of discretization with N = 40 and N = 80 elements, where the solid domain is set to be the same where ∆ = 0.05, indicating ∆ = h and ∆ = 2h for each case, respectively. The polynomial order and penalization parameter are set to P = 3 and η = 1 × 10 −3 . The viscosity parameter η v is varied to study the effect of additional dissipation. This viscous term is discretized by the Local DG (LDG) scheme [79], where two parameters in the scheme are set to β = 0.5 and τ = 0.1. Eigensolution analysis with different η v is firstly performed for the case with N = 40, as shown in Figure 11. Firstly, it indicates that as the viscous term is added, the dispersion behavior is only affected in high wavenumber region. Secondly, from Figure 11c, adding the second-order term will lead to smaller dissipation than the case with only volume penalization. This indicates the IBM treatment can be further enhanced by adding artificial viscosity inside the solid, making the solution further converge to the expected value. Thirdly, it should also be noted that although increasing η v will lead to larger total dissipation (as shown by short-term diffusion in Figure 11d), the dissipation of the physical mode does not increasing accordingly. This indicates that an optimal η v may exist which gives the largest diffusion to the physical mode (thus the solution).
After the eigensolution analysis, we perform numerical simulation for both test cases under different η v , where the final time is set to 1.1 and the initial condition is a sinusoidal wave with wavenumber kh/(P + 1) = 0.3142. Results for both cases are shown in Figure 12 and Figure 13, respectively. In Figure 12, it is observed that for N = 40, when η v is increasing, the solution is firstly damped until about η v = 5 × 10 −2 , where the performance is nearly optimal and the solution is quite close to zero. However, if η v continues to increase, the solution starts to grow with another phase. This is reflected by the eigensolution as well, where the wavenumber should be shifted by π/N before scaling. This means that for the present case, there is no need to increase η v when η v reach 5 × 10 −2 . For the second group of results shown in Figure 13, it can be seen that when very small η v is set, the solution can be improved, but further increasing η v will make the solution less accurate, e.g., when η v > 3 × 10 −2 . These results indicate that there exists an optimal combination of η and η v that gives the best numerical accuracy. These optima will be explored in what follows.    To further investigate the relationship between accuracy and viscosity parameter, the errors are compared. This error is defined as the error between the solution in computational domain from the right solid boundary to the right boundary of the whole computational domain (x ∈ [∆, 1]) and the penalized value u s = 0. Defining the number of  solution points inside the domain of interest as L, we have the mean squared error The simulation error at different η and η v is compared in Figure 14. It is clear that adding the second-order term will lead to lower error when the diffusivity parameter η v is chosen properly. For a given penalization parameter and problem setting, an optimal η v exists to provide the lowest error. However, a big gain is only seen for η = 1 × 10 −3 and η = 1 × 10 −4 , where the optimal η v will lead to a large reduction of the error. When the penalization parameter is sufficiently small, e.g., η = 1 × 10 −5 , the optimal performance requires a relatively larger η v and the gain becomes marginal. In addition, the critical time step ∆t for different η v is compared in Figure 15. This time step is the largest ∆t allowed for a stable time integration based on the present third-order Runge-Kutta time-marching scheme. It can be concluded that by adding the second-order term, less restrictions on time step are required for similar accuracy, which is a main advantage of this strategy. For example, for the first test case N = 40, the optimal accuracy for η = 1 × 10 −4 and η = 1 × 10 −5 is almost the same, while the time step can be about ten times larger (3.1 × 10 −5 against 3.1 × 10 −6 ). In general, the strategy proposed here is useful both in terms of improving the accuracy and relaxing the time step restriction which arises from the stiff source term associated to the volume penalization. Finally, the relationship between optimal η v (defined as η v,opt ) and η is explored. From Figure 14, the optimal η v (indicated by black symbols) increases as η is decreasing, indicating that there exists a dependency between these parameters. To further study this behavior, we consider simulations for different penalization parameters, polynomial orders and solid ratios, as shown in Figure 16a. We fix one solid element in the middle and choose three ratios r including 1/20, 1/40 and 1/80 (representing different solid widths), three P ranging from 2 to 4 and three penalties η including 1 × 10 −3 , 1 × 10 −4 and 1 × 10 −5 . For each parameter pair, an initial condition with a low wavenumber is chosen to represent the resolved wavenumber region of the scheme. Numerical simulations across a group of η v are performed to numerically find η v,opt . From Figure 16a, a nearly linear relation is seen between η v,opt and η. For a given r and η, η v,opt remains approximately constant for polynomial order 2 and 3, while it reduces by half for polynomial order 4. Therefore, scaling relationships of ηη v,opt /(r 2 ) = 0.09 and ηη v,opt /(r 2 ) = 0.05 can be fitted for P = 2/3 and P = 4, respectively. These fitted scaling curves are shown in Figure 16 b. The fitted scaling relation can provide guidelines for the selection of η v,opt based on η and P.

Flow past a circular cylinder
To further investigate the observations from the proposed analyses, volume penalization for Navier-Stokes equations based on high-order FR discretization is considered [53,54]. As a benchmark test case, steady flow past a static circular cylinder at Reynolds number 40 is chosen. For compressible Navier-Stokes equations, the following source term is added to the right-hand-side of the momentum equations [46]: where ρ, u and v denote the density and velocities in x and y directions, respectively. u s and v s refer to the penalized velocities inside the solid, which are fixed to zero in the present study to mimic the no-slip wall boundary condition of a cylinder. An additional source term added to the energy equation is also needed following [46] More details about the implementation, e.g., surface geometry representation and computation of integral aerodynamic coefficients, are given in [53,54]. To apply the scaling law proposed in the previous section, the solid ratio r should be defined for two-dimensional flows. Here we generalize the solid ratio for two-dimensional flows as where S solid and S domain are areas of the solid region and the whole computational domain, respectively. five-stage fourth-order explicit Runge-Kutta method (LSERK) [80,63] is used to march the solution in time. The polynomial order is set to P = 2 and the time step ∆t is varied according different combinations of parameters. At first, the volume penalization is applied, where we set η = ∆t from the traditional guideline [47,51], which gives η = ∆t = 2 × 10 −4 . Based on this case, we find another case where similar time step / computational cost is required but may lead to a smaller error, and set a smaller penalization parameter η = 1 × 10 −3 . From the relationship obtained in the previous section for P = 2, i.e., ηη v,opt /(r 2 ) = 0.09, we can obtain η v,opt ≈ 1.5 × 10 −2 , where r ≈ 0.0128 is obtained from Equation 29. Therefore, three additional simulations with η = 1 × 10 −3 and η v = 1.5 × 10 −2 , 1.0 × 10 −2 and 2.0 × 10 −2 are included for comparison, where the same time step ∆t = 2 × 10 −4 is used for the first two cases, while ∆t = 1.5 × 10 −4 is needed for the last case to ensure stability.
In Figure 17, the temporal evolution of the lift and drag coefficients is compared. Due to the symmetric geometry of the cylinder and the selected Reynolds number, the expected lift coefficient should approach zero. As shown in Figure 17a, the optimal prediction is given by Case 2. This is achieved by the combination of η and η v , where η v is given by the approximated relationship. This highlights that it is possible to utilize the proposed relationship to help the selection of penalization and diffusivity parameters. In addition, different combinations of these parameters also have an impact on the predicted drag coefficient, as shown in Figure 17b. Note that the predicted drag coefficients (around 1.48) agree with the results reproduced from existing literature, e.g., [45,81]. To further evaluate the proposed approach and the approximated relationship, more simulations are conducted, and the predicted lift coefficient at different combinations of penalization and diffusivity parameters are compared in Figure 18. Variations of η v at three penalization parameters, including 1 × 10 −2 , 1 × 10 −3 and 2 × 10 −4 are considered. It can be seen that for small η, there exists an optimal η v that gives the best lift coefficient. When η = 1 × 10 −3 , the proposed relationship accurately predicted the η v,opt , as validated in Figure 17. At large η where the solution is not sufficiently penalized, adding the second-order term will lead to a slightly larger error, while the predicted η v,opt still gives reasonable result since a relatively small η v = 1.5 × 10 −3 is obtained. When η = 2 × 10 −4 , although the optimal value is not well predicted, η v,opt = 7 × 10 −2 still gives a local minimum error. These results indicate that the proposed relationship has the potential to be used as a guideline for Navier-Stokes simulations, however, further study is still needed. In general, the results in this section show the efficacy of the proposed approach to include second-order term to improve the representation of the boundary condition.

Conclusions
This paper applies eigensolution and non-modal analyses for the immersed boundary method in the advection equation, where a high-order flux reconstruction scheme based on volume penalization is investigated. Three main conclusions are drawn: 1) The semi-discrete eigensolution analysis and non-modal analysis show that volume penalization does not incur in significant changes in the dispersion behavior of the underlying high-order scheme in the resolved wavenumber range. It imposes the boundary condition through introducing numerical dissipation inside the solid, in order to damp and drive the solution to the prescribed boundary conditions.
2) Fully-discrete analyses show that the cause of numerical instability is associated to the secondary modes, in particular the modes coming from the IBM treatment. This analysis further reaches the stability condition of penalization parameter for explicit time marching based on third-order Runge-Kutta scheme: 0.4∆t < η critical < 0.5∆t. This analysis provides guidelines to the selection of the penalization parameter, where a smaller η is preferred but should be limited to avoid increased system stiffness and associated numerical stabilities.
3) The proposed analysis shows that the IBM works as a porous medium to absorb and damp the solution inside the solid body. As a result, we propose to add a second order term inside the body. Results show that this approach, in combination with the volume penalization, achieves improved accuracy and reduces the time step restrictions due to stiff penalization terms. For a given penalization parameter, there exists an optimal viscosity coefficient that has the lowest numerical error. The relationship between volume penalization and optimal viscosity coefficient is obtained by curve-fitting. 4) We test our results on a Navier-Stokes flow and find that the proposed second-order damping inside the solid helps to achieve more accurate solutions. Furthermore, the relationship derived using linear analysis, is used as a guideline for choosing the penalty parameters. The cylinder results show improved accuracy, proving the potential of the presented methodology for Navier-Stokes simulations.

Acknowledgments
This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement (MSCA ITN-EID-GA ASIMIA No 813605).

Appendix A. Discretization of the Advection Equation
Appendix A. 1

. Preliminaries
As the first step of space discretization for a one-dimensional problem, the computational domain Ω is divided into N nonoverlapping cells, defined as Ω n = {x|x n < x < x n+1 } with element size h n = x n+1 − x n . To implement the space discretization, each cell is transformed from the physical space within the global domain to a common reference space Ω re f = {r| − 1 < r < 1}. The mapping between physical space and reference space can be defined as the function Γ r = Γ n (x) = 2 x − x n x n+1 − x n − 1, (A.1) and its inversion This transformation leads to the transformed conservation law: where the Jacobian J n is defined as J n = (x n+1 − x n )/2. For each cell, the solution is represented by a polynomial of degree P defined at P + 1 solution points (quadrature points). The polynomial can be represented by a vector of modal coefficientsû j of a Legendre basis, or a vector of nodal coefficientsũ j defined at solution points. Here we retain the former set of basis.
Appendix A.2. Spatial discretization The standard FR workflow for a general conservation law includes seven stages [82], whereas for the advection equation it can be reduced to five stages. In the first and the second stage, the polynomial representation in terms of a nodal basis (ũ(r i ), i = 0, 1, ..., P) with order P for both solution and flux is defined at P + 1 solution points: u n (r) = where d is the gradient operator to obtain the discontinuous gradient at position r, which can be generalized to gradient matrix D at all solution points D[i, j] = ∂l j (r i ) ∂r . In the third stage of FR, the transformed numerical fluxf I at all element interfaces (either ends of the cell Ω r , r = ±1) is computed, which is a common value for the two adjacent cells of each interface. To obtain this flux, the interpolated solution at either end of the interface is needed. For example, for the nth cell, we have the left and right interpolated solutionsũ L n = l T (−1)ũ n andũ R n = l T (1)ũ n , as well as the interpolated fluxesf D,L n = cũ L n andf D,R n = cũ R n . The exact method for numerical flux calculationf I =f I (ũ L ,ũ R ) depends on the nature of the equation being solved, like a Roe type approximate Riemann solver [83] or other upwind flux solvers for the Euler equations. Here we use the Lax-Friedrich type numerical flux defined as where λ is an upwinding parameter that controls the numerical flux. λ = 1 defines the upwind flux and λ = 0 defines the central flux. The upwind flux is considered here. The numerical interface fluxes for both boundaries of the nth cell are defined asf I,L n andf I,R n , respectively. The forth stage of FR process is to correct the transformed discontinuous flux to be continuous across the cell interface, while approximating the flux function as a polynomial of degree P + 1. The transformed correction fluxf C based on a polynomial of degree P + 1 is defined to make the sum betweenf C and the interpolated discontinuous flux equal to the transformed numerical interface flux at r ± 1: To define this correction flux, considering first the correction functions for the left and the right interface, defined as g L = g L (r) and g R = g R (r). This correction satisfies the following condition g L (−1) = 1 , g L (1) = 0 , g R (−1) = 0 , g R (1) = 1, (A.11) with a symmetry consideration g L (r) = g R (−r). (A.12) The specific form of correction function is determined by the above-mentioned condition and the stability criterion. Here the left and right Radau polynomials are used for the right and left boundaries, respectively. This type of polynomial belongs to a class of high-order energy stable flux reconstruction (ESFR) schemes that satisfy energy stability criterion. As mentioned in Vincent et al. [69], the resulting scheme is equivalent to a differential formulation of a nodal DG method. The correction function is then written in terms of g L and g R as f C n (r) = (f I,L n −f D,L n )g L (r) + (f I,R n −f D,R n )g R (r). (A.13) Therefore the total transformed fluxf with a polynomial of degree P + 1 is formed from the discontinuous and correction flux as followsf n =f D n +f C n . (A.14) In the final stage of FR process, the divergence of the transformed total flux ∂f ∂r at each solution point is calculated as follows ∂f n ∂r (r i ) = Finally, the divergence of total flux can be used to advance the solution u in time with some efficient time integration methods, like the Runge-Kutta method. For viscous fluxes, the Local Discontinuous Galerkin (LDG) [79] formulation can be adopted. In summary, it can be concluded that the performance of the present FR scheme depends on three factors, including the position of solution points, the method for solving the transformed interface flux (the Riemann solver) and the form of correction function.