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

Jiaqing Kou<EMAIL_ADDRESS>Aurelio Hurtado-de-Mendoza Saumitra Joshi
Soledad Le Clainche Esteban Ferrer ETSIAE-UPM-School of Aeronautics,
Universidad Politécnica de Madrid, Plaza Cardenal Cisneros 3, E-28040 Madrid,
Spain  NUMECA International S.A., Chaussee de la Hulpe 187, Brussels, B-1170,
Belgium  Center for Computational Simulation, Universidad Politécnica de
Madrid, Campus de Montegancedo, Boadilla del Monte, 28660 Madrid, Spain

###### Abstract

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.

###### keywords:

volume penalization, flux reconstruction, immersed boundary method, high-order
methods, von Neumann analysis, non-modal analysis

††journal: Journal of Computational Physics

###### Contents

  1. 1 Introduction
  2. 2 Methodology
    1. 2.1 The governing equation
    2. 2.2 Immersed boundary method based on volume penalization
  3. 3 Modal and non-modal analysis
    1. 3.1 Eigensolution analysis
    2. 3.2 Non-modal analysis
    3. 3.3 Semi-discrete analysis
    4. 3.4 Fully-discrete analysis
  4. 4 Numerical experiments
    1. 4.1 Advection equation with volume penalization
    2. 4.2 Improvement of the volume penalization through adding second-order terms inside the solid
    3. 4.3 Flow past a circular cylinder
  5. 5 Conclusions
  6. A Discretization of the Advection Equation
    1. A.1 Preliminaries
    2. A.2 Spatial discretization
    3. A.3 Local matrix formulation

##  1 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 $\Delta
t<\IBMparam$ 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 fully-discrete 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.

##  2 Methodology

###  2.1 The governing equation

We consider a one-dimensional advection equation defined in space $x$ and time
$t$:

$\frac{\partial\velx}{\partial t}+\frac{\partial f}{\partial x}=0,$ (1)

where $\velx(x,t)$ is the transported unknown solution and $f=\advcoef\velx$
is the flux, with advection speed $\advcoef$. 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)=\text{exp}(ikx)$:

$u(x,t)=\text{exp}[i(kx-\omega t)],$ (2)

where $k$ is the wavenumber, $\omega=\advcoef k$ is the (angular) frequency
and $i=\sqrt{-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 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.

###  2.2 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 $\mask(x,t)$ that distinguishes between the fluid region $\Omega_{f}$
and solid region $\Omega_{s}$ is defined first:

$\mask(x,t)=\left\\{\begin{matrix}1,\,\,\text{if}\,\,x\in\Omega_{s}\\\
0,\,\,\text{otherwise}\end{matrix}\right..$ (3)

For the one-dimensional advection equation with penalized Dirichlet boundary
conditions, this problem is formulated as follows:

$\frac{\partial\velx}{\partial t}+\advcoef\frac{\partial\velx}{\partial
x}+\frac{\mask}{\IBMparam}(\velx-\velx_{s})=0,$ (4)

where $\velx_{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 $\IBMparam$, which drives
the solution to $\velx_{s}$ as $\IBMparam\rightarrow 0$. For solution inside
the solid, this results in the following governing equation similar with
convection-diffusion-reaction equations [71]

$\frac{\partial\velx}{\partial t}+\advcoef\frac{\partial\velx}{\partial
x}+\frac{1}{\IBMparam}\velx=0.$ (5)

The analytical solution of this equation is

$\velx(x,t)=\text{exp}[i(kx-(\advcoef k-\frac{i}{\IBMparam})t)],$ (6)

indicating that the penalization term introduces additional damping to the
solution, controlled by the penalization parameter $\IBMparam$. 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 $\velx_{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 $\eta$ 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]

$\left\|\velx^{exact}-\velx_{\eta}^{N}\right\|\leq\left\|\velx^{exact}-\velx_{\eta}^{exact}\right\|+\left\|\velx_{\eta}^{exact}-\velx_{\eta}^{N}\right\|,$
(7)

where the first part is the penalization error and the second part is the
discretization error. $u^{exact}$, $u_{\eta}$ and $u_{\eta}^{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 $\left\|\cdot\right\|$ between solutions. The penalization error
depends on the penalization parameter [47]:

$\left\|\velx^{exact}-\velx_{\eta}^{exact}\right\|\propto\eta^{\alpha},$ (8)

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_{\eta\rightarrow
0}\left\|u^{exact}-u_{\eta}^{exact}\right\|\rightarrow 0$. It has been proved
that the volume penalization gives $\alpha=\frac{1}{2}$ for the Dirichlet
boundary conditions, indicating the penalization error has a decay rate of
$\mathcal{O}(\sqrt{\eta})$. For Neumann boundary condition,
$\mathcal{O}(\eta)$ 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 ($\beta>0$):

$\left\|u_{\eta}^{exact}-u_{\eta}^{N}\right\|\propto N^{-\beta}.$ (9)

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_{\eta}$ at the boundary of the penalized
equation. This regularity can be improved by designing proper $\velx_{s}$
inside the solid [73, 74]. However, in this study we will discuss the
classical volume penalization method with $\velx_{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 $\eta$ which
is small enough to ensure low penalization error (high penalization since
$\eta$ is in the denominator). In the meantime, $\eta$ 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 $\Delta t<\eta$ [49] or $\Delta t\approx\eta$ [75] to ensure stability. An
interpretation for this choice is that the penalty term acts as a strong
damping term with order $\eta$ 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.

##  3 Modal and non-modal analysis

###  3.1 Eigensolution analysis

In the present temporal eigensolution analysis, we assume periodic boundary
condition and consider solution $u(x,t)=\text{exp}[i(kx-\omega t)]$ with real
wavenumber $k$. Therefore, for each element with index $n$, we have
$\boldsymbol{u}_{n-1}=e^{-ikh}\boldsymbol{u}_{n}$ and
$\boldsymbol{u}_{n+1}=e^{+ikh}\boldsymbol{u}_{n}$. The space discretization
results in

$\frac{d\boldsymbol{u}_{n}}{dt}=\left(\boldsymbol{L}e^{-ikh}+\boldsymbol{C}+\boldsymbol{R}e^{+ikh}\right)\boldsymbol{u}_{n}=\boldsymbol{M}\boldsymbol{u}_{n}.$
(10)

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

$-i\omega^{*}\boldsymbol{u}_{n}=\boldsymbol{M}\boldsymbol{u}_{n},$ (11)

where $\omega^{*}=\advcoef k^{*}$ becomes a complex value due to the
dispersion and dissipation errors of the space discretization. This wavenumber
$\omega^{*}$ relates to the eigendecomposition of the coefficient matrix
$\boldsymbol{M}$ (with $P+1$ solutions)

$-i\omega^{*}_{m}=\lambda_{m}\\\ ,\
\lambda_{m}\boldsymbol{v}_{m}=\boldsymbol{M}\boldsymbol{v}_{m},$ (12)

where $\lambda_{m}$ and $\boldsymbol{v}_{m}$ are the $m$th eigenvalues and
eigenvectors of matrix $\boldsymbol{M}$, respectively. By defining the
modified wavenumber $k^{*}_{m}$ for each $\omega^{*}_{m}$ and $\lambda_{m}$,
we have:

$\text{Real}(k^{*}_{m})=\frac{\text{Real}(\omega^{*}_{m})}{\advcoef}=-\frac{\text{Imag}(\lambda_{m})}{\advcoef}\\\
,\
\text{Imag}(k^{*}_{m})=\frac{\text{Imag}(\omega^{*}_{m})}{\advcoef}=\frac{\text{Real}(\lambda_{m})}{\advcoef}.$
(13)

For the advection equation, the difference between $\text{Real}(k^{*})$ and
$k$ is due to the dispersion error, indicating the change in wavenumber of the
solution. The difference between $\text{Imag}(k^{*})$ and $0$ corresponds to
the dissipation (diffusion) error, where $\text{Imag}(k^{*})\leq 0$ holds for
a stable scheme. Note that for the analytical solution we have
$\text{Real}(k^{*})=k$ and $\text{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, $\boldsymbol{M}$ will
have $P+1$ eigenmodes, where the one that recovers $k^{*}=k$ as $k\rightarrow
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
$\boldsymbol{U}=[\boldsymbol{u}_{1},\boldsymbol{u}_{2},...,\boldsymbol{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, $\boldsymbol{u}_{0}$ and $\boldsymbol{u}_{N+1}$, derived as
$\boldsymbol{u}_{0}=\text{exp}[-2ikT]\boldsymbol{u}_{N}$ and
$\boldsymbol{u}_{N+1}=\text{exp}[2ikT]\boldsymbol{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:

$\frac{d}{dt}\begin{pmatrix}\boldsymbol{\velx}_{1}\\\
\boldsymbol{\velx}_{2}\\\ \boldsymbol{\velx}_{3}\\\ \vdots\\\
\boldsymbol{\velx}_{N}\end{pmatrix}=\begin{pmatrix}\boldsymbol{C}-\frac{\chi}{\eta}\boldsymbol{I}&\boldsymbol{R}&\boldsymbol{0}&\cdots&\text{exp}[-2ikT]\boldsymbol{L}\\\
\boldsymbol{L}&\boldsymbol{C}-\frac{\chi}{\eta}\boldsymbol{I}&\boldsymbol{R}&\cdots&\boldsymbol{0}&\\\
\boldsymbol{0}&\boldsymbol{L}&\boldsymbol{C}-\frac{\chi}{\eta}\boldsymbol{I}&\cdots&\boldsymbol{0}\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
\text{exp}[2ikT]\boldsymbol{R}&\boldsymbol{0}&\boldsymbol{0}&\cdots&\boldsymbol{C}-\frac{\chi}{\eta}\boldsymbol{I}\end{pmatrix}\begin{pmatrix}\boldsymbol{\velx}_{1}\\\
\boldsymbol{\velx}_{2}\\\ \boldsymbol{\velx}_{3}\\\ \vdots\\\
\boldsymbol{\velx}_{N}\end{pmatrix},$ (14)

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 ($\mask=1$). For example, assuming we have one solid element in the
middle, the global matrix becomes

$\frac{d}{dt}\begin{pmatrix}\boldsymbol{\velx}_{1}\\\
\boldsymbol{\velx}_{2}\\\ \vdots\\\ \boldsymbol{\velx}_{IBM}\\\ \vdots\\\
\boldsymbol{\velx}_{N}\end{pmatrix}=\begin{pmatrix}\boldsymbol{C}&\boldsymbol{R}&\boldsymbol{0}&\cdots&\cdots&\cdots&\text{exp}[-2ikT]\boldsymbol{L}\\\
\boldsymbol{L}&\boldsymbol{C}&\boldsymbol{R}&\cdots&\cdots&\cdots&\boldsymbol{0}&\\\
\vdots&\vdots&\vdots&\ddots&\ddots&\ddots&\vdots\\\
\vdots&\vdots&\vdots&\boldsymbol{L}&\boldsymbol{C}-\frac{1}{\eta}\boldsymbol{I}&\boldsymbol{R}&\vdots\\\
\vdots&\vdots&\vdots&\ddots&\ddots&\ddots&\vdots\\\
\text{exp}[2ikT]\boldsymbol{R}&\boldsymbol{0}&\boldsymbol{0}&\cdots&\cdots&\cdots&\boldsymbol{C}\end{pmatrix}\begin{pmatrix}\boldsymbol{\velx}_{1}\\\
\boldsymbol{\velx}_{2}\\\ \vdots\\\ \boldsymbol{\velx}_{IBM}\\\ \vdots\\\
\boldsymbol{\velx}_{N}\end{pmatrix}.$ (15)

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\rightarrow 0$ [11]). All remaining modes are referred
to as secondary modes.

###  3.2 Non-modal analysis

Non-modal analysis [21] is an approach to analyze the numerical diffusion,
which is more relevant in under-resolved 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
short-term diffusion parameter $\tilde{\omega}^{*}$ is defined as [21]

$\tilde{\omega}^{*}:=\left(\frac{d\log||\boldsymbol{\velx}_{n}||}{d\tau^{*}}\right)_{\tau^{*}=0},$
(16)

where $||\cdot||$ is the $L_{2}$ norm and
$\tau^{*}=\tau(P+1)=t\advcoef(P+1)/h$ is a non-dimensional time based on the
non-dimensional convection time per degree of freedom. For
$\tilde{\omega}^{*}\leq 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:

$\tilde{\omega}^{*}=\frac{1}{P+1}\text{Real}\left[{\frac{(\boldsymbol{\velx}_{n,0})^{\dagger}\,\boldsymbol{M}\,\boldsymbol{\velx}_{n,0}}{(\boldsymbol{\velx}_{n,0})^{\dagger}\boldsymbol{\velx}_{n,0}}}\right],$
(17)

where $\boldsymbol{\velx}_{n,0}$ is the initial condition, and $\dagger$
denotes the complex conjugate. $\boldsymbol{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 $\boldsymbol{M}$.

###  3.3 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 $\advcoef=1$ and a
computational domain defined in $x\in[-1,1]$ discretized by equispaced
elements. All the cases considered here are based on the fully upwind flux
with $\lambda=1$. The solid region lies in the middle, starting from $x=0$,
whose width is defined as $\Delta$:

$\mask(x,t)=\left\\{\begin{matrix}1,\,\,\text{if}\,\,0<x<\Delta\\\
0,\,\,\text{otherwise}\end{matrix}\right..$ (18)

In the present study, we will consider $\Delta$ as integer multiples of
element size $h$ ($\Delta=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 $\advcoef=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$

$\hat{k}=k/(1-r).$ (19)

This definition is consistent with the standard advection equation where we
have $r=0$ thus $\hat{k}=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 $\pi/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 wavelength $\hat{k}$ 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.

Fig. 1: Schematic illustration of the advection problem with IBM.

As the first example, we consider a coarse discretization, where the
computational domain is discretized into $N=10$ elements, with the solid width
$\Delta=h=0.2$ and $r=1/10$. The polynomial order is set to $P=3$, and the
penalization parameter is set to $\IBMparam=1\times 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 $\velx_{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 $\IBMparam$. 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 $\IBMparam$ which determines the decay of the solution inside the solid.
Following this analysis, we can define the IBM-induced dissipation
$\gamma_{IBM}$ for the primary mode as

$\gamma_{IBM}=Imag[k^{*}(0)]h/(P+1),$ (20)

where $k^{*}(0)$ refers to the modified wavenumber $k^{*}$ at $k=0$.
Therefore, the total dissipation $\gamma$ 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)

$\gamma=\gamma_{IBM}+\gamma_{FR}^{*},$ (21)

which also applies to short-term dissipation based on the non-modal analysis.
By comparing $\gamma_{FR}^{*}$ with the dissipation of FR for standard
advection equation, the influence of IBM on the dissipation of the FR scheme
can be studied.

(a)

(b)

(c)

(d)

Fig. 2: Eigensolution and non-modal analysis of advection equation with IBM
($N=10$, $r=1/10$, $P=3$ and $\IBMparam=1\times 10^{-4}$). a) Dispersion. b)
Dissipation. c) Zoom-in dissipation. d) Short-term dissipation. Red and yellow
curves represent the physical mode and the mode induced by the IBM source
term. Deshed line is the refernce curve $Real(k^{*})=k$, $Imag(k^{*})=0$ and
$\omega^{*}=0$.

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 $\gamma_{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 ($\gamma_{FR}^{*}$).

(a)

(b)

(c)

(d)

Fig. 3: Eigensolution and non-modal analysis of advection equation with IBM
($N=40$, $r=1/40$, $P=3$) and different penalization parameter $\eta$. a)
Dispersion. b) Dissipation subtracted by the IBM-induced dissipation
$\gamma_{IBM}$ (dissipation at $k=0$). c) Zoom-in dissipation. d) Short-term
dissipation subtracted by the dissipation at $k=0$. The reason for subtracting
the initial value is because this value depends on $\eta$ and is only an
effect of IBM source term.

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 $\Delta=h$ with $r=1/40$ and $P=3$. Three penalization
parameters $1\times 10^{-3}$, $1\times 10^{-4}$ and $1\times 10^{-5}$ are
selected, which are relatively small because we need $\IBMparam$ 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 ($\gamma_{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 $\eta\rightarrow 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 $\IBMparam=1\times 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 $\gamma_{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.

(a)

(b)

(c)

(d)

(e)

(f)

Fig. 4: Eigensolution and non-modal analysis of advection equation with IBM
($N=40$, $r=1/40$, $\IBMparam=1\times 10^{-4}$) and comparison with those of
standard advection equation. a) Dispersion for $P=1,2,3$. b) Dissipation
subtracted by the IBM-induced dissipation $\gamma_{IBM}$ (dissipation at
$k=0$) 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 $\gamma_{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$.

###  3.4 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]

$\boldsymbol{A}=\boldsymbol{I}+\Delta t\boldsymbol{M}+\frac{1}{2}(\Delta
t\boldsymbol{M})^{2}+\frac{1}{6}(\Delta t\boldsymbol{M})^{3}.$ (22)

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\rightarrow 0$, non-modal
analysis will not be considered. For a given $k$ the eigenvalues related to
the modified wavenumber are given by $\lambda_{m}=e^{-i\omega^{*}_{m}\Delta
t}$, therefore we have the numerical dispersion and dissipation defined as

$\text{Real}(k^{*}_{m})=\frac{\text{Real}(\omega^{*}_{m})}{\advcoef}=\frac{\text{Real}[i\ln(\lambda_{m})]}{\advcoef\Delta
t}\\\ ,\
\text{Imag}(k^{*}_{m})=\frac{\text{Imag}(\omega^{*}_{m})}{\advcoef}=\frac{\text{Imag}[i\ln(\lambda_{m})]}{\advcoef\Delta
t}.$ (23)

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 $\IBMparam=\Delta 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 $\Delta t$ should
be taken into account, which is determined by the Courant–Friedrichs–Lewy
(CFL) condition. The CFL number is defined as [77, 78]

$CFL=\advcoef\Delta t(2P+1)/h.$ (24)

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.

(a)

(b)

Fig. 5: Fully-discrete analysis of advection equation with IBM ($N=40$,
$r=1/40$, $P=3$, $\IBMparam=1\times 10^{-4}$) and comparison with those of
standard advection equation, under different CFL numbers. Third-order Runge-
Kutta method is considered for time integration. a) Dispersion. b) Dissipation
subtracted by the IBM-induced dissipation $\gamma_{IBM}$ (dissipation at
$k=0$).

Furthermore it has been shown that even though the time step is chosen below
the critical CFL number, instability can still occur when $\IBMparam$ is too
small since the problem becomes very stiff in the limit $\IBMparam\rightarrow
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 $\IBMparam$ around $\IBMparam=\Delta 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 $\IBMparam$ becomes smaller than the critical value
$\IBMparam_{critical}$, these modes will show positive dissipation, indicating
the numerical instability occurs. In order to investigate the range of
$\IBMparam_{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 $\IBMparam$ is
selected with very large negative dissipation, and it will appear when
$\IBMparam$ is approaching $\Delta t$. As shown in the figures, when
$\IBMparam$ decreases from $0.5\Delta t$ to $0.4\Delta 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.

(a)

(b)

Fig. 6: Fully-discrete analysis of advection equation with IBM ($N=40$,
$r=1/40$, $P=3$, $\IBMparam=1\times 10^{-4}$, CFL $=0.1$ CFL(max)) and
comparison with those of standard advection equation. Third-order Runge-Kutta
method is considered for time integration. a) Dispersion. b) Dissipation of
both primary mode and the most unstable solid mode.

(a)

(b)

Fig. 7: Fully-discrete analysis of advection equation with IBM ($r=1/40$,
$N=40$, $P=3$, $\IBMparam=1\times 10^{-4}$, CFL $=0.7$ CFL(max)) and
comparison with those of standard advection equation. Third-order Runge-Kutta
method is considered for time integration. a) Dispersion. b) Dissipation of
both primary mode and the most unstable solid mode. Fig. 8: Comparison of the
dissipation of the most unstable solid mode under different $\eta$ and CFL
number.

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 $\IBMparam$ for stability based on the
present time integration is always about $0.4\Delta
t<\IBMparam_{critical}<0.5\Delta t$. As the CFL number increases, this value
will slightly increase within this range. Compared with previous guidance of
$\IBMparam=\Delta t$, this gives a more relaxed guidance to select the
penalization parameter $\IBMparam$.

##  4 Numerical experiments

###  4.1 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$, $\Delta=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\in[\Delta,1]$ is sufficiently
penalized. Since we consider the no-slip wall boundary condition
$\velx_{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 $\Delta t=1\times 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 $\IBMparam$ are shown in Figure 9. It
can be seen that as $\IBMparam$ 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 $\IBMparam$ 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 $\IBMparam\rightarrow 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 $\IBMparam=0$ but to
use a small $\IBMparam$ instead that balances the accuracy and stability.
Moreover, as we go from low to high wavenumber, the amplitude of the solution
for a given $\IBMparam$ decreases. This is due to the dissipation effect of
the FR scheme, since IBM only has a constant damping across all wavenumbers.

(a)

(b)

(c)

Fig. 9: Simulation of advection equation with IBM under different penalization
parameters ($N=40$, $r=1/40$, $P=3$). a) $kh/(P+1)=0.3927$. b)
$kh/(P+1)=0.7854$. c) $kh/(P+1)=1.9635$.

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 $\Delta t=6.5\times 10^{-4}$, $0.00325$ and $0.00455$,
respectively. The penalization parameter is gradually reduced from $2\Delta t$
to $0.4\Delta t$, with some typical solutions shown in Figure 10. From the
figure it can be seen that the numerical method remains stable until
$\IBMparam=0.5\Delta$ for all CFL numbers, which agree with the fully-discrete
analysis. Instabilities are observed when $\IBMparam$ approaches $0.4\Delta
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 $\IBMparam$ does not lead to a reduction of the error,
which means for the present problem that there is no need to set $\IBMparam$
too small when the time step is large.

(a)

(b)

(c)

Fig. 10: Simulation of advection equation with IBM under different CFL numbers
and penalization parameters ($N=40$, $r=1/40$, $P=3$, $kh/(P+1)=0.3927$). a)
CFL = 0.1CFL(max). b) CFL = 0.5CFL(max). c) CFL = 0.7CFL(max).

###  4.2 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

$\frac{\partial\velx}{\partial t}+\advcoef\frac{\partial\velx}{\partial
x}=-\frac{\mask}{\IBMparam}\velx+\mask\IBMparam_{v}\frac{\partial^{2}\velx}{\partial
x^{2}},$ (25)

where $\IBMparam_{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 $\Delta=0.05$, indicating $\Delta=h$ and $\Delta=2h$ for
each case, respectively. The polynomial order and penalization parameter are
set to $P=3$ and $\IBMparam=1\times 10^{-3}$. The viscosity parameter
$\IBMparam_{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 $\beta=0.5$ and $\tau=0.1$. Eigensolution
analysis with different $\IBMparam_{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 $\IBMparam_{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 $\IBMparam_{v}$ may exist which gives the
largest diffusion to the physical mode (thus the solution).

(a)

(b)

(c)

(d)

Fig. 11: Eigensolution and non-modal analysis of advection equation with IBM
($N=40$, $r=1/40$, $P=3$, $\IBMparam=1\times 10^{-3}$) and different
diffusivity $\IBMparam_{v}$. a) Dispersion. b) Dissipation subtracted by the
IBM-induced dissipation $\gamma_{IBM}$ (dissipation at $k=0$). c) Zoom-in
dissipation. d) Short-term dissipation subtracted by the dissipation at $k=0$.

After the eigensolution analysis, we perform numerical simulation for both
test cases under different $\IBMparam_{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
$\IBMparam_{v}$ is increasing, the solution is firstly damped until about
$\IBMparam_{v}=5\times 10^{-2}$, where the performance is nearly optimal and
the solution is quite close to zero. However, if $\IBMparam_{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 $\pi/N$
before scaling. This means that for the present case, there is no need to
increase $\IBMparam_{v}$ when $\IBMparam_{v}$ reach $5\times 10^{-2}$. For the
second group of results shown in Figure 13, it can be seen that when very
small $\IBMparam_{v}$ is set, the solution can be improved, but further
increasing $\IBMparam_{v}$ will make the solution less accurate, e.g., when
$\IBMparam_{v}>3\times 10^{-2}$. These results indicate that there exists an
optimal combination of $\IBMparam$ and $\IBMparam_{v}$ that gives the best
numerical accuracy. These optima will be explored in what follows.

(a)

(b)

Fig. 12: Simulation of advection equation with IBM under different diffusivity
coefficients ($N=40$, $r=1/40$, $P=3$, $\IBMparam=1\times 10^{-3}$, initial
wavenumber $kh/(P+1)=0.3142$). a) Global view. b) Zoom-in view.

(a)

(b)

Fig. 13: Simulation of advection equation with IBM under different diffusivity
coefficients ($N=80$, $r=1/40$, $P=3$, $\IBMparam=1\times 10^{-3}$, initial
wavenumber $kh/(P+1)=0.3142$). a) Global view. b) Zoom-in view.

(a)

(b)

Fig. 14: Simulation error versus diffusivity coefficient ($r=1/40$, $P=3$,
initial wavenumber $kh/(P+1)=0.3142$). a) $N=40$. b) $N=80$. Black symbols
refer to the lowest error with respect to a given penalization parameter
$\IBMparam$.

(a)

(b)

Fig. 15: Critical time step versus diffusivity coefficient ($r=1/40$, $P=3$,
initial wavenumber $kh/(P+1)=0.3142$). a) $N=40$. b) $N=80$. Black symbols
refer to the lowest error with respect to a given penalization parameter
$\IBMparam$.

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\in[\Delta,1]$) and the
penalized value $\velx_{s}=0$. Defining the number of solution points inside
the domain of interest as $L$, we have the mean squared error

$error=\sqrt{\frac{\sum_{i=1}^{L}[\velx(x_{i})-\velx^{exact}(x_{i})]^{2}}{L}}\
,\ x_{i}\in[\Delta,1],\ \velx^{exact}=0.$ (26)

The simulation error at different $\IBMparam$ and $\IBMparam_{v}$ is compared
in Figure 14. It is clear that adding the second-order term will lead to lower
error when the diffusivity parameter $\IBMparam_{v}$ is chosen properly. For a
given penalization parameter and problem setting, an optimal $\IBMparam_{v}$
exists to provide the lowest error. However, a big gain is only seen for
$\IBMparam=1\times 10^{-3}$ and $\IBMparam=1\times 10^{-4}$, where the optimal
$\IBMparam_{v}$ will lead to a large reduction of the error. When the
penalization parameter is sufficiently small, e.g., $\IBMparam=1\times
10^{-5}$, the optimal performance requires a relatively larger $\IBMparam_{v}$
and the gain becomes marginal. In addition, the critical time step $\Delta t$
for different $\IBMparam_{v}$ is compared in Figure 15. This time step is the
largest $\Delta 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 $\IBMparam=1\times
10^{-4}$ and $\IBMparam=1\times 10^{-5}$ is almost the same, while the time
step can be about ten times larger ($3.1\times 10^{-5}$ against $3.1\times
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.

(a)

(b)

Fig. 16: Relation between optimal $\IBMparam_{v}$ (defined as
$\IBMparam_{v,opt}$) and $\IBMparam$ for various cases with different solid
ratio $r$ (one solid element) and polynomial order $P$. a) $\IBMparam_{v,opt}$
versus $\IBMparam$. b) Scaling relationship between $\IBMparam_{v,opt}$ and
$\IBMparam$.

Finally, the relationship between optimal $\IBMparam_{v}$ (defined as
$\IBMparam_{v,opt}$) and $\IBMparam$ is explored. From Figure 14, the optimal
$\IBMparam_{v}$ (indicated by black symbols) increases as $\IBMparam$ 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 $\IBMparam$
including $1\times 10^{-3}$, $1\times 10^{-4}$ and $1\times 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 $\IBMparam_{v}$ are performed to numerically find
$\IBMparam_{v,opt}$. From Figure 16a, a nearly linear relation is seen between
$\IBMparam_{v,opt}$ and $\IBMparam$. For a given $r$ and $\IBMparam$,
$\IBMparam_{v,opt}$ remains approximately constant for polynomial order $2$
and $3$, while it reduces by half for polynomial order $4$. Therefore, scaling
relationships of $\IBMparam\IBMparam_{v,opt}/(r^{2})=0.09$ and
$\IBMparam\IBMparam_{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
$\IBMparam_{v,opt}$ based on $\IBMparam$ and $P$.

###  4.3 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]:

$\bm{S}_{1}=\frac{\chi}{\eta}\times\begin{pmatrix}\rho u_{s}-\rho u\\\ \rho
v_{s}-\rho
v\end{pmatrix}+\rho\chi\IBMparam_{v}\times\begin{pmatrix}\frac{\partial^{2}u}{\partial
x^{2}}+\frac{\partial^{2}u}{\partial x\partial y}\\\
\frac{\partial^{2}v}{\partial x^{2}}+\frac{\partial^{2}v}{\partial y\partial
x}\end{pmatrix},$ (27)

where $\rho$, $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]

$\bm{S}_{2}=\frac{\chi}{\eta}\times\begin{pmatrix}\frac{\rho}{2}(u_{s}^{2}+v_{s}^{2})-\frac{\rho}{2}(u^{2}+v^{2})\end{pmatrix}.$
(28)

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

$r=\sqrt{\frac{S_{solid}}{S_{domain}}},$ (29)

where $S_{solid}$ and $S_{domain}$ are areas of the solid region and the whole
computational domain, respectively.

(a)

(b)

Fig. 17: Time evolution of aerodynamic coefficients for flow past a cylinder
at $Re=40$. Case 1: $\IBMparam=2\times 10^{-4}$, $\IBMparam_{v}=0$. Case 2:
$\IBMparam=1\times 10^{-3}$, $\IBMparam_{v}=1.5\times 10^{-2}$ (predicted by
the relation $\IBMparam\IBMparam_{v,opt}/(r^{2})=0.09$). Case 3:
$\IBMparam=1\times 10^{-3}$, $\IBMparam_{v}=1\times 10^{-2}$. Case 4:
$\IBMparam=1\times 10^{-3}$, $\IBMparam_{v}=2\times 10^{-2}$. a) Lift
coefficient. b) Drag coefficient.

A structured Cartesian mesh with local refinement region near the wall is
considered for the present numerical experiment. The size of rectangular
computational domain is $x\in[-30D,50D]$ and $y\in[-30D,30D]$, where $D$ is
the diameter of the cylinder. In the square region $x\in[-0.6D,0.6D]$ and
$y\in[-D,D]$, uniform grid with size $0.03D$ is chosen. This results in the
total number of elements $89\times 86$. Explicit time integration based on the
low-storage 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 $\Delta t$ is varied according different combinations
of parameters. At first, the volume penalization is applied, where we set
$\eta=\Delta t$ from the traditional guideline [47, 51], which gives
$\eta=\Delta t=2\times 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 $\IBMparam=1\times
10^{-3}$. From the relationship obtained in the previous section for $P=2$,
i.e., $\IBMparam\IBMparam_{v,opt}/(r^{2})=0.09$, we can obtain
$\IBMparam_{v,opt}\approx 1.5\times 10^{-2}$, where $r\approx 0.0128$ is
obtained from Equation 29. Therefore, three additional simulations with
$\IBMparam=1\times 10^{-3}$ and $\IBMparam_{v}=1.5\times 10^{-2}$, $1.0\times
10^{-2}$ and $2.0\times 10^{-2}$ are included for comparison, where the same
time step $\Delta t=2\times 10^{-4}$ is used for the first two cases, while
$\Delta t=1.5\times 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 $\IBMparam$ and $\IBMparam_{v}$, where $\IBMparam_{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].

Fig. 18: Predicted absolute value of lift coefficient across different
diffusivity coefficients. Black symbols represent the optimal $\IBMparam_{v}$
obtained from the relationship $\IBMparam\IBMparam_{v,opt}/(r^{2})=0.09$.

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 $\IBMparam_{v}$ at three penalization parameters,
including $1\times 10^{-2}$, $1\times 10^{-3}$ and $2\times 10^{-4}$ are
considered. It can be seen that for small $\IBMparam$, there exists an optimal
$\IBMparam_{v}$ that gives the best lift coefficient. When $\IBMparam=1\times
10^{-3}$, the proposed relationship accurately predicted the
$\IBMparam_{v,opt}$, as validated in Figure 17. At large $\IBMparam$ where the
solution is not sufficiently penalized, adding the second-order term will lead
to a slightly larger error, while the predicted $\IBMparam_{v,opt}$ still
gives reasonable result since a relatively small $\IBMparam_{v}=1.5\times
10^{-3}$ is obtained. When $\IBMparam=2\times 10^{-4}$, although the optimal
value is not well predicted, $\IBMparam_{v,opt}=7\times 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.

##  5 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\Delta t<\IBMparam_{critical}<0.5\Delta t$. This analysis
provides guidelines to the selection of the penalization parameter, where a
smaller $\IBMparam$ 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

###  A.1 Preliminaries

As the first step of space discretization for a one-dimensional problem, the
computational domain $\Domain$ is divided into $N$ nonoverlapping cells,
defined as $\Domain_{n}=\left\\{x|x_{n}<x<x_{n+1}\right\\}$ 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 $\Domain_{ref}=\left\\{r|-1<r<1\right\\}$. The mapping between
physical space and reference space can be defined as the function $\Gamma$

$r=\Gamma_{n}(x)=2\left(\frac{x-x_{n}}{x_{n+1}-x_{n}}\right)-1,$ (30)

and its inversion

$x=\Gamma_{n}^{-1}(r)=\left(\frac{1-r}{2}\right)x_{n}+\left(\frac{1+r}{2}\right)x_{n+1}.$
(31)

This transformation leads to the transformed conservation law:

$\frac{\partial\tilde{u}}{\partial t}+\frac{\partial\tilde{f}}{\partial
r}=\frac{\partial u(\Gamma_{n}^{-1}(r),t)}{\partial
t}+\frac{1}{J_{n}}\frac{f(\Gamma_{n}^{-1}(r),t)}{\partial r}=0,$ (32)

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 $\hat{u}_{j}$ of a Legendre basis, or a
vector of nodal coefficients $\tilde{u}_{j}$ defined at solution points. Here
we retain the former set of basis.

###  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 ($\tilde{u}(r_{i}),i=0,1,...,P$) with order $P$ for both solution
and flux is defined at $P+1$ solution points:

$\tilde{u}_{n}(r)=\sum_{i=0}^{P}l_{i}(r)\tilde{u}_{n}(r_{i})=\boldsymbol{l}^{T}\tilde{\boldsymbol{u}}_{n},$
(33)
$\tilde{f}^{D}_{n}(r)=\sum_{i=0}^{P}l_{i}(r)\tilde{f}^{D}_{n}(r_{i})=\sum_{i=0}^{P}l_{i}(r)\advcoef\tilde{u}_{n}(r_{i})=\advcoef\boldsymbol{l}^{T}\tilde{\boldsymbol{u}}_{n},$
(34)

where $\tilde{f}^{D}$ denotes the transformed discontinuous flux computed from
the value at solution points for the $n$th cell. $\boldsymbol{l}$ is a
Lagrange interpolation operator and $l_{i}(r)$ denotes the Lagrange polynomial
defined at a solution point $r_{i}$. This polynomial is defined as the
interpolating polynomial through $P+1$ solution points, which takes value $1$
at $r_{i}$ and zero at all other solution points:

$l_{i}(r)=\prod^{P}_{j=0,j\neq i}\left(\frac{r-r_{j}}{r_{i}-r_{j}}\right)$
(35)

In this work, Gauss quadrature points are used as the solution points, while
another popular choice is based on Gauss-Lobatto quadrature rule (including
element end-points). With this representation, the gradient of transformed
discontinuous flux at all solution points is obtained as

$\frac{\partial\tilde{f}^{D}_{n}(r)}{\partial
r}=\sum_{i=0}^{P}\tilde{f}^{D}_{n}(r_{i})\frac{\partial l_{i}(r)}{\partial
r}=\sum_{i=0}^{P}\advcoef\tilde{u}_{n}(r_{i})\frac{\partial l_{i}(r)}{\partial
r}=\advcoef\boldsymbol{d}^{T}\tilde{\boldsymbol{u}}_{n},$ (36)

where $\boldsymbol{d}$ is the gradient operator to obtain the discontinuous
gradient at position $r$, which can be generalized to gradient matrix
$\boldsymbol{D}$ at all solution points $\boldsymbol{D}[i,j]=\frac{\partial
l_{j}(r_{i})}{\partial r}$.

In the third stage of FR, the transformed numerical flux $\tilde{f}^{I}$ at
all element interfaces (either ends of the cell $\Domain_{r}$, $r=\pm 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 $n$th cell, we have the left and
right interpolated solutions
$\tilde{u}^{L}_{n}=\boldsymbol{l}^{T}(-1)\tilde{\boldsymbol{u}}_{n}$ and
$\tilde{u}^{R}_{n}=\boldsymbol{l}^{T}(1)\tilde{\boldsymbol{u}}_{n}$, as well
as the interpolated fluxes $\tilde{f}^{D,L}_{n}=\advcoef\tilde{u}^{L}_{n}$ and
$\tilde{f}^{D,R}_{n}=\advcoef\tilde{u}^{R}_{n}$. The exact method for
numerical flux calculation
$\tilde{f}^{I}=\tilde{f}^{I}(\tilde{u}_{L},\tilde{u}_{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

$\tilde{f}^{I,L}_{n}=\frac{\tilde{f}^{D,L}_{n}+\tilde{f}^{D,R}_{n-1}}{2}-\frac{1}{2}\lambda\left|\advcoef\right|(\tilde{u}^{L}_{n}-\tilde{u}^{R}_{n-1})=\frac{1}{2}(\advcoef+\lambda\left|\advcoef\right|)\boldsymbol{l}^{T}(1)\tilde{\boldsymbol{u}}_{n-1}+\frac{1}{2}(\advcoef-\lambda\left|\advcoef\right|)\boldsymbol{l}^{T}(-1)\tilde{\boldsymbol{u}}_{n},$
(37)

$\tilde{f}^{I,R}_{n}=\frac{\tilde{f}^{D,R}_{n}+\tilde{f}^{D,L}_{n+1}}{2}-\frac{1}{2}\lambda\left|\advcoef\right|(\tilde{u}^{R}_{n}-\tilde{u}^{L}_{n+1})=\frac{1}{2}(\advcoef+\lambda\left|\advcoef\right|)\boldsymbol{l}^{T}(1)\tilde{\boldsymbol{u}}_{n}+\frac{1}{2}(\advcoef-\lambda\left|\advcoef\right|)\boldsymbol{l}^{T}(-1)\tilde{\boldsymbol{u}}_{n+1},$
(38)

where $\lambda$ is an upwinding parameter that controls the numerical flux.
$\lambda=1$ defines the upwind flux and $\lambda=0$ defines the central flux.
The upwind flux is considered here. The numerical interface fluxes for both
boundaries of the $n$th cell are defined as $\tilde{f}^{I,L}_{n}$ and
$\tilde{f}^{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 flux
$\tilde{f}^{C}$ based on a polynomial of degree $P+1$ is defined to make the
sum between $\tilde{f}^{C}$ and the interpolated discontinuous flux equal to
the transformed numerical interface flux at $r\pm 1$:

$\tilde{f}^{C}_{n}(-1)=\tilde{f}^{I,L}_{n}-\tilde{f}^{D,L}_{n}\ ,\
\tilde{f}^{C}_{n}(1)=\tilde{f}^{I,R}_{n}-\tilde{f}^{D,R}_{n}.$ (39)

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,$ (40)

with a symmetry consideration

$g^{L}(r)=g^{R}(-r).$ (41)

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

$\tilde{f}^{C}_{n}(r)=(\tilde{f}^{I,L}_{n}-\tilde{f}^{D,L}_{n})g^{L}(r)+(\tilde{f}^{I,R}_{n}-\tilde{f}^{D,R}_{n})g^{R}(r).$
(42)

Therefore the total transformed flux $\tilde{f}$ with a polynomial of degree
$P+1$ is formed from the discontinuous and correction flux as follows

$\tilde{f}_{n}=\tilde{f}^{D}_{n}+\tilde{f}^{C}_{n}.$ (43)

In the final stage of FR process, the divergence of the transformed total flux
$\frac{\partial\tilde{f}}{\partial{r}}$ at each solution point is calculated
as follows

$\frac{\partial\tilde{f}_{n}}{\partial{r}}(r_{i})=\sum_{j=0}^{P}\tilde{f}_{n}^{D}(r_{j})\frac{\partial
l_{j}}{\partial
r}(r_{i})+(\tilde{f}^{I,L}_{n}-\tilde{f}^{D,L}_{n})\frac{\partial
g^{L}}{\partial
r}(r_{i})+(\tilde{f}^{I,R}_{n}-\tilde{f}^{D,R}_{n})\frac{\partial
g^{R}}{\partial r}(r_{i}).$ (44)

This results in the semi-discrete formulation:

$\frac{du_{n}}{dt}=-\frac{1}{J_{n}}\frac{\partial\tilde{f}_{n}}{\partial{r}}.$
(45)

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.

###  A.3 Local matrix formulation

The eigensolution analysis relies on the eigendecomposition of the
discretization operator in matrix form. Consider uniform element space $h$ for
the computational domain, we can obtain the local matrix formulation for the
$n$th element as follows,

$\frac{d\boldsymbol{u}_{n}}{dt}=\boldsymbol{L}\boldsymbol{u}_{n-1}+\boldsymbol{C}\boldsymbol{u}_{n}+\boldsymbol{R}\boldsymbol{u}_{n+1},$
(46)

where the left, middle and right operators for the current element
$\boldsymbol{u}_{n}$ are defined as $\boldsymbol{L}$, $\boldsymbol{C}$, and
$\boldsymbol{R}$, respectively. These matrices can be derived based on Eq.
(36), Eq. (37), Eq. (38) and Eq. (44):

$\boldsymbol{L}=-\frac{1}{h}\boldsymbol{g}^{L}_{r}\left[(\advcoef+\lambda\left|\advcoef\right|)\boldsymbol{l}^{T}(1)\right],$
(47)

$\boldsymbol{C}=-\frac{2}{h}\left[c\boldsymbol{D}-\frac{1}{2}\boldsymbol{g}^{L}_{r}(\advcoef+\lambda\left|\advcoef\right|)\boldsymbol{l}^{T}(-1)+\frac{1}{2}\boldsymbol{g}^{R}_{r}(\lambda\left|\advcoef\right|-\advcoef)\boldsymbol{l}^{T}(1)\right],$
(48)

$\boldsymbol{R}=-\frac{1}{h}\boldsymbol{g}^{R}_{r}\left[(c-\lambda\left|\advcoef\right|)\boldsymbol{l}^{T}(-1)\right],$
(49)

where the gradient vector of correction function at all solution points is
defined as

$\boldsymbol{g}^{dir}_{r}=\left[\frac{\partial g^{dir}}{\partial
r}(r_{0}),\frac{\partial g^{dir}}{\partial r}(r_{1}),...,\frac{\partial
g^{dir}}{\partial r}(r_{P})\right]^{T}\ ,\ dir=L\ or\ R.$ (50)

## References

  * Mahesh [1998] K. Mahesh,  A family of high order finite difference schemes with good spectral resolution,  Journal of Computational Physics 145 (1998) 332–358.
  * Hirsch [2007] C. Hirsch, Numerical computation of internal and external flows: The fundamentals of computational fluid dynamics, Elsevier, 2007.
  * LeVeque et al. [2002] R. J. LeVeque, et al., Finite volume methods for hyperbolic problems, volume 31, Cambridge university press, 2002.
  * Hughes [2012] T. J. Hughes, The finite element method: linear static and dynamic finite element analysis, Courier Corporation, 2012.
  * Chávez-Modena et al. [2018] M. Chávez-Modena, E. Ferrer, G. Rubio,  Improving the stability of multiple-relaxation lattice boltzmann methods with central moments,  Computers & Fluids 172 (2018) 397–409.
  * Chávez-Modena et al. [2020] M. Chávez-Modena, A. Martínez-Cava, G. Rubio, E. Ferrer,  Optimizing free parameters in the d3q19 multiple-relaxation lattice boltzmann methods to simulate under-resolved turbulent flows,  Journal of Computational Science 45 (2020) 101170.
  * Hu et al. [1999] F. Q. Hu, M. Hussaini, P. Rasetarinera,  An analysis of the discontinuous galerkin method for wave propagation problems,  Journal of Computational Physics 151 (1999) 921–946.
  * Van den Abeele et al. [2007] K. Van den Abeele, T. Broeckhoven, C. Lacor,  Dispersion and dissipation properties of the 1d spectral volume method and application to a p-multigrid algorithm,  Journal of Computational Physics 224 (2007) 616–636.
  * Gassner and Kopriva [2011] G. Gassner, D. A. Kopriva,  A comparison of the dispersion and dissipation errors of gauss and gauss–lobatto discontinuous galerkin spectral element methods,  SIAM Journal on Scientific Computing 33 (2011) 2560–2579.
  * Vincent et al. [2011] P. E. Vincent, P. Castonguay, A. Jameson,  Insights from von neumann analysis of high-order flux reconstruction schemes,  Journal of Computational Physics 230 (2011) 8134–8154.
  * Moura et al. [2015] R. C. Moura, S. J. Sherwin, J. Peiró,  Linear dispersion–diffusion analysis and its application to under-resolved turbulence simulations using discontinuous galerkin spectral/hp methods,  Journal of Computational Physics 298 (2015) 695–710.
  * Alhawwary and Wang [2018] M. Alhawwary, Z. J. Wang,  Fourier analysis and evaluation of dg, fd and compact difference methods for conservation laws,  Journal of Computational Physics 373 (2018) 835–862.
  * Manzanero et al. [2018] J. Manzanero, G. Rubio, E. Ferrer, E. Valero,  Dispersion-dissipation analysis for advection problems with nonconstant coefficients: Applications to discontinuous galerkin formulations,  SIAM Journal on Scientific Computing 40 (2018) A747–A768.
  * Moura et al. [2016] R. C. Moura, S. J. Sherwin, J. Peiró,  Eigensolution analysis of spectral/hp continuous galerkin approximations to advection–diffusion problems: Insights into spectral vanishing viscosity,  Journal of Computational Physics 307 (2016) 401–422.
  * Manzanero et al. [2020] J. Manzanero, E. Ferrer, G. Rubio, E. Valero,  Design of a smagorinsky spectral vanishing viscosity turbulence model for discontinuous galerkin methods,  Computers & Fluids 200 (2020) 104440.
  * Solán-Fustero et al. [2021] P. Solán-Fustero, A. Navas-Montilla, E. Ferrer, J. Manzanero, P. García-Navarro,  Application of approximate dispersion-diffusion analyses to under-resolved burgers turbulence using high resolution weno and uwc schemes,  Journal of Computational Physics 435 (2021) 110246.
  * Trojak et al. [2020] W. Trojak, R. Watson, A. Scillitoe, P. G. Tucker,  Effect of mesh quality on flux reconstruction in multi-dimensions,  Journal of Scientific Computing 82 (2020) 1–36.
  * Mengaldo et al. [2018a] G. Mengaldo, R. Moura, B. Giralda, J. Peiró, S. Sherwin,  Spatial eigensolution analysis of discontinuous galerkin schemes with practical insights for under-resolved computations and implicit les,  Computers & Fluids 169 (2018a) 349–364.
  * Mengaldo et al. [2018b] G. Mengaldo, D. De Grazia, R. C. Moura, S. J. Sherwin,  Spatial eigensolution analysis of energy-stable flux reconstruction schemes and influence of the numerical flux on accuracy and robustness,  Journal of Computational Physics 358 (2018b) 1–20.
  * Kou et al. [2021] J. Kou, S. Le Clainche, E. Ferrer,  Data-driven eigensolution analysis based on a spatio-temporal koopman decomposition, with applications to high order methods,  (Under Review) (2021).
  * Fernandez et al. [2019] P. Fernandez, R. C. Moura, G. Mengaldo, J. Peraire,  Non-modal analysis of spectral element methods: Towards accurate and robust large-eddy simulations,  Computer Methods in Applied Mechanics and Engineering 346 (2019) 43–62.
  * Yang et al. [2013] H. Yang, F. Li, J. Qiu,  Dispersion and dissipation errors of two fully discrete discontinuous galerkin methods,  Journal of Scientific Computing 55 (2013) 552–574.
  * Vermeire and Vincent [2017] B. Vermeire, P. Vincent,  On the behaviour of fully-discrete flux reconstruction schemes,  Computer Methods in Applied Mechanics and Engineering 315 (2017) 1053–1079.
  * Mittal and Iaccarino [2005] R. Mittal, G. Iaccarino,  Immersed boundary methods,  Annu. Rev. Fluid Mech. 37 (2005) 239–261.
  * Sotiropoulos and Yang [2014] F. Sotiropoulos, X. Yang,  Immersed boundary methods for simulating fluid-structure interaction,  Progress in Aerospace Sciences 65 (2014) 1–21.
  * Griffith and Patankar [2020] B. E. Griffith, N. A. Patankar,  Immersed methods for fluid–structure interaction,  Annual Review of Fluid Mechanics 52 (2020) 421–448.
  * Peskin [1972] C. S. Peskin,  Flow patterns around heart valves: a numerical method,  Journal of computational physics 10 (1972) 252–271.
  * Taira and Colonius [2007] K. Taira, T. Colonius,  The immersed boundary method: a projection approach,  Journal of Computational Physics 225 (2007) 2118–2137.
  * Wu and Shu [2009] J. Wu, C. Shu,  Implicit velocity correction-based immersed boundary-lattice boltzmann method and its applications,  Journal of Computational Physics 228 (2009) 1963–1979.
  * Iaccarino and Verzicco [2003] G. Iaccarino, R. Verzicco,  Immersed boundary technique for turbulent flow simulations,  Appl. Mech. Rev. 56 (2003) 331–347.
  * Tamaki et al. [2017] Y. Tamaki, M. Harada, T. Imamura,  Near-wall modification of spalart–allmaras turbulence model for immersed boundary method,  AIAA journal 55 (2017) 3027–3039.
  * Huang et al. [2007] W.-X. Huang, S. J. Shin, H. J. Sung,  Simulation of flexible filaments in a uniform flow by the immersed boundary method,  Journal of computational physics 226 (2007) 2206–2228.
  * Yang and Stern [2012] J. Yang, F. Stern,  A simple and efficient direct forcing immersed boundary framework for fluid–structure interactions,  Journal of Computational Physics 231 (2012) 5029–5061.
  * Tian et al. [2014] F.-B. Tian, H. Dai, H. Luo, J. F. Doyle, B. Rousseau,  Fluid–structure interaction involving large deformations: 3d simulations and applications to biological systems,  Journal of computational physics 258 (2014) 451–469.
  * Seo and Mittal [2011] J. H. Seo, R. Mittal,  A high-order immersed boundary method for acoustic wave scattering and low-mach number flow-induced sound in complex geometries,  Journal of computational physics 230 (2011) 1000–1019.
  * Sun et al. [2012] X. Sun, Y. Jiang, A. Liang, X. Jing,  An immersed boundary computational model for acoustic scattering problems with complex geometries,  The Journal of the Acoustical Society of America 132 (2012) 3190–3199.
  * Wang et al. [2017] L. Wang, G. M. Currao, F. Han, A. J. Neely, J. Young, F.-B. Tian,  An immersed boundary method for fluid–structure interaction with compressible multiphase flows,  Journal of Computational Physics 346 (2017) 131–151.
  * O’Brien and Bussmann [2018] A. O’Brien, M. Bussmann,  A volume-of-fluid ghost-cell immersed boundary method for multiphase flows with contact line dynamics,  Computers & Fluids 165 (2018) 43–53.
  * Ye et al. [1999] T. Ye, R. Mittal, H. Udaykumar, W. Shyy,  An accurate cartesian grid method for viscous incompressible flows with complex immersed boundaries,  Journal of computational physics 156 (1999) 209–240.
  * Udaykumar et al. [2001] H. Udaykumar, R. Mittal, P. Rampunggoon, A. Khanna,  A sharp interface cartesian grid method for simulating flows with complex moving boundaries,  Journal of computational physics 174 (2001) 345–380.
  * Goldstein et al. [1993] D. Goldstein, R. Handler, L. Sirovich,  Modeling a no-slip flow boundary with an external force field,  Journal of computational physics 105 (1993) 354–366.
  * Angot et al. [1999] P. Angot, C.-H. Bruneau, P. Fabrie,  A penalization method to take into account obstacles in incompressible viscous flows,  Numerische Mathematik 81 (1999) 497–520.
  * Fadlun et al. [2000] E. Fadlun, R. Verzicco, P. Orlandi, J. Mohd-Yusof,  Combined immersed-boundary finite-difference methods for three-dimensional complex flow simulations,  Journal of computational physics 161 (2000) 35–60.
  * Luo et al. [2012] H. Luo, H. Dai, P. J. F. de Sousa, B. Yin,  On the numerical oscillation of the direct-forcing immersed-boundary method for moving boundaries,  Computers & Fluids 56 (2012) 61–76.
  * Brown-Dymkoski et al. [2014] E. Brown-Dymkoski, N. Kasimov, O. V. Vasilyev,  A characteristic based volume penalization method for general evolution problems applied to compressible viscous flows,  Journal of Computational Physics 262 (2014) 344–357.
  * Abgrall et al. [2014] R. Abgrall, H. Beaugendre, C. Dobrzynski,  An immersed boundary method using unstructured anisotropic mesh adaptation combined with level-sets and penalization techniques,  Journal of Computational Physics 257 (2014) 83–101.
  * Schneider [2015] K. Schneider,  Immersed boundary methods for numerical simulation of confined fluid and plasma turbulence in complex geometries: a review,  arXiv preprint arXiv:1508.04593 (2015).
  * Kadoch et al. [2012] B. Kadoch, D. Kolomenskiy, P. Angot, K. Schneider,  A volume penalization method for incompressible flows and scalar advection–diffusion with moving obstacles,  Journal of Computational Physics 231 (2012) 4365–4383.
  * Kolomenskiy and Schneider [2009] D. Kolomenskiy, K. Schneider,  A fourier spectral method for the navier–stokes equations with volume penalization for moving solid obstacles,  Journal of Computational Physics 228 (2009) 5687–5709.
  * Horgue et al. [2014] P. Horgue, M. Prat, M. Quintard,  A penalization technique applied to the “volume-of-fluid” method: Wettability condition on immersed boundaries,  Computers & Fluids 100 (2014) 255–266.
  * Engels et al. [2015] T. Engels, D. Kolomenskiy, K. Schneider, J. Sesterhenn,  Numerical simulation of fluid–structure interaction with the volume penalization method,  Journal of Computational Physics 281 (2015) 96–115.
  * Cui et al. [2018] X. Cui, X. Yao, Z. Wang, M. Liu,  A coupled volume penalization-thermal lattice boltzmann method for thermal flows,  International Journal of Heat and Mass Transfer 127 (2018) 253–266.
  * Kou et al. [2021a] J. Kou, S. Joshi, A. Hurtado-de Mendoza, K. Puri, C. Hirsch, E. Ferrer,  High-order flux reconstruction based on immersed boundary method,  in: 14th WCCM-ECCOMAS Congress 2020, volume 700, 2021a.
  * Kou et al. [2021b] J. Kou, S. Joshi, A. Hurtado-de Mendoza, K. Puri, C. Hirsch, E. Ferrer,  Immersed boundary method for high-order flux reconstruction based on volume penalization,  (Under Review) (2021b).
  * Tu and Peskin [1992] C. Tu, C. S. Peskin,  Stability and instability in the computation of flows with moving immersed boundaries: a comparison of three methods,  SIAM Journal on Scientific and Statistical Computing 13 (1992) 1361–1376.
  * Stockie and Wetton [1999] J. M. Stockie, B. R. Wetton,  Analysis of stiffness in the immersed boundary method and implications for time-stepping schemes,  Journal of Computational Physics 154 (1999) 41–64.
  * Gong et al. [2008] Z. Gong, H. Huang, C. Lu,  Stability analysis of the immersed boundary method for a two-dimensional membrane with bending rigidity,  Commun. Comput. Phys. 3 (2008) 704–723.
  * Kolomenskiy et al. [2014] D. Kolomenskiy, K. Schneider, et al.,  Approximation of the laplace and stokes operators with dirichlet boundary conditions through volume penalization: a spectral viewpoint,  Numerische Mathematik 128 (2014) 301–338.
  * Kolomenskiy et al. [2015] D. Kolomenskiy, K. Schneider, et al.,  Analysis and discretization of the volume penalized laplace operator with neumann boundary conditions,  Applied Numerical Mathematics 95 (2015) 238–249.
  * Lee et al. [2017] S. Lee, D. Y. Kwak, I. Sim,  Immersed finite element method for eigenvalue problem,  Journal of Computational and Applied Mathematics 313 (2017) 410–426.
  * Zhou and Balachandar [2021] K. Zhou, S. Balachandar,  An analysis of the spatio-temporal resolution of the immersed boundary method with direct forcing,  Journal of Computational Physics 424 (2021) 109862.
  * Wang et al. [2013] Z. J. Wang, K. Fidkowski, R. Abgrall, F. Bassi, D. Caraeni, A. Cary, H. Deconinck, R. Hartmann, K. Hillewaert, H. T. Huynh, et al.,  High-order cfd methods: current status and perspective,  International Journal for Numerical Methods in Fluids 72 (2013) 811–845.
  * Hesthaven and Warburton [2007] J. S. Hesthaven, T. Warburton, Nodal discontinuous Galerkin methods: algorithms, analysis, and applications, Springer Science & Business Media, 2007\.
  * Karniadakis and Sherwin [2013] G. Karniadakis, S. Sherwin, Spectral/hp element methods for computational fluid dynamics, Oxford University Press, 2013\.
  * Kopriva and Kolias [1996] D. A. Kopriva, J. H. Kolias,  A conservative staggered-grid chebyshev multidomain method for compressible flows,  Journal of computational physics 125 (1996) 244–261.
  * Liu et al. [2006] Y. Liu, M. Vinokur, Z. J. Wang,  Spectral difference method for unstructured grids i: basic formulation,  Journal of Computational Physics 216 (2006) 780–801.
  * Huynh [2007] H. T. Huynh,  A flux reconstruction approach to high-order schemes including discontinuous galerkin methods,  in: 18th AIAA Computational Fluid Dynamics Conference, 2007, p. 4079.
  * Wang and Gao [2009] Z. J. Wang, H. Gao,  A unifying lifting collocation penalty formulation including the discontinuous galerkin, spectral volume/difference methods for conservation laws on mixed grids,  Journal of Computational Physics 228 (2009) 8161–8186.
  * Vincent et al. [2011] P. E. Vincent, P. Castonguay, A. Jameson,  A new class of high-order energy stable flux reconstruction schemes,  Journal of Scientific Computing 47 (2011) 50–72.
  * Sakurai et al. [2019] T. Sakurai, K. Yoshimatsu, N. Okamoto, K. Schneider,  Volume penalization for inhomogeneous neumann boundary conditions modeling scalar flux in complicated geometry,  Journal of Computational Physics 390 (2019) 452–469.
  * Sengupta et al. [2020] S. Sengupta, T. K. Sengupta, J. K. Puttam, K. S. Vajjala,  Global spectral analysis for convection-diffusion-reaction equation in one and two-dimensions: Effects of numerical anti-diffusion and dispersion,  Journal of Computational Physics 408 (2020) 109310.
  * Carbou and Fabrie [2003] G. Carbou, P. Fabrie,  Boundary layer for a penalization method for viscous incompressible flow,  Advances in Differential equations 8 (2003) 1453–1480.
  * Gautier et al. [2014] R. Gautier, S. Laizet, E. Lamballais,  A dns study of jet control with microjets using an immersed boundary method,  International Journal of Computational Fluid Dynamics 28 (2014) 393–410.
  * Stein et al. [2016] D. B. Stein, R. D. Guy, B. Thomases,  Immersed boundary smooth extension: a high-order method for solving pde on arbitrary smooth domains using fourier spectral methods,  Journal of Computational Physics 304 (2016) 252–274.
  * Jause Labert [2012] C. Jause Labert, Simulation numerique d’ecoulements turbulents en rotation, confinement et forcage a l’aide d’une methode de penalisation, Ph.D. thesis, Ecully, Ecole centrale de Lyon, 2012\.
  * He et al. [2020] X. He, D. Yang, C. Qiu,  Dispersion-dissipation analysis of triangular numerical-flux-based discontinuous galerkin method for elastic wave equations,  Journal of Computational Physics 418 (2020) 109630.
  * Cockburn and Shu [1989] B. Cockburn, C.-W. Shu,  Tvb runge-kutta local projection discontinuous galerkin finite element method for conservation laws. ii. general framework,  Mathematics of computation 52 (1989) 411–435.
  * Chalmers et al. [2014] N. Chalmers, L. Krivodonova, R. Qin,  Relaxing the cfl number of the discontinuous galerkin method,  SIAM Journal on Scientific Computing 36 (2014) A2047–A2075.
  * Cockburn and Shu [1998] B. Cockburn, C.-W. Shu,  The local discontinuous galerkin method for time-dependent convection-diffusion systems,  SIAM Journal on Numerical Analysis 35 (1998) 2440–2463.
  * Carpenter and Kennedy [1994] M. H. Carpenter, C. A. Kennedy,  Fourth-order 2n-storage runge-kutta schemes (1994).
  * de Tullio et al. [2007] M. D. de Tullio, P. De Palma, G. Iaccarino, G. Pascazio, M. Napolitano,  An immersed boundary method for compressible flows using local grid refinement,  Journal of Computational Physics 225 (2007) 2098–2117.
  * Williams et al. [2013] D. M. Williams, P. Castonguay, P. E. Vincent, A. Jameson,  Energy stable flux reconstruction schemes for advection–diffusion problems on triangles,  Journal of Computational Physics 250 (2013) 53–76.
  * Roe [1981] P. L. Roe,  Approximate riemann solvers, parameter vectors, and difference schemes,  Journal of computational physics 43 (1981) 357–372.