A sharp interface immersed boundary method for solving flow with arbitrarily irregular and changing geometry
\normalsize In this paper, a sharp interface immersed boundary method is developed for efficiently and robustly solving flow with arbitrarily irregular and changing geometries. The proposed method employs a three-step prediction-correction flow reconstruction scheme for boundary treatment and enforces Dirichlet, Neumann, Robin, and Cauchy boundary conditions in a straightforward and consistent manner. Numerical experiments concerning flow of two and three space dimensions, stationary and moving objects, convex and concave geometries, no-slip and slip wall boundary conditions are conducted to demonstrate the proposed method.
arxiv_huangrui_mo.tex
\documentclass{article}
\pdfoutput=1
\usepackage[affil-it]{authblk}
\usepackage{graphicx}
\usepackage{appendix}
\usepackage{amsmath}% math support
\usepackage[square,comma,numbers,sort&compress]{natbib}% natbib support
\usepackage[usenames,dvipsnames]{xcolor}% color support
\usepackage{mathptmx}% times font
\usepackage{enumitem}% customizing enumerated lists
\usepackage{subcaption}% subfigures
\usepackage[section]{placeins}% to prevent floats from being moved over section

\graphicspath{{./}{./img/}}

\bibliographystyle{elsarticle-num-names}
%%%%%%%%%%%%%%%%%%%%%%%

\begin{document}

\title{A sharp interface immersed boundary method for solving flow with arbitrarily irregular and changing geometry}

\author[1]{Huangrui Mo%
\thanks{Email: \texttt{huangrui.mo@uwaterloo.ca}; Corresponding author}}

\author[1]{Fue-Sang Lien}

\author[2]{Fan Zhang}

\author[1]{Duane S. Cronin}
\affil[1]{Department of Mechanical Engineering, University of Waterloo, 200 University Avenue West, Waterloo, ON N2L 3G1, Canada}
\affil[2]{Defence R$\&$D Canada -- Suffield, P.O. Box 4000, Medicine Hat, AB T1A 8K6, Canada}

\maketitle

\begin{abstract}
\normalsize
In this paper, a sharp interface immersed boundary method is developed for efficiently and robustly solving flow with arbitrarily irregular and changing geometries. The proposed method employs a three-step prediction-correction flow reconstruction scheme for boundary treatment and enforces Dirichlet, Neumann, Robin, and Cauchy boundary conditions in a straightforward and consistent manner. Numerical experiments concerning flow of two and three space dimensions, stationary and moving objects, convex and concave geometries, no-slip and slip wall boundary conditions are conducted to demonstrate the proposed method. 
\end{abstract}

\section{Introduction}

Particle jets, a rapid and nonuniform dispersal of granular media under impulsive energy release, are observed in many physical processes such as explosions with explosive charges surrounded by or mixed with solid particles \cite{frost2012particle, zhang2014large}, explosive volcanic eruptions \cite{kedrinskiy2009hydrodynamic}, and impact of a solid projectile on granular media \cite{lohse2004impact, pacheco2011impact}.

In order to understand this jet phenomenon, many studies have been conducted over the past decades \cite{zhang2009shock}. However, the formation mechanism of particle jets is still unidentified \cite{rodriguez2013solid}, which is primarily due to the complexity of shock-particle interactions \cite{rodriguez2013solid} and random nature of the force-chain networks in granular materials \cite{liu1995force, jaeger1996granular}.

To develop a numerical solver for fully investigating the particle jet phenomenon, a computational fluid dynamics solver for simulating shock-particle interactions under strong shocks and complex particle configurations is a prerequisite, of which the capability of dealing efficiently and robustly with arbitrarily irregular and changing geometries is crucial.

In the immersed boundary method introduced by \citet{peskin1972flow}, numerically solving flow with irregular geometries is conducted on generic Cartesian grids, which can greatly simplify the grid generation process and take advantages of some important features of modern high-performance computing architecture \cite{jung2001two, peskin2002immersed, mittal2005immersed}. For instance, reduced memory requirements and the capability to use linear arrays with linear indexing techniques as main data structures can benefit from modern processor prefetcher and hierarchal-cache architecture to achieve a highly efficient numerical solver.

Extensions of the immersed boundary method have been continuously developed to increase interface resolution and relax stability constraints \citep{fadlun2000combined, tseng2003ghost, uhlmann2005immersed, mori2008implicit, kapahi2013three, yang2015non, schwarz2016immersed}. \citet{yusof1997combined} and \citet{fadlun2000combined} developed the direct forcing immersed boundary method. In the direct forcing approach, boundary forces are implicitly imposed via flow reconstruction, which simplifies the numerical discretization procedure considerably and is well-suited for problems with rigid boundaries. \citet{balaras2004modeling} later improved the reconstruction procedure of direct forcing and applied to large-eddy simulations. Integrating ideas from the ghost fluid method \cite{fedkiw1999non, fedkiw2002coupling} and the direct forcing immersed boundary treatment \cite{fadlun2000combined, iaccarino2003immersed}, \citet{tseng2003ghost} systematically developed a polynomial reconstruction based ghost-cell immersed boundary method to further increase implementation flexibility while maintaining sharp interfaces \cite{mittal2005immersed, tseng2003ghost}. \citet{kapahi2013three} proposed a least square interpolation approach and applied to solving high velocity impact problems. The development of immersed boundary methods was comprehensively reviewed by \citet{peskin2002immersed}, \citet{mittal2005immersed}, and \citet{sotiropoulos2014immersed}.

The robustness of a direct forcing immersed boundary method highly depends on the numerical stability and stencil adaption capability of the employed interpolation method \citep{tseng2003ghost, gao2007improved, kapahi2013three}. Polynomial reconstruction based methods frequently involve constructing linear systems on neighbouring stencils, including a nearby boundary point, of the interpolated node. When one of the stencils is very close to the boundary point, the resulting linear systems may suffer from numerical singularities \citep{tseng2003ghost, gao2007improved}. Additionally, a fixed minimum number of stencils is always required to avoid under-determined linear systems. Therefore, special treatments are required when strongly concave or convex geometries exist \citep{gao2007improved, kapahi2013three}. To enhance numerical stability and stencil adaption capability, the idea of using inverse distance weighting interpolation is firstly introduced by \citet{tseng2003ghost}; and a hybrid Taylor series expansion / inverse distance weighting approach was later developed by \citet{gao2007improved}.

In addition to numerical stability and stencil adaption capability, correctly enforcing different types of boundary conditions in a straightforward and consistent manner is another vital factor in obtaining an efficient, accurate, and robust immersed boundary method, since a variety of boundary conditions require to be repeatedly enforced on numerical boundaries with a large number of computational nodes \cite{gibou2013high}. In solving Navier-Stokes equations, constant temperature at a wall and velocity at a no-slip wall have Dirichlet boundary conditions; pressure at a wall and temperature at an adiabatic wall have Neumann boundary conditions; and velocity at a slip wall has a type of Cauchy boundary conditions.

Excluding the Dirichlet boundary conditions in which boundary values are determined and known, the enforcement of other types of boundary conditions, particularly Cauchy boundary conditions, for immersed boundaries demands considerable efforts \cite{crockett2011cartesian, kempe2015imposing, schwarz2016immersed}. Recently, \citet{kempe2015imposing} firstly devised a numerical implementation of slip wall boundary conditions in the context of immersed boundary methods. However, the realization is not straightforward due to its complexity \citep{kempe2015imposing}. Therefore, to enforce a variety of boundary conditions in a straightforward and consistent manner is a considerably challenging task.

To achieve an efficient and robust boundary treatment method for solving flow with arbitrarily irregular and changing geometries on Cartesian grids, this work develops a sharp interface immersed boundary method. By the development of an inverse distance weighting based three-step prediction-correction flow reconstruction scheme for boundary treatment, the proposed method enforces Dirichlet, Neumann, Robin, and Cauchy boundary conditions in a straightforward and consistent manner. The developed method serves the objective of solving flow interacting with multiple objects involving collision, agglomeration, penetration, and fragmentation processes. 

This paper is structured as below. Section~\ref{sec:framework} presents the sharp interface immersed boundary method under a generalized framework of ghost-cell immersed boundary treatment. Section~\ref{sec:numeric} describes the employed three-dimensional Navier-Stokes solver of this paper. Section~\ref{sec:validity} validates the proposed immersed boundary method. Section~\ref{sec:conclusion} draws conclusions. 

\section{A sharp interface immersed boundary method}\label{sec:framework}

\subsection{Generalized framework}

Fig.~\ref{fig:gcibm_demo} shows 2D and 3D schematic diagrams of a computational domain with an immersed boundary. $G$ denotes a ghost node, a computational node that locates at the numerical boundaries but outside the physical domain. $O$ denotes a boundary point with $\mathbf{GO}$ as the outward normal vector. $I$ is the image point of ghost node $G$ reflected by the physical boundary. 
\begin{figure}[!htbp]
    \centering
    \begin{subfigure}[b]{0.48\textwidth}
        \includegraphics[width=\textwidth]{gcibm_demo_2D}
        \caption{}
        \label{fig:gcibm_demo_2D}
    \end{subfigure}%
    ~
    \begin{subfigure}[b]{0.48\textwidth}
        \includegraphics[width=\textwidth]{gcibm_demo_3D}
        \caption{}
        \label{fig:gcibm_demo_3D}
    \end{subfigure}%
    \caption{Schematic diagrams of a computational domain with an immersed boundary. $G$, ghost node; $O$, boundary point; $I$, image point. (a) 2D space. (b) 3D Space.}
    \label{fig:gcibm_demo}
\end{figure}

In the ghost-cell immersed boundary method \cite{tseng2003ghost}, the numerical boundary treatment is a reconstruction process of variable values at numerical boundaries via physical boundary conditions and variable values at interior physical domain. To construct the flow at numerical boundaries while admitting the existence of physical boundaries, the method of images \cite{chu1974boundary} is an effective way \cite{colella1990multidimensional, tseng2003ghost}. Therefore, the reconstruction of a generic flow variable $\psi$ at a ghost node $G$ is a two-step approach:
\begin{subequations} \label{eq:gcibm}
    \begin{align}
        \psi_G  &= 2\psi_O - \psi_I \label{eq:gcibma} \\
        \psi_I  &= f(x_I, y_I, z_I) \label{eq:gcibmb}
    \end{align}
\end{subequations}
where $f(x, y, z)$ is a local reconstruction function of $\psi$ at spatial point $I$. 

Generally, the reconstruction function needs to be determined by physical boundary conditions and known values of $\psi$ at nearby fluid nodes. That is,
\begin{equation} \label{eq:gcibmb_another}
   \psi_I = f(\{\psi_N\}, \psi_O) 
\end{equation}
where, $\psi_O$ is the value of $\psi$ at the boundary point $O$ where physical boundary conditions are enforced; and $\{\psi_N\}$ represent values of $\psi$ at fluid nodes $\{N\}$ that satisfy:
\begin{equation} \label{eq:domain}
    d_N = ||\mathbf{r}_I - \mathbf{r}_N|| \le R_I
\end{equation}
where, $\mathbf{r}_I$ and $\mathbf{r}_N$ are the position vectors of point $I$ and $N$, respectively; $R_I$, referred as the domain of dependence of point $I$ and illustrated in Fig.~\ref{fig:gcibm_demo}, is the maximum distance from the point $I$ to nearby fluid nodes that are employed for flow reconstruction. 

The benefit of the incorporation of physical boundary conditions in the reconstruction function of $\psi_I$ can be demonstrated by the following equations:
\begin{subequations} \label{eq:limit}
    \begin{align}
        ||\mathbf{r}_G - \mathbf{r}_O|| &= ||\mathbf{r}_I - \mathbf{r}_O|| \\
        \lim_{||\mathbf{r}_I - \mathbf{r}_O|| \to 0}\psi_I  &= \lim_{||\mathbf{r}_I - \mathbf{r}_O|| \to 0}f(\{\psi_N\}, \psi_O) = \ \psi_O \\
        \lim_{||\mathbf{r}_G - \mathbf{r}_O|| \to 0}\psi_G  &= 2\psi_O - \lim_{||\mathbf{r}_I - \mathbf{r}_O|| \to 0}\psi_I = \psi_O
    \end{align}
\end{subequations}

Hence, the constructed $\psi_G$ converges to the exact physical boundary conditions when $G$ converges to $O$. This convergence property is helpful to alleviate unphysical flux over the immersed boundary, an issue resulting from using non-body conformal Cartesian grids \cite{mittal2005immersed, mark2008derivation, seo2011sharp} and examined in the numerical results of this paper.

Several approaches are available to construct flow at image points \cite{tseng2003ghost, gilmanov2003general, gilmanov2005hybrid, gao2007improved, kapahi2013three}. In this study, we develop an inverse distance weighting based flow reconstruction function to achieve efficient and robust boundary treatment for flow with arbitrarily irregular and changing geometries and to enforce Dirichlet, Neumann, Robin, and Cauchy boundary conditions in a straightforward and consistent manner.

\subsection{Inverse distance weighting interpolation}

As a convex combination of candidate stencils, inverse distance weighting is a popular interpolation method for the approximation of scattered data sets \cite{shepard1968two, junkins1973weighting}. The inverse distance weighting \cite{shepard1968two} for interpolating the value of a variable $\psi$ at a spatial point $c$ is expressed as the following:
\begin{equation}\label{eq:weighting}
    \psi_c = \frac{\sum w(d_n) \psi_n }{\sum w(d_n)}, \text{~~~~} d_n \ne 0 \ \text{and} \ \ d_n \le R_c
\end{equation}
where $\psi_c$ is the interpolated value; $\{\psi_n\}$ are the known values of $\psi$ at stencil points $\{n\}$; $\{d_n\}$ are the distance from $\{n\}$ to the interpolated point $c$; $\{w(d_n)\}$ are weighting functions; $R_c$ is the size of the domain of dependence for interpolating $\psi$ at point $c$.

As discussed in \cite{shepard1968two}, the desired $\lim_{d_n \to 0}\psi_c=\psi_n$ is mathematically satisfied. However, an overflow problem may arise from calculating an inverse distance. \citet{shepard1968two} suggested using a conditional statement to avoid this issue:
\begin{equation}
    d = max(d, d_{tiny})
\end{equation}
where $d_{tiny}$ is a predefined positive constant to avoid float arithmetic overflow of inversing distance. In this paper, this value is set as a function of mesh sizes:
\begin{equation}
    d_{tiny} = \epsilon_0 * min(\Delta_x, \Delta_y, \Delta_z)
\end{equation}
where, $\Delta_x, \Delta_y, \Delta_z$ are mesh sizes in $x$, $y$, $z$ directions respectively; $\epsilon_0$ is a constant representing the proportion of $d_{tiny}$ to the smallest mesh size, for instance, $\epsilon_0=1.0\times10^{-6}$.

Generally, the weighting function $w(d)$ employs an inverse-power law $1/d^q$, and the typical value of $q$ is $2$ \cite{shepard1968two, franke1982scattered}. Our numerical experiments on $q=1,\ 2$ with $R_c / \Delta_{max}=2,\ 4,\ 6$ ($\Delta_{max} = max(\Delta_x, \Delta_y, \Delta_z)$) indicate that computational results are not sensitive to the choice of $q$ and $R_c$. Hence, $q=2$ with $R_c / \Delta_{max} = 2$ are used for the numerical results of this paper.

\subsection{A three-step prediction-correction flow reconstruction scheme}

The proposed three-step prediction-correction scheme of this paper for constructing $\psi_I=f(\{\psi_N\}, \psi_O)$ with $\psi$ representing a generic field variable is presented as below:

\begin{enumerate}
    \item Prediction step: pre-estimate the value of $\psi_I$ by applying inverse distance weighting on the fluid nodes that locate in the domain of dependence of the image point $I$. Denote the predicted value as $\psi_I^*$.
\begin{equation}\label{eq:prediciton}
    \psi_I^* = \frac{\sum w(d_N) \psi_N}{\sum w(d_N)}
\end{equation}
    \item Physical boundary condition enforcement step: determine the value of $\psi_O$ via the physical boundary conditions that $\psi$ needs to satisfy at the boundary point $O$ and the values of $\psi$ at interior physical domain.
    \item Correction step: solve the value of $\psi_I$ by adding the boundary point $O$ as a stencil node for the inverse distance weighting of $\psi_I$.
\begin{equation}\label{eq:correction}
    \psi_I = \frac{\sum w(d_N) \psi_N + w(d_O) \psi_O}{\sum w(d_N) + w(d_O)} = \frac{\psi_I^* + \frac{w(d_O)}{\sum w(d_N)}\psi_O}{1+\frac{w(d_O)}{\sum w(d_N)}}
\end{equation}
It is beneficial to note that there is no need to re-do calculations on fluid nodes when the sum of weights and sum of weighted values in Eq.~\eqref{eq:prediciton} are preserved.
\end{enumerate}

The physical boundary condition enforcement step is described below by the implementation of practical boundary conditions.

\paragraph{Dirichlet boundary condition} If $\psi$ satisfies Dirichlet boundary condition, the value of $\psi_O$ is purely determined by the specified boundary condition:
\begin{equation} \label{eq:dirichlet}
    \psi_O = g
\end{equation}
where $g$ is a given value or function.

\paragraph{Neumann boundary condition} $\psi$ satisfies the following equation:
\begin{equation} \label{eq:neumann}
    \left. \frac{\partial \psi}{\partial n} \right|_O = \frac{\partial \psi_O}{\partial n}
\end{equation}
where ${\partial \psi_O}/{\partial n}$ is a given value or function. Rewritten Eq.~\eqref{eq:neumann} as the following:
\begin{equation} \label{eq:reneumann}
    \lim_{l \to 0} \frac{\psi(\mathbf{r}_O + l \mathbf{n}) - \psi(\mathbf{r}_O)}{l} = \frac{\partial \psi_O}{\partial n}
\end{equation}
where $\mathbf{r}_O$ is the position vector and $\mathbf{n}$ is the unit normal vector at boundary point $O$.

Since point $I$ is on the normal direction of point $O$, we have:
\begin{equation} \label{eq:direction}
    \mathbf{n} = \frac{\mathbf{r}_I - \mathbf{r}_O}{||\mathbf{r}_I - \mathbf{r}_O||}
\end{equation}

Therefore,
\begin{equation} \label{eq:approxneumann}
    \frac{\psi_I - \psi_O}{||\mathbf{r}_I - \mathbf{r}_O||} - \left. \frac{\partial^2 \psi}{\partial n^2} \right|_O ||\mathbf{r}_I - \mathbf{r}_O|| + \mathrm{O}(||\mathbf{r}_I - \mathbf{r}_O||^2)= \frac{\partial \psi_O}{\partial n}
\end{equation}

Due to Eq.~\eqref{eq:gcibma}, the second order derivative term is negligible:
\begin{equation} \label{eq:secondderivative}
    \left. \frac{\partial^2 \psi}{\partial n^2} \right|_O = \frac{\psi_I + \psi_G - 2\psi_O}{2||\mathbf{r}_I - \mathbf{r}_O||^2} + \mathrm{O}(||\mathbf{r}_I - \mathbf{r}_O||^2)
\end{equation}

Hence, $\psi_O$ is determined as:
\begin{equation} \label{eq:resultneumann}
    \psi_O = \psi_I - ||\mathbf{r}_I - \mathbf{r}_O||\frac{\partial \psi_O}{\partial n}
\end{equation}

\paragraph{Robin boundary condition} A linear combination of the values of $\psi$ and its normal derivative on the boundary point $O$ is specified:
\begin{equation} \label{eq:robin}
    \alpha \psi_O + \beta \left. \frac{\partial \psi}{\partial n} \right|_O = g
\end{equation}
where $\alpha$ and $\beta$ are the linear combination coefficients, $g$ is a given value or function.

After approximating the normal derivative, we have:
\begin{equation} \label{eq:approxrobin}
    \alpha \psi_O + \beta \frac{\psi_I - \psi_O}{||\mathbf{r}_I - \mathbf{r}_O||} = g
\end{equation}

Then,
\begin{equation} \label{eq:resultrobin}
    \psi_O = \frac{\beta \psi_I - ||\mathbf{r}_I - \mathbf{r}_O|| g}{\beta - ||\mathbf{r}_I - \mathbf{r}_O|| \alpha}
\end{equation}

\paragraph{Cauchy boundary condition}

For illustration purpose, $\psi$ is replaced by the velocity $\mathbf{V}=(u,v,w)$ that satisfies the slip wall boundary condition:
\begin{subequations} \label{eq:slipwall}
    \begin{align}
        \left. (\mathbf{V} \cdot \mathbf{n}) \right|_{\mathbf{r}=\mathbf{r}_O} &= \mathbf{V}_{S} \cdot \mathbf{n} \\
        \left. \frac{\partial (\mathbf{V} \cdot \hat{\mathbf{t}})}{\partial n} \right|_{\mathbf{r}=\mathbf{r}_O} &= 0 \\
        \left. \frac{\partial (\mathbf{V} \cdot \tilde{\mathbf{t}})}{\partial n} \right|_{\mathbf{r}=\mathbf{r}_O} &= 0
    \end{align}
\end{subequations}
where $\mathbf{n}$, $\hat{\mathbf{t}}$, and $\tilde{\mathbf{t}}$ are the unit normal vector, unit tangent vector, and unit bitangent vector at boundary point $O$, respectively. $\mathbf{V}_{S}$ is the velocity of the boundary surface.

After approximating normal derivatives, we have:
\begin{subequations} \label{eq:velocity}
    \begin{align}
        u_O n_x + v_O n_y + w_O n_z &= u_{S}  \hat{n}_x + v_{S}  \hat{n}_y + w_{S}  \hat{n}_z \\
        u_O \hat{t}_x + v_O \hat{t}_y + w_O \hat{t}_z &= u_I  \hat{t}_x + v_I  \hat{t}_y + w_I  \hat{t}_z \\
        u_O \tilde{t}_x + v_O \tilde{t}_y + w_O \tilde{t}_z &= u_I  \tilde{t}_x + v_I  \tilde{t}_y + w_I  \tilde{t}_z
    \end{align}
\end{subequations}

Since the coefficient matrix is orthogonal, $\mathbf{V}_O$ is determined as below: 
\begin{equation} \label{eq:determined}
    \begin{pmatrix}
        u_O \\
        v_O \\
        w_O
    \end{pmatrix}
    =
    \begin{bmatrix}
        n_x & n_y & n_z \\
        \hat{t}_x & \hat{t}_y & \hat{t}_z \\
        \tilde{t}_x & \tilde{t}_y & \tilde{t}_z
    \end{bmatrix}^{T}
    \begin{pmatrix}
        u_{S}  \hat{n}_x + v_{S}  \hat{n}_y + w_{S}  \hat{n}_z \\
        u_I  \hat{t}_x + v_I  \hat{t}_y + w_I  \hat{t}_z \\
        u_I  \tilde{t}_x + v_I  \tilde{t}_y + w_I  \tilde{t}_z
    \end{pmatrix}
\end{equation}

All the solution equations of $\psi_O$ now can be written in a unified form:
\begin{equation} \label{eq:unified}
    \psi_O = C \psi_I + R.R.H.S.
\end{equation}
where, the value of the coefficient $C$ and the rest right hand side $R.R.H.S.$ are in Table~\ref{tab:bcmap}.
\begin{table}[!htbp]
    \centering
    \caption{Value map of $C$ and $R.R.H.S.$ for different boundary conditions.}
    \label{tab:bcmap}
    \scriptsize
    \begin{tabular}{lccc}
        \hline
        Type & Example Form & $C$ & $R.R.H.S.$\\
        \hline
        Dirichlet & $\psi_O = g$ & $0$ & $g$ \\
        Neumann & $\left. \frac{\partial \psi}{\partial n} \right|_O = \frac{\partial \psi_O}{\partial n}$ & $1$ & $- ||\mathbf{r}_I - \mathbf{r}_O||\frac{\partial \psi_O}{\partial n}$ \\
        Robin & $\alpha \psi_O + \beta \left. \frac{\partial \psi}{\partial n} \right|_O = g$ & $\frac{\beta}{\beta - ||\mathbf{r}_I - \mathbf{r}_O|| \alpha}$ & $\frac{- ||\mathbf{r}_I - \mathbf{r}_O|| g}{\beta - ||\mathbf{r}_I - \mathbf{r}_O|| \alpha}$ \\
    Cauchy & $\begin{aligned} \left. (\mathbf{V} \cdot \mathbf{n}) \right|_{\mathbf{r}=\mathbf{r}_O} &= \mathbf{V}_{S} \cdot \mathbf{n} \\ \left. \frac{\partial (\mathbf{V} \cdot \hat{\mathbf{t}})}{\partial n} \right|_{\mathbf{r}=\mathbf{r}_O} &= 0 \\ \left. \frac{\partial (\mathbf{V} \cdot \tilde{\mathbf{t}})}{\partial n} \right|_{\mathbf{r}=\mathbf{r}_O} &= 0 \end{aligned}$ & $\begin{bmatrix} n_x & n_y & n_z \\ \hat{t}_x & \hat{t}_y & \hat{t}_z \\ \tilde{t}_x & \tilde{t}_y & \tilde{t}_z \end{bmatrix}^{T} \begin{bmatrix} 0 & 0 & 0 \\ \hat{t}_x & \hat{t}_y & \hat{t}_z \\ \tilde{t}_x & \tilde{t}_y & \tilde{t}_z \end{bmatrix}$ & $\begin{bmatrix} n_x & n_y & n_z \\ \hat{t}_x & \hat{t}_y & \hat{t}_z \\ \tilde{t}_x & \tilde{t}_y & \tilde{t}_z \end{bmatrix}^{T} \begin{bmatrix} n_x & n_y & n_z \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} \cdot \mathbf{V}_S$ \\
        \hline
    \end{tabular}
\end{table}

Due to the unknown $\psi_I$ in Eq.~\eqref{eq:unified}, the solution equation of $\psi_O$ is coupled with the solution equation of $\psi_I$ in the correction step. To solve this problem, one method is a synchronous solving approach to solve $\psi_O$ and $\psi_I$ simultaneously:
\begin{equation} \label{eq:synchronous}
    \begin{cases}
        \psi_O = \ C \psi_I + R.R.H.S. \\
        \psi_I = \ \frac{\psi_I^* + \frac{w(d_O)}{\sum w(d_N)}\psi_O}{1+\frac{w(d_O)}{\sum w(d_N)}}
    \end{cases}
\end{equation}

The other is an asynchronous solving approach: first, solve $\psi_O$ via approximating the unknown $\psi_I$ with the pre-estimated $\psi_I^*$; then, solve $\psi_I$ in the correction step.
\begin{equation} \label{eq:asynchronous}
    \begin{cases}
        \psi_O = C \psi_I^* + R.R.H.S. \\
        \psi_I = \frac{\psi_I^* + \frac{w(d_O)}{\sum w(d_N)}\psi_O}{1+\frac{w(d_O)}{\sum w(d_N)}}
    \end{cases}
\end{equation}

The enforcement of Dirichlet and trivial Neumann boundary conditions is equivalent in these two approaches. When the asynchronous solving approach is adopted, the physical boundary condition enforcement step and the correction step can be iteratively implemented to improve the accuracy of enforcing other boundary conditions.

\subsection{Method discussion}

The proposed three-step prediction-correction flow reconstruction scheme enables the developed immersed boundary method to enforce a wide class of boundary conditions in a straightforward and consistent manner. In polynomial reconstruction based methods, different linear systems require to be constructed and solved for flow variables satisfying different boundary conditions. In the current method, the enforcement of boundary conditions at boundary points is a separate step of flow reconstruction. This feature provides efficient and uniform boundary treatment for flow with an arbitrary number of field variables that satisfy different types of boundary conditions.

Moreover, the proposed immersed boundary method here is scalable to the number of stencils used in the flow reconstruction. In contrast, for polynomial reconstruction based methods, a fixed minimum number of stencils is always required to avoid under-determined linear systems, and special treatments are required for strongly concave or convex geometries \citep{gao2007improved, kapahi2013three}. Therefore, the scalable property is where the robustness of the method herein lies: it leads to the automatic adaption to a varying number of stencil nodes; and it guarantees uniform validity when at least one fluid node exists in the domain of dependence of the image point, a condition that is ensured by the definition of a ghost node. In addition, the proposed immersed boundary method can be applied to multiple layers of ghost nodes without extra constraints.

The asynchronous solving approach without iterative implementation is currently adopted and examined in this paper, since the validity of the synchronous solving approach is established when the validity of the asynchronous solving approach is proved. A code that features the implementation of the three-step prediction-correction flow reconstruction scheme is provided in the \ref{app:source_code} to illustrate the simplicity and efficiency of our proposed immersed boundary method.

\section{Numerical implementation} \label{sec:numeric}

The governing equations employed in the numerical solver of this paper is the nondimensionalized conservative form of the three-dimensional Navier-Stokes equations in Cartesian coordinates \cite{anderson1995computational}. The temporal derivatives in the governing equations are discretized by using the third-order TVD Runge-Kutta method \cite{shu1988efficient, gottlieb2001strong}. The second order upwind TVD scheme \cite{harten1984class} and the fifth order WENO scheme \cite{jiang1996efficient} are both implemented for the discretization of convective fluxes. Meanwhile, central differencing scheme is used for the discretization of diffusive fluxes \cite{ferziger2002computational}. Strang splitting \cite{strang1964accurate, strang1968construction} is employed for dimensional splitting to relax stability constraints and decrease numerical complexities \cite{woodward1984numerical}. The fluid-solid coupling pattern used in this paper is described in \ref{app:fluid_solid_coupling}.

\section{Numerical experiments} \label{sec:validity}

Due to the using of the Navier-Stokes equations, the no-slip wall condition with an adiabatic assumption is enforced in the test cases unless otherwise stated. Numerical results are computed by the fifth order WENO scheme unless otherwise stated.

\subsection{Shock diffraction over a cylinder} \label{case:1_cyn}

A Mach $2.81$ planar incident shock interacting with a stationary circular cylinder is considered. Comprehensive descriptions of this problem are available in \cite{bryson1961diffraction, bazhenova1984unsteady, kaca1988interferometric,ripley2006numerical, sambasivan2009ghostb}. In the current study, a circular cylinder with diameter $D = 1$ is positioned at the center of a $6D \times 6D$ square domain while an initial shock is positioned $0.5D$ upstream of the cylinder. This configuration is similar to \cite{ripley2006numerical} except that a full domain size without symmetrical boundary assumption is used in this paper.

\subsubsection{Grid sensitivity test}

A consistent numerical method should provide numerical solutions that become less sensitive to the grid size as the mesh is refined. However, the level of numerical error depends on the features of the flow resolved by the grid \cite{oberkampf2002verification}. The presence of discontinuities, such as shocks, slip surfaces, and interfaces, develops numerical errors on a grid in a complex way \cite{oberkampf2002verification} and affects the error estimation process, especially for the evaluation of local grid convergence behavior \cite{roache1998verification}. While admitting these difficulties, this part provides an examination of the grid convergence behavior of the presented immersed boundary method for solving flow with strong discontinuities.
\begin{figure}[!htbp]
    \centering
    \begin{subfigure}[b]{0.48\textwidth}
        \includegraphics[width=\textwidth]{grid_global_distribution}
        \caption{}
        \label{fig:grid_global_distribution}
    \end{subfigure}%
    ~
    \begin{subfigure}[b]{0.48\textwidth}
        \includegraphics[width=\textwidth]{grid_global_position}
        \caption{}
        \label{fig:grid_global_position}
    \end{subfigure}%
    \\
    \begin{subfigure}[b]{0.48\textwidth}
        \includegraphics[width=\textwidth]{grid_local_value}
        \caption{}
        \label{fig:grid_local_value}
    \end{subfigure}%
    ~
    \begin{subfigure}[b]{0.48\textwidth}
        \includegraphics[width=\textwidth]{grid_local_error}
        \caption{}
        \label{fig:grid_local_error}
    \end{subfigure}%
    \caption{Grid sensitivity test. (a) Global grid convergence at a $33^\circ$ line. (b) Global grid convergence on shock position. (c) Local grid convergence at a point with $Arc\ Length = 2.96$. (d) Local relative error at the sample point.}
    \label{fig:grid_sensitivity}
\end{figure}

A series of successively refined grids are employed to study the grid sensitivity of the developed immersed boundary method, and the numerical results are obtained at $t=1.0$. Due to the presence of complex discontinuity patterns, the global grid convergence behavior is represented by and examined on a line segment from point $(-0.27232, 0.41934)$ to point $(2.5, 2.21970)$, which is located on a $33^\circ$ tangent line of the cylinder and is plotted in Fig.~\ref{fig:triplepoint}.

As shown in Fig.~\ref{fig:grid_global_distribution}, when grid resolution changes significantly from the coarsest $251\times251$ grid to the finest $2001\times2001$ grid, excellent overall agreements are achieved on global grid convergence with main discrepancies occurring near the region of boundary interfaces and flow discontinuities. However, these discrepancies are effectively reduced when the grid resolution is sufficiently increased. By observing that the discrepancy of the predicted discontinuity is mainly about its position rather than magnitude, a global grid convergence behavior on discontinuity position is examined and shown in Fig.~\ref{fig:grid_global_position}, in which a well performed grid convergence is presented.

Since local errors are transported throughout the computational region and are strongly affected by flow discontinuities, an examination of the local grid convergence behavior has limited implications \cite{roache1998verification}. For instance, the predicted position of discontinuities, which highly depends on grid resolution, has major influence on the values of flow quantities near discontinuities. Nonetheless, the local grid convergence at a sample point, which corresponds to the peak value of the $2001\times2001$ grid solution, is examined in this paper. Fig.~\ref{fig:grid_local_value} shows the predicted density at this sample point, and Fig.~\ref{fig:grid_local_error} shows the numerical error relative to the predicted value of the finest grid. As shown in Fig.~\ref{fig:grid_local_error}, the local relative error of the predicted density at this sample point is effectively reduced as the grid is refined.

According to the discussed results, the developed immersed boundary method of this paper has well-behaved global and local grid convergence properties over a wide range of grid resolution.

\subsubsection{Numerical results}

Numerical results of three grids with $501 \times 501$, $1001 \times 1001$, and $2001 \times 2001$ nodes, denoted as grid $C$, grid $M$, and grid $F$ respectively, are discussed below.
\begin{figure}[!htbp]
    \centering
    \begin{subfigure}[b]{0.32\textwidth}
        \includegraphics[trim = 60mm 0mm 60mm 0mm, clip, width=\textwidth]{1_cyn_nomv_28}
        \caption{t=0.47445}
        \label{fig:1_cyn_nomv_28}
    \end{subfigure}%
    ~
    \begin{subfigure}[b]{0.32\textwidth}
        \includegraphics[trim = 60mm 0mm 60mm 0mm, clip, width=\textwidth]{1_cyn_nomv_36}
        \caption{t=0.60999}
        \label{fig:1_cyn_nomv_36}
    \end{subfigure}%
    ~
    \begin{subfigure}[b]{0.32\textwidth}
        \includegraphics[trim = 60mm 0mm 60mm 0mm, clip, width=\textwidth]{1_cyn_nomv_40}
        \caption{t=0.67719}
        \label{fig:1_cyn_nomv_40}
    \end{subfigure}%
    \\
    \begin{subfigure}[b]{0.32\textwidth}
        \includegraphics[trim = 60mm 0mm 60mm 0mm, clip, width=\textwidth]{1_cyn_nomv_48}
        \caption{t=0.81157}
        \label{fig:1_cyn_nomv_48}
    \end{subfigure}%
    ~
    \begin{subfigure}[b]{0.32\textwidth}
        \includegraphics[trim = 60mm 0mm 60mm 0mm, clip, width=\textwidth]{1_cyn_nomv_52}
        \caption{t=0.87914}
        \label{fig:1_cyn_nomv_52}
    \end{subfigure}%
    ~
    \begin{subfigure}[b]{0.32\textwidth}
        \includegraphics[trim = 60mm 0mm 60mm 0mm, clip, width=\textwidth]{1_cyn_nomv_60}
        \caption{t=1.00000}
        \label{fig:1_cyn_nomv_60}
    \end{subfigure}%
    \caption{Density contour of shock diffraction over a cylinder. Colored by velocity.}
    \label{fig:1_cyn_nomv}
\end{figure}

Numerical results of grid $M$ are illustrated in Fig.~\ref{fig:1_cyn_nomv} in the form of the time evolution of density contour lines. These snapshots clearly show the diffraction of the curved Mach stems over the cylinder and the formation of wake by the collision of the two opposite diffracting shocks, which physical processes are reported in experimental studies \cite{bryson1961diffraction, bazhenova1984unsteady}. Comparing with the interferometric measurements of \cite{kaca1988interferometric} and numerical results in \cite{sambasivan2009ghostb, ripley2006numerical, ji2010numerical}, the slip line, the reflected and diffracted shocks over the immersed boundary in the results of this paper are all resolved remarkably well, which demonstrates the high validity of the developed method.

To investigate the capability of the developed immersed boundary method on enforcing different physical boundary conditions, two types of wall boundary conditions are studied on grid $M$ and results at $t=1.0$ are shown in Fig.~\ref{fig:bc_compare}. After the formation of collision wakes as in Fig.~\ref{fig:1_cyn_nomv_60_noslip} and Fig.~\ref{fig:1_cyn_nomv_60_slip}, a more wall-adhesive low-speed wake is formed in the case of no-slip wall. Fig.~\ref{fig:1_cyn_nomv_60_vel_gradient} and Fig.~\ref{fig:1_cyn_nomv_60_vel_gradient_slip} show the velocity gradient distribution at the wall region. A thin but gradually growing boundary layer with large velocity gradients is produced along the no-slip wall, while no boundary layer is presented in the slip wall case.
\begin{figure}[!htbp]
    \centering
    \begin{subfigure}[b]{0.48\textwidth}
        \includegraphics[width=\textwidth]{1_cyn_nomv_60}
        \caption{}
        \label{fig:1_cyn_nomv_60_noslip}
    \end{subfigure}%
    ~
    \begin{subfigure}[b]{0.48\textwidth}
        \includegraphics[width=\textwidth]{1_cyn_nomv_60_vel_gradient}
        \caption{}
        \label{fig:1_cyn_nomv_60_vel_gradient}
    \end{subfigure}%
    \\
    \begin{subfigure}[b]{0.48\textwidth}
        \includegraphics[width=\textwidth]{1_cyn_nomv_60_slip}
        \caption{}
        \label{fig:1_cyn_nomv_60_slip}
    \end{subfigure}%
    ~
    \begin{subfigure}[b]{0.48\textwidth}
        \includegraphics[width=\textwidth]{1_cyn_nomv_60_vel_gradient_slip}
        \caption{}
        \label{fig:1_cyn_nomv_60_vel_gradient_slip}
    \end{subfigure}%
    \caption{Numerical results with two types of wall boundary conditions. (a) Density contour colored by velocity, no-slip wall. (b) Velocity gradient color map, no-slip wall. (c) Density contour colored by velocity, slip wall. (d) Velocity gradient color map, slip wall.}
    \label{fig:bc_compare}
\end{figure}
\begin{figure}[!htbp]
    \centering
    \includegraphics[width=0.48\textwidth]{1_cyn_nomv_vel_profile}
    \caption{X-velocity profiles at a line segment from point $(0,\ 0.50)$ to point $(0,\ 0.55)$.}
    \label{fig:1_cyn_nomv_vel_profile}
\end{figure}

Fig.~\ref{fig:1_cyn_nomv_vel_profile} further shows the x-velocity profiles at a vertical line segment from point $(0,\ 0.50)$ to point $(0,\ 0.55)$. A velocity profile that indicates a zero velocity at wall changing to local free-stream value away from the wall is presented in the no-slip wall case. A velocity profile that indicates a maximum velocity at wall decreasing to local free-stream velocity is presented in the slip wall case. These observations agree with the flow physics at this vertical line of flow over a cylinder.

According to the successful solutions of shock diffraction over cylinder with no-slip and slip wall boundary conditions, the developed immersed boundary method is able to correctly enforce different types of physical boundary conditions.
\begin{figure}[!htbp]
    \centering
    \begin{subfigure}[b]{0.32\textwidth}
        \includegraphics[trim = 60mm 0mm 60mm 0mm, clip, width=\textwidth]{1_cyn_nomv_triple_gridc_tvd}
        \caption{}
        \label{fig:1_cyn_nomv_triple_gridc_tvd}
    \end{subfigure}%
    ~
    \begin{subfigure}[b]{0.32\textwidth}
        \includegraphics[trim = 60mm 0mm 60mm 0mm, clip, width=\textwidth]{1_cyn_nomv_triple_gridm_tvd}
        \caption{}
        \label{fig:1_cyn_nomv_triple_gridm_tvd}
    \end{subfigure}%
    ~
    \begin{subfigure}[b]{0.32\textwidth}
        \includegraphics[trim = 60mm 0mm 60mm 0mm, clip, width=\textwidth]{1_cyn_nomv_triple_gridf_tvd}
        \caption{}
        \label{fig:1_cyn_nomv_triple_gridf_tvd}
    \end{subfigure}%
    \\
    \begin{subfigure}[b]{0.32\textwidth}
        \includegraphics[trim = 60mm 0mm 60mm 0mm, clip, width=\textwidth]{1_cyn_nomv_triple_gridc_weno}
        \caption{}
        \label{fig:1_cyn_nomv_triple_gridc_weno}
    \end{subfigure}%
    ~
    \begin{subfigure}[b]{0.32\textwidth}
        \includegraphics[trim = 60mm 0mm 60mm 0mm, clip, width=\textwidth]{1_cyn_nomv_triple_gridm_weno}
        \caption{}
        \label{fig:1_cyn_nomv_triple_gridm_weno}
    \end{subfigure}%
    ~
    \begin{subfigure}[b]{0.32\textwidth}
        \includegraphics[trim = 60mm 0mm 60mm 0mm, clip, width=\textwidth]{1_cyn_nomv_triple_gridf_weno}
        \caption{}
        \label{fig:1_cyn_nomv_triple_gridf_weno}
    \end{subfigure}%
    \\
    \begin{subfigure}[b]{0.32\textwidth}
        \includegraphics[trim = 60mm 0mm 60mm 0mm, clip, width=\textwidth]{1_cyn_nomv_triple_gridm_tvd_slip}
        \caption{}
        \label{fig:1_cyn_nomv_triple_gridm_tvd_slip}
    \end{subfigure}%
    ~
    \begin{subfigure}[b]{0.32\textwidth}
        \includegraphics[trim = 60mm 0mm 60mm 0mm, clip, width=\textwidth]{1_cyn_nomv_triple_gridm_weno_slip}
        \caption{}
        \label{fig:1_cyn_nomv_triple_gridm_weno_slip}
    \end{subfigure}%
    \caption{Superimposition of density contours (white lines are $33^{\circ}$ tangent lines of the cylinder). (a) TVD, grid $C$. (b) TVD, grid $M$. (c) TVD, grid $F$. (d) WENO, grid $C$. (e) WENO, grid $M$. (f) WENO, grid $F$. (g) TVD, grid $M$, slip wall. (h) WENO, grid $M$, slip wall.}
    \label{fig:triplepoint}
\end{figure}

\subsubsection{Shock detachment distance}

As shown in Fig.~\ref{fig:triplepoint}, a bow-shaped shock is reflected from the cylinder when the incident shock hits the cylinder. In \cite{kaca1988interferometric}, the concepts of shock detachment distance and nondimensional time after collision were introduced to describe the reflected incident shock. Fig.~\ref{fig:1_cyn_detach} shows the comparison of predicted detachment distance between present numerical results and the experimental results of \cite{kaca1988interferometric}.
\begin{figure}[!htbp]
    \centering
    \begin{subfigure}[b]{0.48\textwidth}
        \includegraphics[width=\textwidth]{1_cyn_detach_tvd}
        \caption{}
        \label{fig:1_cyn_detach_tvd}
    \end{subfigure}%
    ~
    \begin{subfigure}[b]{0.48\textwidth}
        \includegraphics[width=\textwidth]{1_cyn_detach_weno}
        \caption{}
        \label{fig:1_cyn_detach_weno}
    \end{subfigure}%
    \caption{Comparison of detachment distance. (a) TVD. (b) WENO.}
    \label{fig:1_cyn_detach}
\end{figure}
\begin{figure}[!htbp]
    \centering
    \begin{subfigure}[b]{0.48\textwidth}
        \includegraphics[width=\textwidth]{1_cyn_triple_path_tvd}
        \caption{}
        \label{fig:1_cyn_triple_path_tvd}
    \end{subfigure}%
    ~
    \begin{subfigure}[b]{0.48\textwidth}
        \includegraphics[width=\textwidth]{1_cyn_triple_path_weno}
        \caption{}
        \label{fig:1_cyn_triple_path_weno}
    \end{subfigure}%
    \caption{Comparison of triple-point path. (a) TVD. (b) WENO.}
    \label{fig:1_cyn_triple_path}
\end{figure}

As it is pointed out in \cite{ripley2006numerical} that numerical studies of the detachment distance generally predicts a greater detachment distance, and the detachment distance subjects to a parabolic distribution with respect to the nondimensional time after collision rather than the linear behavior reported in \cite{kaca1988interferometric}. This parabolic behavior of the detachment distance, which is clearly observed in the results of this paper, is also presented in polynomial reconstruction based results of \cite{sambasivan2009ghostb}, cut-cell based results of \cite{ji2010numerical}, and unstructured mesh based results of \cite{ripley2006numerical}.

\subsubsection{Triple-point path}

The incident shock, reflected shock, and diffracted shock intersect and form an upper triple point on each side of the plane of symmetry. As shown in Fig.~\ref{fig:triplepoint}, this triple point travels in space and produces an upper triple-point path. Interferometric measurements of \cite{kaca1988interferometric} predict that this upper triple-point path is tangent to the cylinder at an angle of $33^{\circ}$ for Mach numbers in the range of $1.42 - 5.96$.

The predicted triple-point paths of this paper and the experimental correlation of \cite{kaca1988interferometric} are compared in Fig.~\ref{fig:1_cyn_triple_path}. The least square linear regressions of the predicted triple-point paths of TVD with grid $C$, $M$, and $F$ are about $29.0^{\circ}$, $29.5^{\circ}$, and $29.6^{\circ}$ respectively, of WENO with grid $C$, $M$, and $F$ are about $30.2^{\circ}$, $30.3^{\circ}$, and $30.3^{\circ}$ respectively. For grid $M$ with slip wall boundary condition, TVD and WENO predict $30.1^{\circ}$ and $30.3^{\circ}$ respectively. These results, which agree well with the experimental correlation of \cite{kaca1988interferometric} and very well with polynomial reconstruction based results of \cite{sambasivan2009ghostb}, cut-cell based results of \cite{ji2010numerical}, and unstructured mesh based results of \cite{ripley2006numerical}, demonstrate the accuracy of the developed method.

\subsection{Explosive dispersal of zero-gap particles}

To test the robustness of the proposed method, strongly irregular, concave, and changing geometries formed by initially zero-gap configured particles are used. The zero-gap particles may represent one of the most challenging geometries in fluid-solid interactions. As preliminary cases for studying particle jetting instabilities from multiphase explosive detonation \cite{ripley2014jetting, zhang2014large}, no study case concerning zero-gap particles, which is the common setting when particles are packed in practice, has been found so far. 

As shown in Fig.~\ref{fig:18_cyn_domain}, eighteen identical particles are zero-gap configured in a $1 \times 1$ computational domain. The centers of particles are evenly distributed on a circle whose radius is equal to $0.2$. A flow state $(\rho, u, v, p)^T=(3.67372, 0, 0, 9.04545)^T$ is initially positioned at a circular region centered in the domain, and the radius of the circular region is $0.1$. The flow state at the rest of the region is set to $(\rho, u, v, p)^T=(1, 0, 0, 1)^T$. A grid with $2001 \times 2001$ nodes is used for the numerical solution. Newton's second law of motion is employed for evolving spatial distribution of particles, during which no collision and gravity effects are currently integrated.
\begin{figure}[!htbp]
    \centering
    \includegraphics[width=0.6\textwidth]{18_cyn_domain}
    \caption{Illustration of computational configuration of zero-gap particles.}
    \label{fig:18_cyn_domain}
\end{figure}

Numerical solutions of the explosive dispersal process are illustrated in Fig.~\ref{fig:18_cyn_dispersal}. Complex shock-shock interaction, shock reflection and diffraction behaviors are clearly resolved in the numerical results. In addition, compression waves are formed in front of the moving particles, and high velocity fluid jets are generated at regions between moving particles. These physical and successful solutions of explosive dispersal of zero-gap particles prove the high robustness of the proposed immersed boundary method for strongly irregular, concave, and changing geometries.
\begin{figure}[!htbp]
    \centering
    \begin{subfigure}[b]{0.30\textwidth}
        \includegraphics[trim = 60mm 0mm 60mm 0mm, clip, width=\textwidth]{18_cyn_mv_25}
        \caption{t=0.05124}
        \label{fig:18_cyn_mv_25}
    \end{subfigure}%
    ~
    \begin{subfigure}[b]{0.30\textwidth}
        \includegraphics[trim = 60mm 0mm 60mm 0mm, clip, width=\textwidth]{18_cyn_mv_28}
        \caption{t=0.05752}
        \label{fig:18_cyn_mv_28}
    \end{subfigure}%
    ~
    \begin{subfigure}[b]{0.30\textwidth}
        \includegraphics[trim = 60mm 0mm 60mm 0mm, clip, width=\textwidth]{18_cyn_mv_33}
        \caption{t=0.06789}
        \label{fig:18_cyn_mv_33}
    \end{subfigure}%
    \\
    \begin{subfigure}[b]{0.30\textwidth}
        \includegraphics[trim = 60mm 0mm 60mm 0mm, clip, width=\textwidth]{18_cyn_mv_44}
        \caption{t=0.09031}
        \label{fig:18_cyn_mv_44}
    \end{subfigure}%
    ~
    \begin{subfigure}[b]{0.30\textwidth}
        \includegraphics[trim = 60mm 0mm 60mm 0mm, clip, width=\textwidth]{18_cyn_mv_70}
        \caption{t=0.14389}
        \label{fig:18_cyn_mv_70}
    \end{subfigure}%
    ~
    \begin{subfigure}[b]{0.30\textwidth}
        \includegraphics[trim = 60mm 0mm 60mm 0mm, clip, width=\textwidth]{18_cyn_mv_80}
        \caption{t=0.16431}
        \label{fig:18_cyn_mv_80}
    \end{subfigure}%
    \caption{Explosive dispersal of zero-gap particles. Colored by pressure.}
    \label{fig:18_cyn_dispersal}
\end{figure}

\subsection{Shock diffraction over two partially overlapped spheres}

A Mach $2.81$ planar shock interacting with two stationary and partially overlapped spheres (Fig.~\ref{fig:2_sphere_domain}) is solved to test the proposed immersed boundary method for three-dimensional irregular and concave geometries. The diameter of these two identical spheres is $D = 1$, and the centers of them are located at $(0, 0, 0)$ and $(0.5D, 0, 0)$ respectively. The size of the computational domain is $6D \times 6D \times 6D$, and the computational grid is $251 \times 251 \times 251$ nodes.

Fig.~\ref{fig:2_sphere_rho} shows density contour on two perpendicular semi-planar slices. As shown in the two identical density contour slices, flow discontinuities resulted from the concave region are adequately resolved on this relatively coarse grid. Moreover, the axis symmetry property of this flow problem is able to be preserved.
\begin{figure}[!htbp]
    \centering
    \begin{subfigure}[b]{0.48\textwidth}
        \includegraphics[width=\textwidth]{2_sphere_domain}
        \caption{}
        \label{fig:2_sphere_domain}
    \end{subfigure}%
    ~
    \begin{subfigure}[b]{0.48\textwidth}
        \includegraphics[width=\textwidth]{2_sphere_slices}
        \caption{}
        \label{fig:2_sphere_rho}
    \end{subfigure}%
    \caption{Shock interacting with two partially overlapped spheres. (a) Computational domain. (b) Density contour on two perpendicular semi-planar slices. Colored by velocity.}
    \label{fig:2_sphere}
\end{figure}

Due to using non-body conformal Cartesian grids, unphysical flux over the immersed boundary is a fundamental issue in immersed boundary methods \cite{mittal2005immersed}, and considerable efforts to overcome this issue have been devoted in \cite{mark2008derivation, seo2011sharp}. Since Eq.~\eqref{eq:limit} has established the convergence property of current method that the constructed ghost flow is able to converge to exact physical boundary conditions, we now numerically examine the unphysical flux under practical grid sizes by synthesizing the solved problems.

Fig.~\ref{fig:stream_trace} presents the stream traces of shock diffraction problems. The solved stream traces by the developed immersed boundary method are closely aligned with geometry surfaces, even in the three-dimensional problem where a coarse grid is employed.
\begin{figure}[!htbp]
    \centering
    \begin{subfigure}[b]{0.48\textwidth}
        \includegraphics[width=\textwidth]{1_cyn_nomv_streamtrace}
        \caption{}
        \label{fig:1_cyn_nomv_streamtrace}
    \end{subfigure}%
    ~
    \begin{subfigure}[b]{0.48\textwidth}
        \includegraphics[width=\textwidth]{1_cyn_nomv_slip_streamtrace}
        \caption{}
        \label{fig:1_cyn_nomv_slip_streamtrace}
    \end{subfigure}%
    \\
    \begin{subfigure}[b]{0.48\textwidth}
        \includegraphics[width=\textwidth]{2_sphere_streamtrace}
        \caption{}
        \label{fig:2_sphere_streamtrace}
    \end{subfigure}%
    \caption{Stream traces colored by velocity with corresponding analytical geometry boundaries. (a) Shock diffraction over a cylinder, no-slip wall. (b) Shock diffraction over a cylinder, slip wall. (c) Shock diffraction over two partially overlapped spheres, no-slip wall.}
    \label{fig:stream_trace}
\end{figure}

To quantitatively examine the unphysical flux over immersed boundary, the absolute flux of the shock diffraction over a cylinder problems at time $t=1.0$ is examined and the flux is calculated on the first layer of ghost nodes.
\begin{equation} \label{eq:absolute_flux}
    \text{absolute flux} = \frac{1}{S}\iint_{S} \, |(\mathbf{V} - \mathbf{V}_{S})\ \cdot \ \mathbf{n}| \, \mathrm{d}s
\end{equation}

The absolute flux distribution $|(\mathbf{V} - \mathbf{V}_{S})\ \cdot \ \mathbf{n}|$ over the cylinder with a no-slip wall on grid $M$ is plotted in Fig.~\ref{fig:1_cyn_nomv_flux_distribution}, in which a symmetrical distribution is presented. The main unphysical flux is observed at the angular range of $\pm[60^{\circ},\, 100^{\circ}]$, where flow has a large velocity gradient at near wall region.
\begin{figure}[!htbp]
    \centering
    \begin{subfigure}[b]{0.48\textwidth}
        \includegraphics[width=\textwidth]{1_cyn_nomv_flux_distribution}
        \caption{}
        \label{fig:1_cyn_nomv_flux_distribution}
    \end{subfigure}%
    ~
    \begin{subfigure}[b]{0.48\textwidth}
        \includegraphics[width=\textwidth]{1_cyn_nomv_flux_overall}
        \caption{}
        \label{fig:1_cyn_nomv_flux_overall}
    \end{subfigure}%
    \caption{Absolute flux over analytical geometry boundaries of the shock diffraction over a cylinder problems. (a) Absolute flux distribution over the cylinder with no-slip wall on grid $M$. (b) Absolute flux over cylinder with no-slip or slip wall.}
    \label{fig:absolute_flux}
\end{figure}

Fig.~\ref{fig:1_cyn_nomv_flux_overall} plots the absolute flux of different cases, of which the average velocity and maximum velocity of flow fields are provided in Table~\ref{tab:velocity}. The slip wall and no-slip wall boundary conditions show similar flux values. When the grid is refined from $251\times251$ to $2001\times2001$, the absolute flux of the no-slip wall case is effectively reduced from about $1.44\times10^{-2}$ to about $4.19\times10^{-3}$, in which the latter is about $29.09\%$ of the former and is about $0.21\%$ of the average velocity of the flow field.
\begin{table}[!htbp]
    \centering
    \caption{Average velocity and maximum velocity of shock diffraction over a cylinder.}
    \label{tab:velocity}
    \scriptsize
    \begin{tabular}{lcc}
        \hline
        & No-slip wall & Slip wall \\
        \hline
        Average velocity of flow field  & $2.00$  & $2.01$ \\
        Maximum velocity in flow field  & $3.83$ & $4.41$ \\
        \hline
    \end{tabular}
\end{table}

According to these qualitative and quantitative results, the developed immersed boundary method of this paper retains a sharp interface and is able to effectively alleviate unphysical flux over physical boundaries when grid resolution is improved.

\section{Conclusions} \label{sec:conclusion}

A sharp interface immersed boundary method is developed. This method enforces Dirichlet, Neumann, Robin, and Cauchy boundary conditions in a straightforward and consistent manner and is able to provide efficient and robust boundary treatment for numerically solving flow with arbitrarily irregular and changing geometries while maintaining accuracy. The effectiveness of the proposed method is conformed by numerical experiments concerning flow of two and three space dimensions, stationary and moving objects, convex and concave geometries, no-slip and slip wall boundary conditions.

\appendix

\section{A sample code for boundary treatment implementation} \label{app:source_code}

\subsection{Ghost flow reconstruction function} \label{app:reconstruction}

\begin{footnotesize}
\begin{verbatim}
/* 
 * Flow reconstruction of Field_G[N] for N field variables at ghost node G.
 */
/* pre-estimate Field_I[N] in domain of dependence of image point I. */
compute weightedSum[N] and sumOfWeights by Appendix A.2.;
for (int n = 0; n < N; ++n) {
    Field_Istar[n] = weightedSum[n] / sumOfWeights;
}
/* enforce physical boundary conditions to determine Field_O[N] */
Field_O[0] = C[0] * Field_Istar[0] + R.R.H.S.[0];
.
.
.
Field_O[N-1] = C[N-1] * Field_Istar[N-1] + R.R.H.S.[N-1];
/* correction step to solve Field_I[N] */
for (int n = 0; n < N; ++n) {
    Field_I[n] = (weightedSum[n] + Field_O[n] * weight_O) /
                 (sumOfWeights + weight_O);
}
/* apply the method of images to construct Field_G[N] at ghost node G */
for (int n = 0; n < N; ++n) {
    Field_G[n] = 2 * Field_O[n] - Field_I[n];
}
\end{verbatim}
\end{footnotesize}

\subsection{A search function for inverse distance weighting} \label{app:inverse_weighting}

\begin{footnotesize}
\begin{verbatim}
/* 
 * Search fluid nodes around node(kI, jI, iI) in the domain of dependence R
 * to apply inverse distance weighting. node(kI, jI, iI) is a computational
 * node whose node coordinates are derived from the corresponding spatial 
 * coordinates of a image point I(xI, yI, zI).
 */
for (int kh = -R; kh <= R; ++kh) {
    for (int jh = -R; jh <= R; ++jh) {
        for (int ih = -R; ih <= R; ++ih) {
            if (Flag[kI+kh][jI+jh][iI+ih] != FLUID) { /* not a fluid node */
                continue;
            }
            /* a valid stencil node, apply inverse distance weighting */
            compute weight for node(kI+kh, jI+jh, iI+ih) to I(xI, yI, zI);
            sumOfWeights = sumOfWeights + weight;
            for (int n = 0; n < N; ++n) {
                weightedSum[n] = weightedSum[n] +
                                 Field[kI+kh][jI+jh][iI+ih][n] * weight;
            }
        }
    }
}
\end{verbatim}
\end{footnotesize}

\section{Fluid-solid coupling} \label{app:fluid_solid_coupling}

Currently, a simple fluid-solid coupling pattern is used. The interactions between fluid and solid are explicitly coupled by applying Strang splitting \cite{strang1964accurate, strang1968construction} for physical process splitting as the following:

\begin{enumerate}
    \item Evolve particle dynamics for $\frac{1}{2}\Delta t$ 
    \item Evolve fluid dynamics for $\Delta t$
    \item Evolve particle dynamics for $\frac{1}{2}\Delta t$ 
\end{enumerate}

The evolution of particle dynamics can be expanded as the following procedures:

\begin{enumerate}
    \item Integrate surface forces.
    \item Update the spatial distribution of particles by particle dynamics models.
    \item Detect ghost nodes that fall out the regions of their corresponding particles. These ghost nodes now change to fluid nodes, of which values of flow variables can be reconstructed by inverse distance weighting.
    \item Re-mesh computational domain, and apply boundary treatment.
\end{enumerate}

\section*{Acknowledgements}

Financial support of this work was provided by Natural Sciences and Engineering Research Council of Canada (NSERC) and Defence Research and Development Canada (DRDC). This work was made possible by the facilities of the Shared Hierarchical Academic Research Computing Network (SHARCNET:www.sharcnet.ca) and Compute/Calcul Canada. The first author of this paper is grateful to Prof. Deliang Zhang for introducing Computational Fluid Dynamics and is thankful to Dr. Deyong Wen for discussions of flow visualization.

\bibliography{ref}

\end{document}

arxiv_huangrui_mo.bbl
\begin{thebibliography}{57}
\providecommand{\natexlab}[1]{#1}
\providecommand{\url}[1]{\texttt{#1}}
\providecommand{\urlprefix}{URL }
\expandafter\ifx\csname urlstyle\endcsname\relax
  \providecommand{\doi}[1]{doi:\discretionary{}{}{}#1}\else
  \providecommand{\doi}[1]{doi:\discretionary{}{}{}\begingroup
  \urlstyle{rm}\url{#1}\endgroup}\fi
\providecommand{\bibinfo}[2]{#2}

\bibitem[{Frost et~al.(2012)Frost, Gr{\'e}goire, Petel, Goroshin, and
  Zhang}]{frost2012particle}
\bibinfo{author}{D.~L. Frost}, \bibinfo{author}{Y.~Gr{\'e}goire},
  \bibinfo{author}{O.~Petel}, \bibinfo{author}{S.~Goroshin},
  \bibinfo{author}{F.~Zhang}, \bibinfo{title}{Particle jet formation during
  explosive dispersal of solid particles}, \bibinfo{journal}{Phys. Fluids}
  \bibinfo{volume}{24}~(\bibinfo{number}{9}) (\bibinfo{year}{2012})
  \bibinfo{pages}{091109}.

\bibitem[{Zhang et~al.(2014)Zhang, Ripley, Yoshinaka, Findlay, Anderson, and
  von Rosen}]{zhang2014large}
\bibinfo{author}{F.~Zhang}, \bibinfo{author}{R.~Ripley},
  \bibinfo{author}{A.~Yoshinaka}, \bibinfo{author}{C.~Findlay},
  \bibinfo{author}{J.~Anderson}, \bibinfo{author}{B.~von Rosen},
  \bibinfo{title}{Large-scale spray detonation and related particle jetting
  instability phenomenon}, \bibinfo{journal}{Shock Waves}
  (\bibinfo{year}{2014}) \bibinfo{pages}{1--16}.

\bibitem[{Kedrinskiy(2009)}]{kedrinskiy2009hydrodynamic}
\bibinfo{author}{V.~Kedrinskiy}, \bibinfo{title}{Hydrodynamic aspects of
  explosive eruptions of volcanoes: simulation problems},
  \bibinfo{journal}{Shock Waves} \bibinfo{volume}{18}~(\bibinfo{number}{6})
  (\bibinfo{year}{2009}) \bibinfo{pages}{451--464}.

\bibitem[{Lohse et~al.(2004)Lohse, Bergmann, and Mikkelsen}]{lohse2004impact}
\bibinfo{author}{D.~Lohse}, \bibinfo{author}{R.~Bergmann},
  \bibinfo{author}{R.~Mikkelsen}, \bibinfo{title}{Impact on soft sand: void
  collapse and jet formation}, \bibinfo{journal}{Phys. Rev. Lett.}
  \bibinfo{volume}{93}~(\bibinfo{number}{19}) (\bibinfo{year}{2004})
  \bibinfo{pages}{198003}.

\bibitem[{Pacheco-V{\'a}zquez and Ruiz-Su{\'a}rez(2011)}]{pacheco2011impact}
\bibinfo{author}{F.~Pacheco-V{\'a}zquez}, \bibinfo{author}{J.~Ruiz-Su{\'a}rez},
  \bibinfo{title}{Impact craters in granular media: grains against grains},
  \bibinfo{journal}{Phys. Rev. Lett.}
  \bibinfo{volume}{107}~(\bibinfo{number}{21}) (\bibinfo{year}{2011})
  \bibinfo{pages}{218001}.

\bibitem[{Zhang(2009)}]{zhang2009shock}
\bibinfo{author}{F.~Zhang}, \bibinfo{title}{Shock Wave Science and Technology
  Reference Library: Heterogeneous Detonation}, vol.~\bibinfo{volume}{4},
  \bibinfo{publisher}{Springer Science \& Business Media},
  \bibinfo{year}{2009}.

\bibitem[{Rodriguez et~al.(2013)Rodriguez, Saurel, Jourdan, and
  Houas}]{rodriguez2013solid}
\bibinfo{author}{V.~Rodriguez}, \bibinfo{author}{R.~Saurel},
  \bibinfo{author}{G.~Jourdan}, \bibinfo{author}{L.~Houas},
  \bibinfo{title}{Solid-particle jet formation under shock-wave acceleration},
  \bibinfo{journal}{Phys. Rev. E} \bibinfo{volume}{88}~(\bibinfo{number}{6})
  (\bibinfo{year}{2013}) \bibinfo{pages}{063011}.

\bibitem[{Liu et~al.(1995)Liu, Nagel, Schecter, Coppersmith, Majumdar, Narayan,
  and Witten}]{liu1995force}
\bibinfo{author}{C.~Liu}, \bibinfo{author}{S.~R. Nagel},
  \bibinfo{author}{D.~Schecter}, \bibinfo{author}{S.~Coppersmith},
  \bibinfo{author}{S.~Majumdar}, \bibinfo{author}{O.~Narayan},
  \bibinfo{author}{T.~Witten}, \bibinfo{title}{Force fluctuations in bead
  packs}, \bibinfo{journal}{SCIENCE}
  \bibinfo{volume}{269}~(\bibinfo{number}{5223}) (\bibinfo{year}{1995})
  \bibinfo{pages}{513--513}.

\bibitem[{Jaeger et~al.(1996)Jaeger, Nagel, and Behringer}]{jaeger1996granular}
\bibinfo{author}{H.~M. Jaeger}, \bibinfo{author}{S.~R. Nagel},
  \bibinfo{author}{R.~P. Behringer}, \bibinfo{title}{Granular solids, liquids,
  and gases}, \bibinfo{journal}{Rev. Mod. Phys.}
  \bibinfo{volume}{68}~(\bibinfo{number}{4}) (\bibinfo{year}{1996})
  \bibinfo{pages}{1259}.

\bibitem[{Peskin(1972)}]{peskin1972flow}
\bibinfo{author}{C.~S. Peskin}, \bibinfo{title}{Flow patterns around heart
  valves: a numerical method}, \bibinfo{journal}{J. Comput. Phys.}
  \bibinfo{volume}{10}~(\bibinfo{number}{2}) (\bibinfo{year}{1972})
  \bibinfo{pages}{252--271}.

\bibitem[{Jung and Peskin(2001)}]{jung2001two}
\bibinfo{author}{E.~Jung}, \bibinfo{author}{C.~S. Peskin},
  \bibinfo{title}{Two-dimensional simulations of valveless pumping using the
  immersed boundary method}, \bibinfo{journal}{SIAM J. Sci. Comput.}
  \bibinfo{volume}{23}~(\bibinfo{number}{1}) (\bibinfo{year}{2001})
  \bibinfo{pages}{19--45}.

\bibitem[{Peskin(2002)}]{peskin2002immersed}
\bibinfo{author}{C.~S. Peskin}, \bibinfo{title}{The immersed boundary method},
  \bibinfo{journal}{Acta Numer.} \bibinfo{volume}{11} (\bibinfo{year}{2002})
  \bibinfo{pages}{479--517}.

\bibitem[{Mittal and Iaccarino(2005)}]{mittal2005immersed}
\bibinfo{author}{R.~Mittal}, \bibinfo{author}{G.~Iaccarino},
  \bibinfo{title}{Immersed boundary methods}, \bibinfo{journal}{Annu. Rev.
  Fluid Mech.} \bibinfo{volume}{37} (\bibinfo{year}{2005})
  \bibinfo{pages}{239--261}.

\bibitem[{Fadlun et~al.(2000)Fadlun, Verzicco, Orlandi, and
  Mohd-Yusof}]{fadlun2000combined}
\bibinfo{author}{E.~A. Fadlun}, \bibinfo{author}{R.~Verzicco},
  \bibinfo{author}{P.~Orlandi}, \bibinfo{author}{J.~Mohd-Yusof},
  \bibinfo{title}{Combined immersed-boundary finite-difference methods for
  three-dimensional complex flow simulations}, \bibinfo{journal}{J. Comput.
  Phys.} \bibinfo{volume}{161}~(\bibinfo{number}{1}) (\bibinfo{year}{2000})
  \bibinfo{pages}{35--60}.

\bibitem[{Tseng and Ferziger(2003)}]{tseng2003ghost}
\bibinfo{author}{Y.-H. Tseng}, \bibinfo{author}{J.~H. Ferziger},
  \bibinfo{title}{A ghost-cell immersed boundary method for flow in complex
  geometry}, \bibinfo{journal}{J. Comput. Phys.}
  \bibinfo{volume}{192}~(\bibinfo{number}{2}) (\bibinfo{year}{2003})
  \bibinfo{pages}{593--623}.

\bibitem[{Uhlmann(2005)}]{uhlmann2005immersed}
\bibinfo{author}{M.~Uhlmann}, \bibinfo{title}{An immersed boundary method with
  direct forcing for the simulation of particulate flows}, \bibinfo{journal}{J.
  Comput. Phys.} \bibinfo{volume}{209}~(\bibinfo{number}{2})
  (\bibinfo{year}{2005}) \bibinfo{pages}{448--476}.

\bibitem[{Mori and Peskin(2008)}]{mori2008implicit}
\bibinfo{author}{Y.~Mori}, \bibinfo{author}{C.~S. Peskin},
  \bibinfo{title}{Implicit second-order immersed boundary methods with boundary
  mass}, \bibinfo{journal}{Comput. Methods Appl. Mech. Engrg.}
  \bibinfo{volume}{197}~(\bibinfo{number}{25}) (\bibinfo{year}{2008})
  \bibinfo{pages}{2049--2067}.

\bibitem[{Kapahi et~al.(2013)Kapahi, Sambasivan, and
  Udaykumar}]{kapahi2013three}
\bibinfo{author}{A.~Kapahi}, \bibinfo{author}{S.~Sambasivan},
  \bibinfo{author}{H.~Udaykumar}, \bibinfo{title}{A three-dimensional sharp
  interface {C}artesian grid method for solving high speed multi-material
  impact, penetration and fragmentation problems}, \bibinfo{journal}{J. Comput.
  Phys.} \bibinfo{volume}{241} (\bibinfo{year}{2013})
  \bibinfo{pages}{308--332}.

\bibitem[{Yang and Stern(2015)}]{yang2015non}
\bibinfo{author}{J.~Yang}, \bibinfo{author}{F.~Stern}, \bibinfo{title}{A
  non-iterative direct forcing immersed boundary method for strongly-coupled
  fluid-solid interactions}, \bibinfo{journal}{J. Comput. Phys.}
  \bibinfo{volume}{295} (\bibinfo{year}{2015}) \bibinfo{pages}{779--804}.

\bibitem[{Schwarz et~al.(2016)Schwarz, Kempe, and
  Fr{\"o}hlich}]{schwarz2016immersed}
\bibinfo{author}{S.~Schwarz}, \bibinfo{author}{T.~Kempe},
  \bibinfo{author}{J.~Fr{\"o}hlich}, \bibinfo{title}{An immersed boundary
  method for the simulation of bubbles with varying shape},
  \bibinfo{journal}{J. Comput. Phys.} .

\bibitem[{Mohd-Yusof(1997)}]{yusof1997combined}
\bibinfo{author}{J.~Mohd-Yusof}, \bibinfo{title}{Combined
  immersed-boundary/{B}-spline methods for simulations of flow in complex
  geometries}, \bibinfo{journal}{CTR-Annual Research Briefs-Stanford Univ./NASA
  Ames} .

\bibitem[{Balaras(2004)}]{balaras2004modeling}
\bibinfo{author}{E.~Balaras}, \bibinfo{title}{Modeling complex boundaries using
  an external force field on fixed {C}artesian grids in large-eddy
  simulations}, \bibinfo{journal}{Comput. Fluids}
  \bibinfo{volume}{33}~(\bibinfo{number}{3}) (\bibinfo{year}{2004})
  \bibinfo{pages}{375--404}.

\bibitem[{Fedkiw et~al.(1999)Fedkiw, Aslam, Merriman, and
  Osher}]{fedkiw1999non}
\bibinfo{author}{R.~P. Fedkiw}, \bibinfo{author}{T.~Aslam},
  \bibinfo{author}{B.~Merriman}, \bibinfo{author}{S.~Osher}, \bibinfo{title}{A
  non-oscillatory {E}ulerian approach to interfaces in multimaterial flows (the
  ghost fluid method)}, \bibinfo{journal}{J. Comput. Phys.}
  \bibinfo{volume}{152}~(\bibinfo{number}{2}) (\bibinfo{year}{1999})
  \bibinfo{pages}{457--492}.

\bibitem[{Fedkiw(2002)}]{fedkiw2002coupling}
\bibinfo{author}{R.~P. Fedkiw}, \bibinfo{title}{Coupling an {E}ulerian fluid
  calculation to a {L}agrangian solid calculation with the ghost fluid method},
  \bibinfo{journal}{J. Comput. Phys.}
  \bibinfo{volume}{175}~(\bibinfo{number}{1}) (\bibinfo{year}{2002})
  \bibinfo{pages}{200--224}.

\bibitem[{Iaccarino and Verzicco(2003)}]{iaccarino2003immersed}
\bibinfo{author}{G.~Iaccarino}, \bibinfo{author}{R.~Verzicco},
  \bibinfo{title}{Immersed boundary technique for turbulent flow simulations},
  \bibinfo{journal}{Appl. Mech. Rev.}
  \bibinfo{volume}{56}~(\bibinfo{number}{3}) (\bibinfo{year}{2003})
  \bibinfo{pages}{331--347}.

\bibitem[{Sotiropoulos and Yang(2014)}]{sotiropoulos2014immersed}
\bibinfo{author}{F.~Sotiropoulos}, \bibinfo{author}{X.~Yang},
  \bibinfo{title}{Immersed boundary methods for simulating fluid--structure
  interaction}, \bibinfo{journal}{Prog. Aerosp. Sci.} \bibinfo{volume}{65}
  (\bibinfo{year}{2014}) \bibinfo{pages}{1--21}.

\bibitem[{Gao et~al.(2007)Gao, Tseng, and Lu}]{gao2007improved}
\bibinfo{author}{T.~Gao}, \bibinfo{author}{Y.-H. Tseng}, \bibinfo{author}{X.-Y.
  Lu}, \bibinfo{title}{An improved hybrid {C}artesian/immersed boundary method
  for fluid-solid flows}, \bibinfo{journal}{Int. J. Numer. Meth. Fluids}
  \bibinfo{volume}{55}~(\bibinfo{number}{12}) (\bibinfo{year}{2007})
  \bibinfo{pages}{1189--1211}.

\bibitem[{Gibou et~al.(2013)Gibou, Min, and Fedkiw}]{gibou2013high}
\bibinfo{author}{F.~Gibou}, \bibinfo{author}{C.~Min}, \bibinfo{author}{R.~P.
  Fedkiw}, \bibinfo{title}{High resolution sharp computational methods for
  elliptic and parabolic problems in complex geometries}, \bibinfo{journal}{J.
  Sci. Comput.} \bibinfo{volume}{54}~(\bibinfo{number}{2-3})
  (\bibinfo{year}{2013}) \bibinfo{pages}{369--413}.

\bibitem[{Crockett et~al.(2011)Crockett, Colella, and
  Graves}]{crockett2011cartesian}
\bibinfo{author}{R.~Crockett}, \bibinfo{author}{P.~Colella},
  \bibinfo{author}{D.~T. Graves}, \bibinfo{title}{A {C}artesian grid embedded
  boundary method for solving the {P}oisson and heat equations with
  discontinuous coefficients in three dimensions}, \bibinfo{journal}{J. Comput.
  Phys.} \bibinfo{volume}{230}~(\bibinfo{number}{7}) (\bibinfo{year}{2011})
  \bibinfo{pages}{2451--2469}.

\bibitem[{Kempe et~al.(2015)Kempe, Lennartz, Schwarz, and
  Fr{\"o}hlich}]{kempe2015imposing}
\bibinfo{author}{T.~Kempe}, \bibinfo{author}{M.~Lennartz},
  \bibinfo{author}{S.~Schwarz}, \bibinfo{author}{J.~Fr{\"o}hlich},
  \bibinfo{title}{Imposing the free-slip condition with a continuous forcing
  immersed boundary method}, \bibinfo{journal}{J. Comput. Phys.}
  \bibinfo{volume}{282} (\bibinfo{year}{2015}) \bibinfo{pages}{183--209}.

\bibitem[{Chu and Sereny(1974)}]{chu1974boundary}
\bibinfo{author}{C.~Chu}, \bibinfo{author}{A.~Sereny}, \bibinfo{title}{Boundary
  conditions in finite difference fluid dynamic codes}, \bibinfo{journal}{J.
  Comput. Phys.} \bibinfo{volume}{15}~(\bibinfo{number}{4})
  (\bibinfo{year}{1974}) \bibinfo{pages}{476--491}.

\bibitem[{Colella(1990)}]{colella1990multidimensional}
\bibinfo{author}{P.~Colella}, \bibinfo{title}{Multidimensional upwind methods
  for hyperbolic conservation laws}, \bibinfo{journal}{J. Comput. Phys.}
  \bibinfo{volume}{87}~(\bibinfo{number}{1}) (\bibinfo{year}{1990})
  \bibinfo{pages}{171--200}.

\bibitem[{Mark and van Wachem(2008)}]{mark2008derivation}
\bibinfo{author}{A.~Mark}, \bibinfo{author}{B.~G. van Wachem},
  \bibinfo{title}{Derivation and validation of a novel implicit second-order
  accurate immersed boundary method}, \bibinfo{journal}{J. Comput. Phys.}
  \bibinfo{volume}{227}~(\bibinfo{number}{13}) (\bibinfo{year}{2008})
  \bibinfo{pages}{6660--6680}.

\bibitem[{Seo and Mittal(2011)}]{seo2011sharp}
\bibinfo{author}{J.~H. Seo}, \bibinfo{author}{R.~Mittal}, \bibinfo{title}{A
  sharp-interface immersed boundary method with improved mass conservation and
  reduced spurious pressure oscillations}, \bibinfo{journal}{J. Comput. Phys.}
  \bibinfo{volume}{230}~(\bibinfo{number}{19}) (\bibinfo{year}{2011})
  \bibinfo{pages}{7347--7363}.

\bibitem[{Gilmanov et~al.(2003)Gilmanov, Sotiropoulos, and
  Balaras}]{gilmanov2003general}
\bibinfo{author}{A.~Gilmanov}, \bibinfo{author}{F.~Sotiropoulos},
  \bibinfo{author}{E.~Balaras}, \bibinfo{title}{A general reconstruction
  algorithm for simulating flows with complex 3{D} immersed boundaries on
  {C}artesian grids}, \bibinfo{journal}{J. Comput. Phys.}
  \bibinfo{volume}{191}~(\bibinfo{number}{2}) (\bibinfo{year}{2003})
  \bibinfo{pages}{660--669}.

\bibitem[{Gilmanov and Sotiropoulos(2005)}]{gilmanov2005hybrid}
\bibinfo{author}{A.~Gilmanov}, \bibinfo{author}{F.~Sotiropoulos},
  \bibinfo{title}{A hybrid {C}artesian/immersed boundary method for simulating
  flows with 3{D}, geometrically complex, moving bodies}, \bibinfo{journal}{J.
  Comput. Phys.} \bibinfo{volume}{207}~(\bibinfo{number}{2})
  (\bibinfo{year}{2005}) \bibinfo{pages}{457--492}.

\bibitem[{Shepard(1968)}]{shepard1968two}
\bibinfo{author}{D.~Shepard}, \bibinfo{title}{A two-dimensional interpolation
  function for irregularly-spaced data}, in: \bibinfo{booktitle}{Proc. 23rd
  Nat. Conf. ACM}, \bibinfo{organization}{ACM}, \bibinfo{pages}{517--524},
  \bibinfo{year}{1968}.

\bibitem[{Junkins et~al.(1973)Junkins, Miller, and
  Jancaitis}]{junkins1973weighting}
\bibinfo{author}{J.~L. Junkins}, \bibinfo{author}{G.~W. Miller},
  \bibinfo{author}{J.~R. Jancaitis}, \bibinfo{title}{A weighting function
  approach to modeling of irregular surfaces}, \bibinfo{journal}{J. Geophys.
  Res.} \bibinfo{volume}{78}~(\bibinfo{number}{11}) (\bibinfo{year}{1973})
  \bibinfo{pages}{1794--1803}.

\bibitem[{Franke(1982)}]{franke1982scattered}
\bibinfo{author}{R.~Franke}, \bibinfo{title}{Scattered data interpolation:
  tests of some methods}, \bibinfo{journal}{Math. Comput.}
  \bibinfo{volume}{38}~(\bibinfo{number}{157}) (\bibinfo{year}{1982})
  \bibinfo{pages}{181--200}.

\bibitem[{Anderson and Wendt(1995)}]{anderson1995computational}
\bibinfo{author}{J.~D. Anderson}, \bibinfo{author}{J.~Wendt},
  \bibinfo{title}{Computational Fluid Dynamics}, vol. \bibinfo{volume}{206},
  \bibinfo{publisher}{Springer}, \bibinfo{year}{1995}.

\bibitem[{Shu and Osher(1988)}]{shu1988efficient}
\bibinfo{author}{C.-W. Shu}, \bibinfo{author}{S.~Osher},
  \bibinfo{title}{Efficient implementation of essentially non-oscillatory
  shock-capturing schemes}, \bibinfo{journal}{J. Comput. Phys.}
  \bibinfo{volume}{77}~(\bibinfo{number}{2}) (\bibinfo{year}{1988})
  \bibinfo{pages}{439--471}.

\bibitem[{Gottlieb et~al.(2001)Gottlieb, Shu, and Tadmor}]{gottlieb2001strong}
\bibinfo{author}{S.~Gottlieb}, \bibinfo{author}{C.-W. Shu},
  \bibinfo{author}{E.~Tadmor}, \bibinfo{title}{Strong stability-preserving
  high-order time discretization methods}, \bibinfo{journal}{SIAM Rev.}
  \bibinfo{volume}{43}~(\bibinfo{number}{1}) (\bibinfo{year}{2001})
  \bibinfo{pages}{89--112}.

\bibitem[{Harten(1984)}]{harten1984class}
\bibinfo{author}{A.~Harten}, \bibinfo{title}{On a class of high resolution
  total-variation-stable finite-difference schemes}, \bibinfo{journal}{SIAM J.
  Numer. Anal.} \bibinfo{volume}{21}~(\bibinfo{number}{1})
  (\bibinfo{year}{1984}) \bibinfo{pages}{1--23}.

\bibitem[{Jiang and Shu(1996)}]{jiang1996efficient}
\bibinfo{author}{G.-S. Jiang}, \bibinfo{author}{C.-W. Shu},
  \bibinfo{title}{Efficient implementation of weighted {ENO} schemes},
  \bibinfo{journal}{J. Comput. Phys.}
  \bibinfo{volume}{126}~(\bibinfo{number}{1}).

\bibitem[{Ferziger and Peri{\'c}(2002)}]{ferziger2002computational}
\bibinfo{author}{J.~H. Ferziger}, \bibinfo{author}{M.~Peri{\'c}},
  \bibinfo{title}{Computational Methods for Fluid Dynamics},
  vol.~\bibinfo{volume}{3}, \bibinfo{publisher}{Springer Berlin},
  \bibinfo{year}{2002}.

\bibitem[{Strang(1964)}]{strang1964accurate}
\bibinfo{author}{G.~Strang}, \bibinfo{title}{Accurate partial difference
  methods}, \bibinfo{journal}{Numer. Math.}
  \bibinfo{volume}{6}~(\bibinfo{number}{1}) (\bibinfo{year}{1964})
  \bibinfo{pages}{37--46}.

\bibitem[{Strang(1968)}]{strang1968construction}
\bibinfo{author}{G.~Strang}, \bibinfo{title}{On the construction and comparison
  of difference schemes}, \bibinfo{journal}{SIAM J. Numer. Anal.}
  \bibinfo{volume}{5}~(\bibinfo{number}{3}) (\bibinfo{year}{1968})
  \bibinfo{pages}{506--517}.

\bibitem[{Woodward and Colella(1984)}]{woodward1984numerical}
\bibinfo{author}{P.~Woodward}, \bibinfo{author}{P.~Colella},
  \bibinfo{title}{The numerical simulation of two-dimensional fluid flow with
  strong shocks}, \bibinfo{journal}{J. Comput. Phys.}
  \bibinfo{volume}{54}~(\bibinfo{number}{1}) (\bibinfo{year}{1984})
  \bibinfo{pages}{115--173}.

\bibitem[{Bryson and Gross(1961)}]{bryson1961diffraction}
\bibinfo{author}{A.~Bryson}, \bibinfo{author}{R.~Gross},
  \bibinfo{title}{Diffraction of strong shocks by cones, cylinders, and
  spheres}, \bibinfo{journal}{J. Fluid Mech.}
  \bibinfo{volume}{10}~(\bibinfo{number}{01}) (\bibinfo{year}{1961})
  \bibinfo{pages}{1--16}.

\bibitem[{Bazhenova et~al.(1984)Bazhenova, Gvozdeva, and
  Nettleton}]{bazhenova1984unsteady}
\bibinfo{author}{T.~Bazhenova}, \bibinfo{author}{L.~Gvozdeva},
  \bibinfo{author}{M.~Nettleton}, \bibinfo{title}{Unsteady interactions of
  shock waves}, \bibinfo{journal}{Prog. Aerosp. Sci.} \bibinfo{volume}{21}
  (\bibinfo{year}{1984}) \bibinfo{pages}{249--331}.

\bibitem[{Kaca(1988)}]{kaca1988interferometric}
\bibinfo{author}{J.~Kaca}, \bibinfo{title}{An interferometric investigation of
  the diffraction of a planar shock wave over a semicircular cylinder},
  \bibinfo{type}{Tech. Rep.}, \bibinfo{institution}{University of Toronto
  Institute for Aerospace Studies}, \bibinfo{note}{269}, \bibinfo{year}{1988}.

\bibitem[{Ripley et~al.(2006)Ripley, Lien, and
  Yovanovich}]{ripley2006numerical}
\bibinfo{author}{R.~Ripley}, \bibinfo{author}{F.-S. Lien},
  \bibinfo{author}{M.~Yovanovich}, \bibinfo{title}{Numerical simulation of
  shock diffraction on unstructured meshes}, \bibinfo{journal}{Comput. Fluids}
  \bibinfo{volume}{35}~(\bibinfo{number}{10}) (\bibinfo{year}{2006})
  \bibinfo{pages}{1420--1431}.

\bibitem[{Sambasivan and Udaykumar(2009)}]{sambasivan2009ghostb}
\bibinfo{author}{S.~K. Sambasivan}, \bibinfo{author}{H.~Udaykumar},
  \bibinfo{title}{Ghost fluid method for strong shock interactions Part 2:
  immersed solid boundaries}, \bibinfo{journal}{AIAA J.}
  \bibinfo{volume}{47}~(\bibinfo{number}{12}) (\bibinfo{year}{2009})
  \bibinfo{pages}{2923--2937}.

\bibitem[{Oberkampf and Trucano(2002)}]{oberkampf2002verification}
\bibinfo{author}{W.~L. Oberkampf}, \bibinfo{author}{T.~G. Trucano},
  \bibinfo{title}{Verification and validation in computational fluid dynamics},
  \bibinfo{journal}{Prog. Aerosp. Sci.}
  \bibinfo{volume}{38}~(\bibinfo{number}{3}) (\bibinfo{year}{2002})
  \bibinfo{pages}{209--272}.

\bibitem[{Roache(1998)}]{roache1998verification}
\bibinfo{author}{P.~J. Roache}, \bibinfo{title}{Verification and validation in
  computational science and engineering}, \bibinfo{publisher}{Hermosa},
  \bibinfo{year}{1998}.

\bibitem[{Ji et~al.(2010)Ji, Lien, and Yee}]{ji2010numerical}
\bibinfo{author}{H.~Ji}, \bibinfo{author}{F.-S. Lien},
  \bibinfo{author}{E.~Yee}, \bibinfo{title}{Numerical simulation of detonation
  using an adaptive {C}artesian cut-cell method combined with a cell-merging
  technique}, \bibinfo{journal}{Comput. Fluids}
  \bibinfo{volume}{39}~(\bibinfo{number}{6}) (\bibinfo{year}{2010})
  \bibinfo{pages}{1041--1057}.

\bibitem[{Ripley and Zhang(2014)}]{ripley2014jetting}
\bibinfo{author}{R.~Ripley}, \bibinfo{author}{F.~Zhang},
  \bibinfo{title}{Jetting instability mechanisms of particles from explosive
  dispersal}, in: \bibinfo{booktitle}{J. Phys. Conf. Ser.}, vol.
  \bibinfo{volume}{500}, \bibinfo{organization}{IOP Publishing},
  \bibinfo{pages}{152012--152019}, \bibinfo{year}{2014}.

\end{thebibliography}
