A modified equation analysis for immersed boundary methods based on volume penalization: applications to linear advection-diffusion and high-order discontinuous Galerkin schemes

The Immersed Boundary Method (IBM) is a popular numerical approach to impose boundary conditions without relying on body-fitted grids, thus reducing the costly effort of mesh generation. To obtain enhanced accuracy, IBM can be combined with high-order methods (e.g., discontinuous Galerkin). For this combination to be effective, an analysis of the numerical errors is essential. In this work, we apply, for the first time, a modified equation analysis to the combination of IBM (based on volume penalization) and high-order methods (based on nodal discontinuous Galerkin methods) to analyze a priori numerical errors and obtain practical guidelines on the selection of IBM parameters. The analysis is performed on a linear advection-diffusion equation with Dirichlet boundary conditions. Three ways to penalize the immerse boundary are considered, the first penalizes the solution inside the IBM region (classic approach), whilst the second and third penalize the first and second derivatives of the solution. We find optimal combinations of the penalization parameters, including the first and second penalizing derivatives, resulting in minimum errors. We validate the theoretical analysis with numerical experiments for one- and two-dimensional advection-diffusion equations.


Introduction
Despite successful applications of Immersed Boundary Method (IBM) in simulating complex flows [41,21] and fluid-structure interaction problems [53,23,17], understanding and controlling numerical errors in the IBM approach remains a challenge.IBM refers to a group of numerical strategies that handle the boundary condition when the solid is immersed in the computational domain, avoiding body-fitted meshes and enabling the use of simple meshes (e.g., Cartesian or Octree).The IBM approach originates from the idea of Peskin [45], where singular forces represented by delta functions were positioned at solid boundaries to mimic the effect of physical boundaries.Since IBM reduces the complexity of mesh generation and handles moving boundaries efficiently, it has received a lot of attention over the past few decades.In general, IBM treatment can be achieved using the cut cell approach [63,58,15,54] or by introducing source terms, e.g., ghost cell [38], projection method [55], direct forcing [14,37] or volume penalization [3,24].Although the cut-cell approach shows better convergence properties, the extension to moving boundaries and the treatment of different types of cut-cells remain challenging.A more flexible approach is the IBM based on Volume Penalization (VP).The latter shows advantages in robustness, ease of implementation, and theoretical convergence estimates [5,3].
Volume penalization is a classic IBM treatment based on modeling the solid as a porous medium with low permeability [22,47].The method imposes the boundary condition by introducing a source term (or penalty term) to the computational nodes located inside the solid.This approach dates back to the work of Courant [11], where a penalty method was used to transform constrained optimization problems into a constraint-free problem.Volume penalization methods for the Navier-Stokes equations were first proposed by Arquis and Caltagirone [5] with a Brinkman-type penalization for the momentum equations.After that, Angot et al. [3] and Carbou and Fabrie [9] proved the convergence of volume penalization, showing that as the penalization parameter η approaches 0, the model error converges if non-slip boundary conditions are considered.Subsequently, the volume penalization was extended to allow Neumann boundary conditions [22,47] and Robin boundary conditions [46], as well as spatially varying Neumann and Robin boundary conditions [56].The volume penalization method was extended to compressible flows by Liu and Vasilyev [35], Brown-Dymkoski et al. [8], Abgrall et al. [2] and Abalakin et al. [1].This method has been applied to a variety of problems, including flapping wings [24], two-phase flows [20], aeroacoustics [26], fluid-structure interactions [13] and thermal flows [12].So far, IBM research for high-order methods has been relatively unexplored, with efforts focused on Poisson problems [32,33] and cut-cell approaches [15,44].In the context of volume penalization, we have recently extended this approach to high-order flux reconstruction schemes [28,29], and now to the high-order discontinuous Galerkin Spectral Element Method (DGSEM) in this work.
There have been several attempts to analyze the errors of the IBM approach.Bever and Leveque [7] analyzed the error of traditional IBM applied to one-dimensional problems, and highlighted the importance of choosing appropriate discrete delta functions to maintain optimal accuracy.Following a similar strategy, the immersed interface method [34], which modifies the finite difference scheme with a jump condition for the immersed boundary, was derived [31].Tornberg and Engquist [57] performed error analyzes of traditional IBM with regularization and found first-order convergence for the standard central difference scheme with smoothing discrete delta functions.This error analysis was then extended to Stokes flows by Mori [42], Chen et al. [10], and Liu and Mori [36] where error estimates for velocity and pressure were reported.Most analyses focus on the traditional IBM method, where the numerical property is based on the selection of the appropriate delta functions.Error analyzes for the direct forcing approach were also explored in [18,65], where the importance of maintaining smoothness in the solution and the choice of a suitable temporal and spatial resolution was highlighted.For these types of approach, the discretization error from space-time discretization and the modeling error from particular IBM treatment are coupled.In contrast, volume penalization has the advantage that the modeling error and the numerical error can be handled separately.The convergence of modeling errors was studied rigorously by Angot et al. [3] and Carbou and Fabrie [9].The modeling errors for the incompressible flow past a cylinder and a sphere were analyzed by Zhang and Zheng [64].Therefore, the main concern is the discretization error, which depends on the numerical scheme used, and where a detailed error analysis is lacking, especially in the context of high-order methods.
In the present work, we perform error analyses of the IBM based on combination of volume penalization and nodal DGSEM, and propose new penalties to cancel spatial errors to improve the accuracy of the solution.In particular, we use modified equation analysis [51,43], which has been extensively used to analyze the stability and accuracy of low-order numerical discretization, and to obtain high-order schemes [61,50].The relationship between the errors introduced by the IBM based on volume penalization and high-order schemes remains unclear and motivates this work.First, using a modified equation analysis for volume penalization using DGSEM, we determine the shape and relationship of the dominant errors (i.e.dissipative/dispersive character).Second, we design the volume penalization scheme by including additional penalty terms that cancel the undesired numerical errors.In recent work, we have attempted to damp the numerical errors that arise from the volume penalization approach, using second-order derivatives [30] or combining it with a frequency damping technique [27].These studies have tried to minimize errors, without explicitly considering the causes of such errors, and therefore can be considered a posteriori palliative treatment.In this work, we consider a different perspective and analyze the source of these errors.By doing so, we are able to cancel the errors at the source.This approach can be considered an a priori error control.Note that we limit our analysis to linear advection-diffusion equations.The results from our analysis can thus be extended to linear systems (or linearized version of nonlinear equations).Examples include acoustics [49] or stability analysis [52].
The article is organized as follows.In Section 2, the volume penalization method and the DGSEM technique are introduced for the governing equation.Next, Section 3 introduces the principal errors in a volume penalization approach to investigate the discretization error using the modified equation analysis in Section 4. The numerical results are shown in Section 5, to validate the conclusions of the analysis, where one and two-dimensional advectiondiffusion equations are investigated.Finally, conclusions are given in Section 6.

The governing equation
Let us introduce the problem by considering the following time-dependent 1D advection-diffusion equation for the transported solution u = u(x, t), where the flow parameters are constant: velocity field c and kinetic viscosity ν ≥ ν min > 0. The PDE is completed with the set of initial and boundary conditions, The transport problem (1) is discretized in space based on a high-order DG method; in time, a Runge-Kutta method; and some of the solution points would be penalized by additional source terms to impose the IBM conditions.

The volume penalization approach
Motivated by the characteristic-based VP [8] and the inclusion of local dissipation [30], we consider the governing equation with penalization terms for the solution (classic volume penalization) and additional first-order and secondorder penalization terms: The additional term in Equation (2a) helps to impose the IBM on a given region of the domain Ω = [0, L].In this work, we consider a boundary condition of homogeneous Dirichlet type, namely u s = u s (x, t) = 0 (that is, a nonslip wall).The other two parameters are penalized terms determined in the modified equation analysis section, where we focus only on the spatial errors of the discretization.The penalization parameters for variable, first and second order derivatives are η 1 , η 2 , and η 3 respectively; and to facilitate the analysis, a continuous mask function represented by a hyperbolic tangent function is defined as which distinguishes between the fluid, Ω f , and the solid, Ω s , regions such that The distance of any solution point from the boundary interface is defined as d = d(x, t).The width of the hyperbolic tangent function is defined as δ, which should be infinitely small to reduce the modeling error and approximate the sharp mask function.This smooth mask ensures that spatial derivatives can be calculated.Note that for standard volume penalization, the mask function is sharp, which is 1 in the solid and 0 in the fluid.Later in Section 5.1, we will show that both sharp and smooth masks result in similar behavior given a sufficiently small δ.

The VP-DG discrete equation
Equation (2a) is discretized using the DG spectral element method.We group the terms in Equation (2a) that leads to the penalized equation: where the second-order differential operator is represented by where c and ν are the VP velocity field and the VP viscosity, respectively.Note that here we consider the element to be either a fully solid or a fully fluid one.This implies that the solid boundary aligns with the element interface, which is a natural choice for the DG method when using the local r-refinement, e.g.[40].The domain Ω is divided into multiple subdomains named elements 2, . . ., K, as can be seen in Figure 1, and mapped to the reference interval ξ ∈ [−1, 1].The global solution is assumed to be approximated by a piecewise polynomial defined as the direct sum ⊕ of the K local polynomial solutions, also for the test function and the VP flux function.On each element, we describe the local solution by the Legendre orthogonal interpolating polynomial, which is written in the Lagrange form, We select Gauss-Lobatto (GL) points, as they are becoming very popular in newly energy-stable and entropy-conserving schemes ( [39,62]).Using GL points, the nodal (grid point) values become u k h,j (t) = u k h (ξ j , t), and l j (ξ) is the N -th order Lagrange interpolating polynomial, which satisfy l j (ξ i ) = δ ij , being δ ij the Kronecker delta.After obtaining weak forms of Equation (3a) and applying the Gaussian quadrature to the inner product in the reference interval (see Appendix A for details), the semi-discrete equation writes as follows: for k = 1, 2, . . ., K and j = 0, 1, . . ., N where is the VP-DG derivative, and the numerical source.The weights of the GL quadrature are {w i } N i=0 and l represents the derivative of the Lagrange polynomial.In the previous formula, we define the numerical flux as F and its numerical counterpart as U.These fluxes are given by: The weights f and g depend on c, ν, and some numerical parameters that determine the advective/diffusive flux scheme used; see Appendix A for details.

Errors in volume penalization
Rigorous proofs of the convergence of modeling errors have been provided in previous work [3,9], showing that the numerical error introduced from the penalization term can be controlled a priori [8].Analysis of volume penalization suggests that the two contributions to total error are modeling and discretization errors [13]: where u exact is the exact solution of the governing equation, u η and u num.
η are the exact and numerical solutions of the penalized equation, respectively, and • is the L p norm used to quantify the error.The modeling error depends on the penalization parameter [48]: This explains that the convergence of the solution to the exact solution requires the error norm to approach zero for small penalization parameter limit, i.e., According to Angot et al. [3] and Carbou and Fabrie [9], the volume penalization gives α = 1/2, indicating that the penalization error has a decay rate of O( √ η 1 ) for Dirichlet boundary conditions.For the Neumann boundary conditions, a decay rate of O(η 1 ) can be obtained [25].
The second part of the overall error is the discretization error, which refers to the error between the exact solution and the numerical solution of the penalized equation.Details are given in the next section.

The modified equation analysis
When we discretize the penalized equation (3a) numerically, we translate it into a semi-discrete system (7a) for each element.This scheme is an approximation of our original equation.A different view is that the discrete system is the solution of modified differential equations but with some extra terms.This equation is named the modified equation (or reduced PDE): and is obtained by expanding the solutions in Equation (7a) with a Taylor series around a point in the mesh x(ξ j ).We omit the superscript k.HOT is the high-order term due to the Taylor series.L h is the operator L applied at x(ξ j ) and s(u h ) is a function of u h that is the solution transported from the element Ω k at the interfaces; see Figure 1.The purpose of the modified equation is to obtain the truncation error, T E, which is defined as the difference between the original equation and the modified equation.With a consistent discretization and a stable numerical scheme, the discretization error or the T E term writes as follows: where ∆x is a geometric discretization parameter representative of the grid spacing, p the order of accuracy of the numerical scheme, and C 0 , C 1 , C 2 , . . .some constants that do not depend on the discretization.However, the discretization error is not only determined by the numerical scheme, but is also limited by the regularity of the solution, as pointed out by Schneider et al. [22,48].For high-order methods, with good regularity at the interface, the high-order convergence property can be recovered, that is, O(∆x N +1 ).This has been shown in the recent work of Kou et al. [29].Here, we further analyze discretization errors to control and reduce these errors and improve accuracy.
Consider that finite-volume/difference methods are local by nature.In contrast, DGSEM is local within the whole domain but global within the element.Due to the non-local character of DG within each element, the solution depends on every point at the GL mesh.To perform the analysis, we center the solution at the same point of the source for each component of the discrete equation.Let us simplify the analysis to three GL points (N = 2), which are located at ξ 0 = −1, ξ 1 = 0 and ξ 2 = 1 with weights w 0 = 1/3, w 1 = 4/3 and w 2 = 1/3 respectively.The Lagrange polynomials are and the VP-DG matrix, comming from Equation (7b).Due to the non-local character of DG within each element, the solution depends on every point at the GL mesh.To perform the analysis, we center the solution at the same point of the source for each component of the discrete equation.For example, the discrete equation for the and the numerical source, Expanding u k h,1 and u k h,2 around u k h,0 , we have where the numerical source now reads: Neither u k−1 h,2 nor u k+1 h,0 can be expanded using Taylor series due to the discontinuous nature of DG.The terms in brackets of Equation (19) simplify to and For the Taylor term of first order, c k 0 = Ж (1)k 0 ; and the second order term, Finally, we find the modified equation at the left-boundary element, where and Additionally, the original PDE at the left-boundary element is and, therefore, the truncation error at the left-boundary element becomes: We can proceed in a similar manner to obtain the modified equations and truncation errors for the inner point, j = 1, and the right-boundary point, j = 2.For j = 1, u k h,0 and u k h,2 are centered on u k h,1 ; for j = 2, u k h,0 and u k h,1 are centered on u k h,2 .Their formulae can be written using equations ( 23), but with differences in the reactive parameter, r k j , the coefficient Ж and the numerical source, S k j , for j = 0, 1, 2, see Table 1 and 2. The source of the DG, s k DG,j , arises from the discontinuous nature of the DG approach (discontinuous boundary values) and the selected diffusive scheme.Now suppose that an element, Ω k , belongs to a solid region, Ω s , then χ k j = 1 for j = 0, 1, 2, and the truncation error still remains inside the solid.If we want to eliminate all the error terms in T E k j for j = 0, 1, 2, we need to solve the system: The reaction parameter and the coefficient Ж in the modified equations for a three-point GL grid.
for j = 0, 1, 2 and m ≥ 3.However, the problem is given by Eq. ( 24) has an infinite number of equations and finite unknowns.In total, there are 10 unknowns, which are the 4 weights f and g, The solution to the system is as follows.
This solution will be referred to as the trivial solution of the problem.At this point, one may wonder if there is any other set, a nontrivial family, that cleans up almost all the errors within the solid region.To investigate this, a determined system should be formed.Ideal errors remaining within the solid region should be: for m ≥ 3 as a representation of high order.However, the investigation of nontrivial solutions did not meet the previous requirement; see more details in the Appendix B. Table 3 summarizes both the trivial and non-trivial solutions that have been found.The trivial solution is the condition for DGSEM to compensate (or kill) spatial truncation errors within the body region, but additional insight can be obtained.If we substitute the values for η 2 and η 3 in our penalized equation and isolate χ, we get the following: If we are in the solid region, χ = 1, then the physical contribution of the PDE is removed, so this region is modeled with only the reaction penalization term, and therefore only time integration methods will lead to errors within the solid.This result agrees with the use of a typical characteristic-based volume penalization approach [8,1], where the RHS vanishes to smooth out the errors in the solid region but without a theoretical explanation, which is provided here.Cancellation of particular terms can reduce the error inside the solid, thus improving the accuracy in the fluid region.
Table 3: Summary of family of solutions for VP-IBM DGSEM, the trivial solution is the last row.* means equivalent to a continuous Galerkin (CG) method.
Finally, let us note that the modified equation analysis, based on Taylor expansions, is more general than the one first anticipated for the considered advection-diffusion problem.Indeed, the Taylor analysis can be reframed into a more general Fourier analysis framework, as detailed in Appendix C, where we show that Fourier analysis links with the Taylor series.

Numerical results
In this section we introduce two numerical experiments to evaluate and validate the trivial solution derived from the modified equation analysis.The first group of cases is the one-dimensional advection-diffusion equation, where the influence of penalization parameters is studied in detail.The optimal parameters obtained from the numerical experiments and the analysis of the modified equations are then applied to the two-dimensional advection-diffusion equation.

One-dimensional advection-diffusion equation
We start from the one-dimensional advection equation, where a non-slip wall is placed in the middle of the computational domain.This problem has been formulated in previous works [30,27], which is illustrated in Figure 2. The advection speed is set to c = 1 and the computational domain is defined in x ∈ [−1, 1], discretized by K equispaced elements with mesh size ∆x.The solution points are selected according to the Gauss-Lobatto quadrature rule, which is consistent with the previous analysis.Periodic boundary conditions are imposed on both sides of the computational domain.An upwind flux for the advection term is selected.The solid region is defined as a non-slip wall, i.e., u s = 0.It lies in the middle of the computational domain and starts from x = 0, whose width is defined as ∆ s , leading to the solid region 0 ≤ x ≤ ∆ s .For consistency with the analysis of the modified equations, we consider ∆ s = ∆x, which means that the solid boundaries lie exactly at the interface between the elements (if we have an even number of elements).This allows us to define the solid ratio r = 1/K as the ratio between the solid region and the computational domain.As shown in Figure 2, the initial wavelike solution passes through the non-slip wall in the middle, and the damped solution moves to the right as time evolves.Since periodic boundary conditions are considered, the solution will eventually become 0. If the non-slip wall boundary condition is exactly imposed, the solutions coming out of the wall will be zero.However, in practice, this solution depends on the damping provided by the volume penalization approach, where a smaller η makes the solution closer to zero.Therefore, the accuracy of the IBM imposition can be evaluated by comparing the exact solution (zero) and the damped solution in both the flow and solid regions (e.g., within a short advection time 0 < x < t ).We first perform the numerical experiment of a linear advection equation with a wavelike initial condition.The initial condition is defined as a sinusoidal wave with wavenumber ω, which is nondimensionalized by the mesh size ∆x and the polynomial order N , defined as ω∆x/(N + 1).Furthermore, due to the existence of a solid wall, the actual fluid domain is shorter than the entire computational domain; therefore, the effective wavenumber in the fluid region is greater than ω.This effective wavenumber is rescaled by the solid ratio r to ω = ω/(1−r) [30].We consider a spatial discretization with K = 40 elements in the computational domain (∆x = 0.05).Based on this mesh, we set ∆ s = ∆x with r = 1/40 and choose N = 3 as a representative order for high-order methods.The initial condition with wavenumber ω∆x/(N + 1) = 0.3223 is considered, which lies in the resolved wavenumber region of the scheme.The time integration is based on the third-order Runge-Kutta scheme.To reduce the temporal error, a sufficiently small time step is set to ∆t = 10 −5 .The final time is set to 1.1 to obtain a sufficiently penalized solution in the right region of the computational domain.Different combinations of parameters (with and without the first-order term) are considered.To evaluate accuracy, the error (in the flow) is defined as the error in x ∈ [∆ s , 1] and the penalized value u s = 0. Defining the number of solution points inside the flow domain of interest as N p = (N + 1)K, we have the L 2 -norm of the error as and the L 2 -norm of the error in the solid is defined as First, a numerical study is performed to justify the equivalence of using sharp or smooth mask functions (given the small width δ for the smooth mask function).Due to the fact that GL points overlap at the element interface, we can manually modify the smooth mask for the interface point corresponding to the solid element being (tanh(1/δ) + 1)/2 and the one in the fluid being (tanh(−1/δ) + 1)/2.We run the simulation at different δ, η 1 = 10 −3 , and final time 1.1, and compare the mean squared error between the results from the sharp and the smooth mask function.This error is compared in Figure 4, where the results based on the sharp or smooth mask are almost identical at small δ, and the difference becomes dominant when δ is sufficiently large δ > 0.5.Similar results are maintained when δ < 0.5, which is sufficient to guarantee the equivalence of the analysis for smooth and sharp masks.As δ further increases, the mask becomes too smooth, resulting in additional modeling errors due to the wrong representation of the interface.
A comparison of the solution at the final time is shown in Figure 3. Four cases are tested, where the first three cases contain only the volume penalization term for the solution, while the first-order penalization term with η 2 = −1/c is added to the last case.The figure shows that as the penalization parameter η 1 decreases, the solution approaches zero, indicating that the boundary condition is imposed more accurately.Note that in the last case, a large penalization parameter (i.e.weaker penalization) η 1 = 10 −3 is used.In addition, when the first-order term is added, improved accuracy is seen as the solution is closer to zero.The errors in the fluid region of the four cases are 3.071 • 10 −2 , 5.385 • 10 −3 , 5.698 • 10 −4 , and 1.022 • 10 −4 , respectively.This indicates that by introducing the first-order term with a proper selection of the penalization parameter, it is possible to improve the accuracy.
Furthermore, to study the effect of η 2 , we run additional simulations for a range of η 2 , and show the errors in Figure 5.The errors in the flow and solid regions for η 1 = 10 −3 , η 1 = 10 −4 , and η 1 = 10 −5 are shown in Figures 5a, 5b, 5c and 5d, respectively.For consistency with the analysis of modified equations, the first group of cases is performed in polynomial order N = 2, and the second group of cases is performed in polynomial order N = 3. Improved accuracy is seen when the penalization parameter is decreased.In addition, there exists an optimal η 2 that leads to minimal errors in both the flow and solid regions, which is the same for all penalization parameters.This optimal value is η 2 = −1/c, indicating that inside the solid the first-order penalization term becomes −∂u/∂x thus the physical advection is canceled out.From Figure 5b and Figure 5d, this cancelation will lead to almost zero error inside the solid, indicating that the boundary condition is satisfied exactly.At a larger η 1 , this optimal value remains valid, but the optimal error increases, as shown in Figure 5c and 5d.Therefore, to reach the optimal accuracy, we need to use a small penalization parameter η 1 , in combination with the optimal η 2 .These findings are consistent with the theory that the volume penalization modeling error converges with η 1 → 0. Furthermore, the conclusion of the modified equation analysis is also validated, since choosing η 2 = −1/c leads to improved accuracy and almost satisfies the boundary conditions exactly.
To investigate the effect of the viscous term, the advection-diffusion equation is investigated.Since the optimal η 2 in the advection equation has been obtained, η 2 = −1/c is selected for all cases.We proceed as for the advection equation, we set K = 40 elements, ∆ s = ∆x, r = 1/40 and N = 3.The initial condition with wavenumber ω∆x/(N + 1) = 0.3223 is used and marched in time to t = 1.5.Taking into account the effect of diffusion, the error in the flow region is limited to For the discretization of the viscous flux, either the BR1 or the LDG scheme is considered.The results for two physical viscosities ν = 0.001 and ν = 0.01 are shown in Figures 6 and Figures 7, respectively.For both cases, it is observed that the optimal second-order coefficient η 3 exists and can lead to a minimal error within the solid (as shown in Figure 6b and 6d and Figure 7b and 7d).This optimal value shows the relationship 1/η 3 = ν, which also indicates the cancelation of the viscous term inside the solid.This agrees with the optimal η 3 derived from the modified equation analysis.However, when looking at the error inside the flow, the optimal second-order penalization term does not lead to the lowest error when the BR1 scheme is used.This highlights the importance of choosing appropriate Riemann solvers to maintain good accuracy in the flow region.When the LDG scheme is selected, the optimal η 3 will reach the lowest error in the flow region, indicating that this flux is more suitable for the present problem.Therefore, when handling the viscous term, the LDG scheme is preferred, which gives consistent results of errors in the solid and in the fluid.In summary, the one-dimensional test case shows that the optimal penalization parameters derived from the modified equation analysis achieve minimal numerical errors in imposing the boundary conditions.

Two-dimensional advection-diffusion equation
In this section, a numerical experiment is performed for the two-dimensional advection-diffusion equation, using the conclusions of the modified equation analysis.The one-dimensional test case in the previous section is extended to two space directions.Therefore, the optimal parameters derived from onedimensional test cases are then dependent on each space direction.Extensions for GL points can be found in [19].The governing equation is (the solid region is the non-slip wall): where the advection flux is f adv = (c x u, c y u) T , the diffusion flux f diff = (−ν x ∂u/∂x, −ν y ∂u/∂y) T , g = (1/η 2,x , 1/η 2,y ) T , and The first-order and second-order penalization parameters in each direction is denoted by the second subscript.Here we set c x = c y = 1 and ν x = ν y = 0.001, therefore, the optimal parameters satisfy η 2,x = η 2,y and η 3,x = η 3,y .Note that, for the present linear equation, the extension to different advection velocities and viscosities in each direction is straightforward, while the optimal penalization parameter (i.e., trivial solution from modified equation analysis) also varies in different directions.As in the one-dimensional test case, the solid wall is considered in the middle of the computational domain.A schematic illustration of the two-dimensional problem for the present study is shown in Figure 8.The solid non-slip region has a L shape, which is centered in the middle of the square domain, making the top right region amplified by the wall.If the advection direction is set appropriately, the initial wave will move towards the wall.After that, we can solve the equation until all the solutions in the top right region have been penalized (which are then expected to be zero), and compute the error in this region.The error is again the difference between the numerical and exact solution (here set to zero).We consider a square computational domain in x ∈ [−0.1, 0.1] and y ∈ [−0.1, 0.1], with periodic boundary conditions.The domain is discretized into 20 equispaced elements in both the x and the y directions, resulting in 400 square elements in total and uniform mesh size ∆x = ∆y = 0.01.The penalization parameter and the explicit time step is set to η 1 = ∆t = 10 −4 .The polynomial order N = 3 is selected.Due to the preset flow advection parameters, the advection moves towards the top right direction.The width of the solid region is set to the size of a uniform grid ∆ s = ∆x, resulting in the solid ratio r = 1/20.We use the wavelike initial condition u(x, y) = sin(ωx + ωy), where a nondimensional wavelength ω∆x/(N + 1) = 0.3307 is selected.The first simulation for only advection is performed when only the first penalization term (η 1 for u) is included.Figure 9 shows three typical solution fields at different times.As shown in the figure, the penalized solution will move towards the top right corner, and finally dominate the entire domain due to the periodic boundary conditions.To compare the accuracy of simulation, the final simulation time is set to 0.11.Two solution fields, without and with the optimal first-order penalization term, are shown in Figure 10, where the values of η 2 , x and η 2 , y are set to −1 to match the physical advection speed.The improved accuracy from adding the optimal first-order term is seen in both the solution field and in the error.The error in the fluid region has been greatly reduced from 0.0207 to 1.4616 • 10 −5 .This numerical experiment extends and validates the conclusions obtained from the modified equation analysis, where the optimal first-order penalization term cancels the advection term and leads to improved accuracy.
Additional numerical experiments are performed for the advection-diffusion equation.The space and time discretizations remain the same as in the ad- Diffusive Flux Scheme vection case.The final time is set to t = 0.15.Three strategies are considered: 1) only volume penalization for the non-slip wall boundary condition, 2) volume penalization for both the value and the first-order term, and 3) volume penalization for all terms.Two types of viscous fluxes with either the BR1 or the LDG scheme are considered.The errors inside the fluid region are compared in Table 4, where conclusions similar to one-dimensional advection can be drawn.When the BR1 scheme is used, adding additional first-order and second-order penalization terms improves the overall accuracy, compared with the standard case (the first strategy).However, the addition of a second-order term does not lead to improved accuracy in the flow region.When the LDG scheme is used, adding first-and second-order terms will lead to a greater reduction of the error.This is consistent with the observations for the onedimensional advection equation, where the LDG scheme is shown to provide more accurate results than the BR1 scheme.This numerical experiment validates the proposed modified equation analysis for the second-order derivative in two-dimensional linear equations.

Conclusions
This study contributes to a better understanding of the numerical errors for Immersed Boundary Methods based on volume penalization, in combination with a high-order nodal discontinuous Galerkin scheme.For this purpose, an analysis of the modified equation is provided.The modified equation is a useful tool to analyze dissipative/dispersive errors related to the numerical discretizations.In this paper, we focus on the spatial errors introduced by the Immersed Boundary Method.Nodal solutions are expanded as Taylor series, and by rearranging the pseudo-differential equation new terms arise.These terms allow us to obtain insight into the dissipative/dispersive character of the errors and guidelines for their minimization.For example, the inclusion of extra penalization terms of the first and second derivatives, in addition to the classic penalization of the variable.Through this analysis, we provide optimal values for the first-and second-order penalization parameters to cancel the advection/diffusive errors inside the solid, which in turn lead to improved errors in the flow.
Numerical experiments validate the theoretical findings obtained from the analysis of modified equations, where optimal penalization parameters can lead to minimal errors (with a sufficiently small penalization parameter η 1 ).When combined with an appropriate numerical scheme (here, Local discontinuous Galerkin for viscous terms), minimal errors in the flow region are reached.
Future work will extend these findings to systems of partial differential equations with non-linearities, and extend the theoretical analysis to multidimensional systems.

A The DGSEM technique
We re-write Equation (3a) in its weak form: where ψ = ψ(x, t) is a local smooth test function.Given that Ω = [0, L] is divided into K elements, the integral is split into the sum of element integrals: Each element, x = x(ξ), is transformed according to: Then, dx = (∆x k /2)dξ and ∂/∂x = (2/∆x k )∂/∂ξ.Thus the weak form becomes: Assuming that global variables are represented by K local polynomial variables and substituting the Lagrange interpolation of the test function into the Galerkin weak form, ψ = ψ j l j , we get the following: for k = 1, 2, . . ., K and j = 0, 1, . . ., N .The first and third integrals are evaluated as follows: whereas the second integral is integrated by parts, being l j = dl j /dξ.The VP flux function in the first term is substituted by a numerical flux, i.e., depending on the normal at the boundary, ±e k x , and the solution at two adjacent elements.We discuss the choice of the numerical flux later.The remaining integral is divided as follows: Substituting the integrals ( 35), ( 36), (37), and (40), we get the following: for k = 1, 2, . . ., K and j = 0, 1, . . ., N where the inner product of the given functions a = a(ξ) and b = b(ξ) is defined as follows: Additionally, the VP diffusive flux involves the derivative of u and must be discretized consistently with the rest of the scheme.If we write Equation (??) in weak form, and repeat the interpolating and integration-by-part procedures, we get: for k = 1, 2, . . ., K and j = 0, 1, . . ., N where U k is an analogy of the numerical flux for the solution, that is, The computation of the inner products is done via Gaussian quadrature: where wm are the Gauss-Lobatto weights ( N m=0 wm = 2) and χ When all of them are combined, Equations ( 7) are obtained.
Finally, the last stage of a DGSEM is the calculation of F and U to reproduce the physics of advection and diffusion.A variety of fluxes are available for DG, and most are summarized by Arnold et al. [4].Here, we use a unifying function: for λ ∈ R. If λ = 0 the discretization becomes a central scheme; λ = −1, upwind; λ = 1, downwind.The subscript "+" means the right boundary element and "−" the left boundary element.We also denote {{•}} as the averaging operator: and [[•]] as the jump operator: Once these operators are defined, we divide the numerical flux into an advective term and a diffusive term: F = F adv + F diff .The computation of the advective numerical flux is as follows: and the diffusive numerical flux is: The values of f diff at the boundary elements are computed with: Finally, U is computed as BR1 is recovered by setting α = −1 and β = γ = δ = 0, while LDG is obtained by setting α = γ = −1 and β = −δ = −1.The weights f and g of ( 8) are obtained by finding the us at x k−1 and x k from the function W described in Equation ( 52).

B Non-trivial solutions
In all the cases, the considered element is inside the solid region.The first case is related to an inviscid problem without a second derivative penalty term or a viscid problem with η 3 = 1/ν.The second case includes second derivatives and is therefore more general.
Case 1: Problem with ν = 0 The parameters T E and HOT are listed in Tables 5 and 6.The main findings include the following: In total, we have five unknowns ( c, ) and, therefore, a determined system would be: ∼ O ∆x 3 k for j = 0, 1, 2 within Ωs.However, the unique solution of the system is the trivial one.If we leave η 2 free, the system that determines the numerical flux weights becomes being T E k j ∼ Ж (3)k j ∼ O ∆x 2 k for j = 0, 1, 2. A non-trivial solution would be: Alternatively, if the upwiding numerical flux is the solution, i.e.
Then the system becomes: whose truncation error leads to: Setting u k−1 h,2 = u k h,0 is very similar to using a Continuous Galerkin (CG) method.A downwind numerical flux mimics the results of upwinding, but for the right-hand boundary element.Other numerical fluxes leave: Table 5: The reaction parameter and the coefficient Ж in the modified equations for a three-point GL grid and a problem with ν = 0.
Table 6: Numerical source in the DG source, s k DG,j = S k j − r k j u k h,j , for a problem with ν = 0.

Case 2: Problem with ν = 0
In this second case (with second derivatives) we consider η 3 free.The parameters of T E and HOT are listed in Tables 7 and 8. Again, we conclude that and Since c k j = c for j = 0, 1, 2, the term c k j − ∆ξ c k j in the truncation error should be suppressed since it is O(∆x −1 k ).The choice is g k 2 = g k 0 = 2.However, ν k j − ∆ξ ν k j = 0 for j = 0, 1, 2 and therefore T E k j ∼ O ∆x 0 k .If we want to find an optimal value of η 3 to increase the order of the scheme, we come to the conclusion that ν = 0, but this case was already discussed previously.Keeping g k 2 = g k 0 = 2 and η 3 free, g k−1 2 = g k+1 0 = 0 kills s k DG,1 , to have s k DG,0 , s k DG,2 = 0, In this case, a solution of the system is as follows: Additionally, if upwind in such a way that the second equation of the system is met, but the first one becomes: To eliminate this term, a relation of ηs is obtained, c + (4/∆x k ) ν = 0, since in a DG method u k−1 h,2 = u k h,0 .Note that in a CG method, it is not necessary to fill this relation, since the solution is continuous between elements.In all cases described previously, T E k j ∼ O(∆x 0 k ) for j = 0, 1, 2.
A summary of all the conditions derived can be found in Table 3. Table 7: The reaction parameter and the coefficient Ж in the modified equations for a three-point GL grid and a problem with ν = 0.

C Connection between Fourier series and Taylor series
Table 8: Numerical source in the DG source, s k DG,j = S k j − r k j u k h,j , for a problem with ν = 0.
with ∆x = x j+1 − x j = x j − x j−1 and u j = u(x j , t).From the point of view of the Taylor series, the nodal solutions are expanded as and, therefore, the modified equation is From the point of view of the Fourier analysis, we start with a trial solution in the form: u(x j , t) = uω(t) exp(iωx j ), where i is the imaginary unit and ω the wavenumber.The semi-discrete equation ( 84 Therefore, substituting Equation (90) into Equation (89) we obtain the modified equation: which is the same as the modified equation ( 86) derived from the Taylor series.A similar and comprehensive understanding of the discrete errors by Fourier analysis can be found in [59].This connection is not coincidental, indeed, for the Taylor series we can obtain the trial solution and vice versa, e.g., Although a Taylor series approximates a function point-wise and a Fourier series approximates a function globally, both methods yield the same modified equation and later the same truncation error.This achievement can facilitate extension to multidimensional problems using the Fourier framework; see [60].

Fig. 1 :
Fig. 1: Domain decomposition and reference interval in the DGSEM technique.

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

Fig. 4 :
Fig. 4: Mean squared error between sharp and smooth mask function with increasing δ.

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

Assume that we analyze the errors of the advection equation,
non-zero constant, c > 0. For example and simplicity, the semidiscretization using a central finite difference is

Table 2 :
Numerical source for the DG source,

Table 4 :
Error comparison of the advection equation with IBM wall under different diffusive flux schemes and different combinations of penalization parameters.