Volume penalization for inhomogeneous Neumann boundary conditions modeling scalar flux in complicated geometry

TL;DR

This paper introduces a volume penalization method for inhomogeneous Neumann boundary conditions, enabling modeling of scalar flux in complex geometries within simple Cartesian domains, and validates it through numerical experiments.

Contribution

The paper generalizes existing flux-based volume penalization to inhomogeneous Neumann conditions, facilitating easier simulation of scalar flux in complex geometries.

Findings

The method accurately models scalar flux in complex geometries.

Convergence analysis shows reliable numerical behavior.

Application to convection demonstrates practical effectiveness.

Abstract

We develop a volume penalization method for inhomogeneous Neumann boundary conditions, generalizing the flux-based volume penalization method for homogeneous Neumann boundary condition proposed by Kadoch et al. [J. Comput. Phys. 231 (2012) 4365]. The generalized method allows us to model scalar flux through walls in geometries of complex shape using simple, e.g. Cartesian, domains for solving the governing equations. We examine the properties of the method, by considering a one-dimensional Poisson equation with different Neumann boundary conditions. The penalized Laplace operator is discretized by second order central finite-differences and interpolation. The discretization and penalization errors are thus assessed for several test problems. Convergence properties of the discretized operator and the solution of the penalized equation are analyzed. The generalized method is then applied to an advection-diffusion equation coupled with the Navier-Stokes equations in an annular domain which is immersed in a square domain. The application is verified by numerical simulation of steady free convection in a concentric annulus heated through the inner cylinder surface using an extended square domain.

Figures

Equations

$$-d_x^{2}w(x)=f(x)\text{mbox}forx\in {\Omega }_f$$

$$f(x)=m^2\cos mx.$$

$$d_{x}w(0)=d_{x}w(\pi )=\alpha ,$$

$$w(x)=\cos mx+\alpha x+A_0,$$

$$-d_x(\theta d_{x}v+\chi \alpha )=(1-\chi )f,$$

$$\theta (x)=1-\chi (x)+\eta \chi (x),$$

$$\displaystyle\chi(x)=\left\{\begin{array}[]{lll}0&{\rm for}&0<x<\pi,\\ 1/2&{\rm at}&x=0,\pi,\\ 1&{\rm for}&\pi<x<2\pi.\end{array}\right.$$

$$d_{x}v=\eta d_{x}v+\alpha \mathrm{at}x=0(=2\pi ),\pi .$$

$$v(x)=\left\{\begin{array}[]{ll}\cos{mx}+A_{1}x+A_{2}&{\rm for}\quad 0\leq x\leq\pi,\\ B_{1}x+B_{2}&{\rm for}\quad\pi\leq x<2\pi.\end{array}\right.$$

$$A_1=\frac{\alpha }{1+\eta }+\eta \frac{1-(-1{)}^m}{\pi (1+\eta )},$$

$$A_2=2\pi \left\{-\frac{\alpha }{1+\eta }+\frac{1-(-1{)}^m}{\pi (1+\eta )}\right\}-1+B_2,$$

$$B_1=-\frac{\alpha }{1+\eta }+\frac{1-(-1{)}^m}{\pi (1+\eta )}.$$

$$B_2=\frac{3\pi }{2}\frac{\alpha }{1+\eta }-\frac{1-(-1{)}^m}{1+\eta }\left(\frac{\eta }{2}+2\right)+1.$$

$${\hat{v}}_k=\frac{1}{2\pi }{\int }_0^{2\pi }v(x)\exp (-\mathrm{i}kx)dx.$$

$$\displaystyle{\hat{v}}_{k}=\left\{\begin{array}[]{ll}\displaystyle{\frac{{\mathrm{i}}}{\pi}\frac{m^{2}}{k(m^{2}-k^{2})}}&{\rm if}\ k\ {\rm even},\\ &\\ \displaystyle{\frac{2}{\pi^{2}m^{2}}\frac{1-\eta}{1+\eta}-\frac{2\alpha}{\pi m^{2}}\frac{1}{1+\eta}+\frac{1}{4}}&{\rm if}\ k\ {\rm odd\ and}\ k=\pm m,\\ &\\ \displaystyle{\frac{2}{\pi^{2}k^{2}}\frac{1-\eta}{1+\eta}-\frac{2\alpha}{\pi k^{2}}\frac{1}{1+\eta}}&{\rm if}\ k\ {\rm odd\ and}\ k\neq\pm m,\end{array}\right.$$

$$\displaystyle{\hat{v}}_{k}=\left\{\begin{array}[]{ll}\displaystyle{-\frac{2\alpha}{\pi k^{2}}\frac{1}{1+\eta}+\frac{{\mathrm{i}}}{\pi}\frac{m^{2}}{k(m^{2}-k^{2})}}&{\rm if}\ k\ {\rm odd},\\ &\\ \displaystyle{\frac{1}{4}}&{\rm if}\ k\ {\rm even\ and}\ k=\pm m,\\ &\\ 0&{\rm if}\ k\ {\rm even\ and}\ k\neq\pm m.\end{array}\right.$$

$$d_x\varphi (x)=\frac{\varphi (x+h/2)-\varphi (x-h/2)}{h}+O(h^2),$$

$$\varphi (x)=\frac{\varphi (x+h/2)+\varphi (x-h/2)}{2}+O(h^2),$$

$$f(x)=m^2\sin (mx)$$

$$d_{x}w(0)=\alpha +m,\mathrm{and}{d}_{x}w(\pi )=\alpha -m.$$

$$w(x)=\sin (mx)+\alpha x-\frac{2}{m\pi }-\frac{\pi }{2}\alpha .$$

$$-d_x(\theta d_{x}v+\beta \chi )=(1-\chi )f-\chi d_x\beta ,$$

$$d_{x}v=\eta d_{x}v+\beta \mathrm{at}x=0(=2\pi ),\pi .$$

$$\displaystyle v(x)=\left\{\begin{array}[]{ll}\sin(mx)+A_{1}^{\prime}x+A_{2}^{\prime}&{\rm for}\quad 0\leq x\leq\pi,\\ B_{1}^{\prime}x+B_{2}^{\prime}&{\rm for}\quad\pi\leq x<2\pi.\end{array}\right.$$

$$A_1^{\prime }$$

$$B_1^{\prime }$$

$$-{\nabla }^{2}w(\mathbf{x})=f(\mathbf{x}),$$

$$w(\mathbf{x})=\sin x\cos 2y+{\alpha }_{x}x+{\alpha }_{y}y-({\alpha }_x+{\alpha }_y)\pi .$$

$$-\nabla \cdot (\theta \nabla v+\chi \mathbf{\alpha })=(1-\chi )f,$$

$$\chi({\bm{x}})=\left\{\begin{array}[]{ll}0&{\rm if}\,\,\,{\bm{x}}\in\Omega_{f},\\ 1/4&{\rm if}\,\,\,{\bm{x}}=(\pi/2,\pi/2),(\pi/2,3\pi/2),\\ &\quad\quad\,\,\,\,\,\,\,(3\pi/2,\pi/2),(3\pi/2,3\pi/2),\\ 1/2&{\rm if}\,\,\,x=\pi/2,\,3\pi/2\,\,{\mathrm{for}}\,\,\pi/2<y<3\pi/2,\,\,{\mathrm{or}}\,\,y=\pi/2,\,3\pi/2\,\,{\mathrm{for}}\,\,\pi/2<x<3\pi/2\\ 1&{\rm otherwise}.\end{array}\right.$$

$$f(\mathbf{x})=f(r)=16\cos 4r+\frac{4\sin 4r}{r},$$

$$\frac{\partial w}{\partial r}=3\alpha \mathrm{at}r=\frac{\pi }{4},\mathrm{and}\frac{\partial w}{\partial r}=\alpha \mathrm{at}r=\frac{3}{4}\pi .$$

Full text

I Introduction

Numerical modeling of multiphysics problems in complicated geometries is still a challenging problem in computational fluid dynamics. For example, the numerical prediction of the exchanges of scalar quantities, e.g., heat and mass, between solid and fluid phases has numerous industrial applications. A non-exhaustive list includes radiators and heat exchangers, for example the fan cooler of CPUs, or crystal growth processes with important applications for semi-conductors and light emission devices. The heat or mass flux through interfaces between solid and fluid phases can be mathematically modeled, and leads typically to Neumann or Robin boundary conditions of the governing advection-diffusion equations.

Immersed boundary methods yield an attractive approach for solving partial differential equations in domains of complex shape. The underlying idea is to consider a simple computational domain for which effective numerical methods are available and the boundary conditions are imposed by so-called penalty or direct forcing terms; for a review we refer to Refs. Peskin ; Mittal ; reviewKai . For Dirichlet boundary conditions, there is an abundant literature available. Penalization techniques for Neumann and Robin boundary conditions are more recent and less developed, see e.g., Refs. KaKoAnSc ; KoNgSc ; Morales ; Vasilyev ; Ren2013 ; BCL2015 . Most immersed boundary methods result in low order, i.e., first or second order, approximation and computational stiffness.

In this paper, we focus on the volume penalization (VP) method. In the pioneering work by Angot et al. Angot1999 , the VP method was developed for imposing Dirichlet boundary conditions on the Navier–Stokes equations. The boundary conditions for velocity of viscous flow are given by no-slip conditions on the surface of the solid, such as walls and obstacles. The VP method models the solid as a porous medium whose permeability $\eta\,(>0)$ is sufficiently small. Mathematically, it was shown that the solution of the penalized Navier–Stokes equations converges towards the solution of the Navier–Stokes equations with no-slip boundary conditions, as $\eta\rightarrow 0$ Angot1999 ; Carbou . This VP method has been applied to various flows, e.g., confined hydrodynamic turbulence KaiMarie , confined magnetohydrodynamic turbulence Morales ; KaiNeffaa , fluid-structure interaction for moving obstacles Kolmenskiy and for flexible beams KoEnSc , and the aerodynamics of insect flight KoMoFaSc .

A VP method for general boundary conditions of Neumann and Robin types was proposed by generalization of the VP method using a weak formulation Ramiere1 , and applied in the context of finite element methods and finite volume methods Ramiere2 . Kadoch et al. KaKoAnSc extended this approach for pseudo-spectral discretizations using a strong formulation in the case of homogeneous Neumann boundary conditions.

The aim of the current work is to extend the VP method, developed in Ref. KaKoAnSc for homogeneous Neumann boundary conditions to inhomogeneous ones. Thus net scalar flux can be modeled keeping the continuity of its flux. The proposed method is analyzed in detail for a one-dimensional (1D) Poisson equation where exact solutions for both the penalized equation and the non-penalized equation are available. The penalization error can be hence computed analytically. For the numerical discretization, second order finite-differences and interpolation are used. The discretization error is determined and a guide for choosing the spatial discretization and $\eta$ in actual implementations is presented.

This paper is organized as follows. In Section II, we develop a flux-based VP representation of inhomogeneous Neumann boundary conditions. In Section III, the VP representation is extended to a two-dimensional (2D) penalized Poisson equation. In Section IV, we apply the VP representation to an advection-diffusion equation coupled with Navier–Stokes equations in an annular domain. Numerical results are given for 2D steady incompressible convection in an annulus subjected to heat flux through its inner wall. Conclusions are given in Section V.

II Poisson equation with inhomogeneous Neumann boundary conditions

We consider a 1D Poisson equation with inhomogeneous Neumann boundary conditions modeled with the VP method. In this case, the exact solution of the non-penalized and the penalized equations can be obtained and thus the penalization error can be assessed. The discretization error of the penalized equation using finite differences is also analyzed.

II.1 A flux-based volume penalization representation of inhomogeneous Neumann boundary conditions

We first consider a 1D Poisson equation,

[TABLE]

with a source term $f(x)$, which is here given by

[TABLE]

Here, $m=1,2,\cdots$, and $d_{x}=d/dx$. The fluid domain is $\Omega_{f}=\{x\ |\ 0<x<\pi\}$ and inhomogeneous Neumann boundary conditions,

[TABLE]

are imposed, where $\alpha$ is a real valued constant. The argument $x$ will be omitted, unless otherwise stated.

The compatibility condition, $-[d_{x}w]_{0}^{\pi}=\int_{0}^{\pi}f(x)dx\,(=0)$, is satisfied. Therefore, exact solutions of Eq. (1) with Eq. (2) exist satisfying the condition (3). For simplicity, we have chosen here the same value $\alpha$ on the left and right boundaries, $x=0$ and $x=\pi$. If $\alpha=0$, the boundary conditions (3) are homogeneous. In Section II.3 we will generalize this case to inhomogeneous Neumann boundary conditions with different values on the boundaries. The exact solutions of Eqs. (1)–(3) are given by

[TABLE]

where $A_{0}$ is an additive constant which reflects the non-uniqueness of the solution. To determine $A_{0}$, we set the condition $\int_{\Omega_{f}}w(x)dx=0$, which yields the value $A_{0}=-\pi\alpha/2$.

Generalizing the VP representation for homogeneous Neumann boundary conditions KaKoAnSc , we obtain the VP representation of Eqs. (1) and (3):

[TABLE]

where

[TABLE]

and $\eta\,(>0)$ is the penalization parameter. The mask function $\chi(x)$ determines the geometry of $\Omega_{f}$, and is given by

[TABLE]

Here, we have chosen to imbed the fluid domain $\Omega_{f}=\{x\ |\ 0<x<\pi\}$ into the larger domain $\Omega=\{x\ |\ 0\leq x<2\pi\}={\bar{\Omega}_{f}}\cup\Omega_{s}$ with the solid domain $\Omega_{s}=\{x\ |\ \pi<x<2\pi\}$, and ${\bar{\Omega}_{f}}$ denoting the closure of $\Omega_{f}$. Then we can impose $2\pi$-periodic boundary conditions. This VP representation (6) is based on the continuity of the flux through the non-penalized boundaries at $x=0\,(=2\pi)$ and $x=\pi$. The flux at the fluid-solid interfaces satisfies

[TABLE]

Now we solve Eq. (5) analytically, and obtain

[TABLE]

Using Eq. (11) and the $C^{0}$-continuity of the solutions at the interface of $\Omega_{f}$ and $\Omega_{s}$, i.e., $v(0^{+})=v(2\pi^{-})$ and $v(\pi^{+})=v(\pi^{-})$, we can determine three out of the four constants.

[TABLE]

Under the condition $\int_{\Omega_{f}}v(x)dx=0$, we get

[TABLE]

In the limit of $\eta\rightarrow 0$, the penalized solution $v(x)$, Eq. (12), converges towards the exact solution $w(x)$, Eq. (4), of the non-penalized problem. For even $m$ and $\alpha=0$, $v(x)$ is identical to $w(x)$. Fig. 1 shows $v(x)$ at $\eta=10^{-1}$ and $10^{-8}$ for $m=1$ together with $w(x)$. It is seen that for decreasing value of $\eta$, $v(x)$ perfectly matches $w(x)$. For the first derivative $d_{x}v$, we observe that no penalization boundary layer is present, which is in good agreement with the findings in Ref. KoNgSc . The second derivative $d_{x}^{2}v$ confirms the absence of the penalization boundary layer.

The penalization error $\left|v(x)-w(x)\right|$ ($=\left|(A_{1}-\alpha)x+A_{2}-A_{0}\right|$ in $\Omega_{f}$) is $O(\eta)$, because $A_{1}=\alpha+O(\eta)$, $A_{2}=A_{0}+O(\eta)$ and $A_{0}=-\pi\alpha/2$. We also find that $\left|d_{x}v-d_{x}w\right|$ is constant and we have again an error of $O(\eta)$. Fig. 2 shows that the difference $\left|v(x)-w(x)\right|$ is symmetric with respect to $x=\pi/2$ and takes maximum values at the boundaries of the fluid domain, i.e., at $x=0$ or $x=\pi$. Figs. 1 and 2 suggest that for the largest value of $\eta$, i.e., $\eta=10^{-1}$, the penalization approach does not well approximate the boundary condition of the non-penalized problem.

To study its regularity, we also analyze the Fourier coefficients ${\hat{v}}_{k}$ of the penalized solution $v(x)$, Eq. (12), where ${\hat{v}}_{k}$ is defined by

[TABLE]

The Fourier coefficients ${\hat{v}}_{k}$ for odd $m$ can be computed explicitly and are given by

[TABLE]

while the Fourier coefficients for even $m$ are given by

[TABLE]

For sufficiently large $|k|$ in the sense that $|k|\gg|m|$, Eqs. (23) and (29) show that $|{\hat{v}}_{k}|$ decays proportional to $k^{-2}$ for odd wave numbers $k$, whereas it decays as $O(k^{-3})$ only for even $k$ and odd $m$. Fig. 3 shows the decay of $|{\hat{v}}_{k}|$ for $\eta=10^{-2}$. There is no intermediate region at low $k$, as it is the case for the VP representation of homogeneous Neumann boundary conditions KoNgSc . The absence of the region is in contrast to what is observed for the VP representation of the Dirichlet boundary condition Romain , and is thus attributed to the absence of any boundary layer due to the penalization of Neumann boundary conditions.

II.2 Discretization error of the second order finite-difference scheme

We discretize the penalized equation, using the second order central finite-differences of the first derivative of $\phi(x)$ given by

[TABLE]

and second order interpolation of $\phi$,

[TABLE]

where $\phi=v,d_{x}v,\,\theta d_{x}v$, $x=x_{i},x_{i}\pm h/2$, $x_{i}=ih$ $(i=0,\cdots,N-1)$, $h\,(=2\pi/N)$ is the grid width, and $N$ is the number of grid points.

We then analyze the discretization and penalization errors. In Figs. 4 and 5, we plot the errors between the numerical solution of the penalized equation (5) and the exact solution (4) of the non-penalized problem, Eqs. (1)–(3), as a function of $h$. The errors are measured in $\ell^{\infty}$-, $\ell^{1}$- and $\ell^{2}$-norms, and are computed only in the fluid domain $\Omega_{f}$. The grid points $x_{0}$ and $x_{N/2}$ are on the boundary of $\Omega_{f}$.

For small $\alpha\,(=0.1)$, Fig. 4 shows that the error measured in each norm has a pronounced minimum, which corresponds to the optimal value of $h$, for a given penalization parameter $\eta$. By decreasing $h$ below this value, we see that the errors increase, and eventually saturate. This saturation is due to the penalization error. For smaller $\eta$, the errors decay approximately as $O(h^{2})$ with decreasing $h$. The minimum value and the optimal value of $h$ decrease, as $\eta$ becomes smaller. These behaviors are the same as those obtained by Ref. KoNgSc for the penalized problem of the homogeneous Neumann boundary condition, i.e., $\alpha=0$. It can be seen that for $h\approx 10^{-2}$, the error has its minimum for $\eta=10^{-5}$, for the three norms shown in Fig. 4.

However, there is no pronounced minimum for larger values $\alpha\,(=1)$. Fig. 5 shows that, for larger $\eta\,(=10^{-2})$, the errors for all considered norms are almost independent of $h$. For $\eta=10^{-3}$, the errors decay approximately as $O(h^{2})$ with decreasing $h$ and then saturate. The level of the saturation becomes smaller, as $\eta$ decreases. The $O(h^{2})$ convergence is again observed at sufficiently small $\eta$. For much larger $\alpha\,(=10)$, we observe the behavior similar to the case of $\alpha=1$ (figure omitted).

It is suggested in Section II.3 and Appendix A that this $O(h^{2})$ convergence is attributed to the cancellation of the errors of $O(h)$. The cancellation is due to the symmetry of Eq. (5) with respect to $x=\pi/2$.

The discrete eigenvalue problem of the penalized Laplace operator results in the same as the problem in the homogeneous case, computed in Ref. KoNgSc , that shows two different behaviors in the spectrum of the penalized operator. The eigenvalues of one part of the spectrum converge to the eigenvalues of the non-penalized Laplace operator. For the other part, corresponding to the small eigenvalues, these values depend on $\eta$ and vanish at $\eta\rightarrow 0$. Let $F$ be the discretized and penalized Laplace operator which is obtained by application of the second order finite-difference scheme, Eqs. (30) and (31), to $-d_{x}(\theta d_{x})$. The corresponding eigenvalue problem has the form $F{\bm{v}}={\bm{b}}$. The circulant matrix $F$ can be diagonalized applying Fourier series expansion, and we have ${\mathrm{ker}}F=1$. Hence it is not invertible. To remove this kernel, we use the constraint, a discretized version of $\int_{\Omega_{f}}v(x)dx=0$ using the midpoint rule. Under periodic boundary conditions, the eigenfunctions can be represented by the Fourier series. Thus, all of the non-zero eigenvalues are positive.

II.3 1D Poisson equation with different inhomogeneous Neumann boundary conditions

In this section, we study a VP representation for two different values of inhomogeneous Neumann boundary conditions. We consider again a 1D Poisson equation Eq. (1), $-d_{x}^{2}w(x)=f(x)$, with the source term given by

[TABLE]

for odd $m$ in the domain $\Omega_{f}=\{x\ |\ 0<x<\pi\}$ and imposing

[TABLE]

The compatibility condition is fulfilled, since we obtain $-[d_{x}w]_{0}^{\pi}=2m$ and $\int_{0}^{\pi}f(x)dx=2m$. Imposing that $\int_{\Omega_{f}}w(x)dx=0$, we find that the exact solution is given by

[TABLE]

The VP representation of Eqs. (1) and (33) reads

[TABLE]

in the extended domain $\Omega=\{x\ |\ 0\leq x<2\pi\}$ under $2\pi$-periodic boundary conditions, where the mask function $\chi(x)$ is given by Eq. (10). Here we set $\beta(x)=\alpha+m\cos(mx)$, which satisfies $\beta(0)=\alpha+m$ and $\beta(\pi)=\alpha-m$. The second term on the right-hand side of Eq. (35), $-\chi d_{x}\beta$, is added such that Eq. (35) reduces to $\eta d_{x}^{2}v=0$ in the solid domain $\Omega_{s}\,(=\Omega\backslash{\bar{\Omega}}_{f})$, and thus the penalized solution converges towards the non-penalized one (34) for $\eta\rightarrow 0$. The penalized representation (35) is again based on the continuity of the flux though the non-penalized boundary at $x=0\,(=2\pi),\pi$. The flux at the interface satisfies

[TABLE]

We obtain the penalized solution of Eq. (35);

[TABLE]

The coefficients can be determined in the same way for Eq. (12), and we obtain

[TABLE]

The penalized solution $v$, Eq. (39), results in the non-penalized solution $w$, Eq. (34), for $\eta\rightarrow 0$. We find $|v-w|=O(\eta)$ in the fluid domain $\Omega_{f}$, using Eq. (40). Again it is found that $\left|d_{x}v-d_{x}w\right|$ is constant and of $O(\eta)$. Equation (39) is independent of $\eta$ for $\alpha=0$. Fig. 6 presents the comparison of the penalized solution $v$ with the non-penalized solution $w$ for $\alpha=1$ and $m=1$ using different values of $\eta$. It is seen that they excellently agree with each other for $\eta=10^{-8}$. The difference $|v-w|$, which characterizes the penalization error at each position $x$, is shown in Fig. 7. For each $\eta$, it is symmetric with respect to $x=\pi/2$. The maximum value is located at the interface $x=0$ or $x=\pi$, and decreases as $\eta$ becomes smaller.

The numerical solutions of the non-penalized problem, Eqs. (1), (32) and (33), and the penalized equation (35) are obtained by using the second order finite-differences and interpolation in Eqs. (30) and (31). The errors measured by the different norms are plotted in Fig. 8. The errors decay with decreasing $\eta$ and the grid width $h(=2\pi/N)$. For smaller $\eta\,(=10^{-3}$, $10^{-8})$, we observe that the errors exhibit first order convergence in terms of $h$. This is in contrast to the case of Section II.1 where identical values of the inhomogeneous boundary conditions are considered on the left and the right boundary. This $O(h)$ convergence is thus due to the different values of $\beta(x)$ at the interfaces $x=0,\pi$, i.e., $\beta(0)\neq\beta(\pi)$ (see Appendix A).

The developed VP method for inhomogeneous Neumann boundary conditions can be also applied to time-dependent problems. A convergence issue in terms of $\eta$ and $h$ is discussed in Appendix B for a 1D penalized thermal diffusion equation.

III Application to a 2D Poisson equation

Now, we extend the VP representation for inhomogeneous Neumann boundary conditions to 2D domains. Let us consider a Poisson equation in the 2D fluid domain $\Omega_{f}=\{(x,y)\ |\ \pi/2<x<3\pi/2\ {\rm and}\ \pi/2<y<3\pi/2\}$;

[TABLE]

imposing inhomogeneous Neumann boundary conditions ${\bm{n}}{\bm{\cdot}}\nabla w={\bm{n}}{\bm{\cdot}}{\bm{\alpha}}$ at the boundary of $\Omega$, i.e., $x=\pi/2,3\pi/2$ for $\pi/2\leq y\leq 3\pi/2$, and $y=\pi/2,3\pi/2$ for $\pi/2\leq x\leq 3\pi/2$, where ${\bm{x}}=(x,y)$, $\nabla=(\partial_{x},\partial_{y})=(\partial/\partial x,\partial/\partial y)$, ${\bm{\alpha}}=(\alpha_{x},\alpha_{y})$, $\alpha_{x}$ and $\alpha_{y}$ are constant, and ${\bm{n}}$ is the outward pointing unit normal vector of $\Omega_{f}$. We set the source function as $f({\bm{x}})=5\sin{x}\cos{2y}$, and verified the compatibility condition of Eq. (42). Then we obtain an exact solution;

[TABLE]

We imposed the condition that $\int\!\!\int_{\Omega_{f}}w({\bm{x}})d{\bm{x}}=0$, to determine the solution uniquely.

A VP representation of Eq. (42) reads

[TABLE]

which may be solved in the double periodic domain $\Omega=\{(x,y)\ |\ 0\leq x<2\pi,0\leq y<2\pi\}$. The mask function $\chi({\bm{x}})$ is defined by

[TABLE]

Imposing $\int\!\!\int_{\Omega_{f}}v({\bm{x}})d{\bm{x}}=0$, we select the penalized solution of Eq. (44) that converges to the exact non-penalized solution (43) for $\eta\rightarrow 0$. Fig. 9 illustrates the numerical solution of Eq. (44) at $\eta=10^{-8}$ for ${\bm{\alpha}}=(2,1)$, which is obtained by the second order finite-differences and interpolation in Eqs. (30) and (31), using $256$ grid points in each Cartesian direction. We observe that the penalized solution $v({\bm{x}})$ is in excellent agreement with the exact solution (43) in $\Omega_{f}$.

Fig. 10 shows the dependence of the $\ell^{2}$ error between the numerical solution $v$ and the exact solution $w$, Eq. (43), in $\Omega_{f}$ on the grid width $h$ for $\eta=10^{-2}$ and $10^{-8}$. For $\eta=10^{-2}$, the $\ell^{2}$ error has the minimum in terms of $h$. The error becomes almost constant for smaller values of $h\,(\lesssim 0.1)$. In contrast, for much smaller $\eta\,(=10^{-8})$, it is seen that the error monotonically decays with decreasing $h$, and behaves approximately as $O(h^{2})$. This $O(h^{2})$ convergence is consistent with what was found in Fig. 5 for the 1D penalized Poisson equation, as well as for the 2D penalized Poisson equation with homogeneous Neumann boundary conditions in a square domain KoNgSc . The errors measured by the $\ell^{1}$- and $\ell^{\infty}$-norms are omitted, because the errors show the same behaviors as those of the $\ell^{2}$ error.

Next, let us consider a 2D Poisson equation (42) with inhomogeneous Neumann boundary conditions in an annular fluid domain $\Omega_{f}=\{(x,y)\ |\ \pi/4<r<3\pi/4\}$, where $r=\sqrt{(x-\pi)^{2}+(y-\pi)^{2}}$. We set the source term $f({\bm{x}})$ as

[TABLE]

and impose inhomogeneous Neumann boundary conditions given by

[TABLE]

The compatibility condition $\int_{\pi/4}^{3\pi/4}rf(r)dr=0$ is fulfilled. Noting $\nabla^{2}w=(1/r)d_{r}(rd_{r})$ and $w({\bm{x}})=w(r)$, we obtain the exact solution;

[TABLE]

where

[TABLE]

imposing the condition that $\int_{\pi/4}^{3\pi/4}rw(r)dr=0$.

The penalized equation of this problem reads

[TABLE]

where $\Omega_{f}$ is imbedded into the $2\pi$ double periodic domain $\Omega=\{(x,y)|0\leq x<2\pi,\,0\leq y<2\pi\}$, ${\bm{\beta}}({\bm{x}})=g({\bm{x}}){\bm{e}}_{r}$, ${\bm{e}}_{r}=(x-\pi,y-\pi)/r$ and

[TABLE]

Note that ${\bm{\beta}}({\bm{x}})=3\alpha{\bm{e}}_{r}$ at $r=\pi/4$, and ${\bm{\beta}}({\bm{x}})=\alpha{\bm{e}}_{r}$ at $r=3\pi/4$, which is consistent with Eq. (47). The mask function $\chi({\bm{x}})$ is given by

[TABLE]

Fig. 11 shows 1D cuts of the numerical solution $v$, which are compared with those of the exact solution $w$, Eq. (48), of the non-penalized problem for $\alpha=1$ and $\eta=10^{-8}$. The numerical solutions were computed in Cartesian coordinates using $256$ grid points in each direction. It is seen that $v$ is in good agreement with $w$. Fig. 12 shows the dependence of the $\ell^{\infty}$ error between the numerical solution $v$ and the exact one $w$ as a function of the grid width $h$ in $\Omega_{f}$. The error decays approximately as $O(h)$. This $O(h)$ convergence is attributed to the following two errors. One is the absence of grid points at the circular interfaces, i.e. for $r=\pi/4$ and $3\pi/4$ (see Ref. KoNgSc ). The other is the use of the different boundary values for the inhomogeneous Neumann boundary conditions, as already discussed in Section II.3 and Appendix A.

IV Application to free convection in two dimensions

We now apply the VP method developed in this paper to an advection-diffusion equation with an inhomogeneous Neumann boundary condition and a Dirichlet boundary condition. The equation is coupled with the incompressible Navier–Stokes equations. Numerical simulations were performed for 2D steady thermal flows in a concentric annulus that is heated from its inner cylinder. The VP method is examined by comparison with numerical results of Refs. Ren2013 ; Yoo2003 .

IV.1 Volume penalized governing equations

We consider steady free convection of an incompressible fluid in a concentric annulus. Its configuration is illustrated in Fig. 13, where $r_{i}$ and $r_{o}$ are respectively the inner and outer cylinder radius, $L=r_{o}-r_{i}$, $T_{o}$ is constant temperature of the outer cylinder, and $Q$ is a constant heat flux from the inner cylinder. In this study, we set $r_{i}/L=1$ and $r_{o}/L=2$, which was used in Refs. Ren2013 ; Yoo2003 . Let ${\bm{u}}$, $T$ and ${\bm{x}}$ be respectively the velocity, the temperature and the position. The dimensionless velocity ${\bm{u}}^{*}$ is given by ${\bm{u}}^{*}={\bm{u}}L/\kappa$, the dimensionless temperature $\Phi^{*}$ is denoted by $\Phi^{*}=(T-T_{o})k/(QL)$, and ${\bm{x}}^{*}={\bm{x}}/{L}$, where $\kappa$ is the thermal diffusivity defined by $\kappa=k/(\rho C_{p})$; $k$ is the thermal conductivity, $\rho$ is the density, and $C_{p}$ is the heat capacity at constant pressure. We will drop the superscript ∗ hereafter, unless otherwise stated.

The non-penalized governing equations under the Boussinesq approximation are given by

[TABLE]

where $p$ is the dimensionless pressure, $\partial_{t}=\partial/\partial t$, ${\bm{e}}_{y}=(0,1)$, ${\mathrm{Pr}}=\nu/\kappa$ is the Prandtl number, ${\mathrm{Ra}}=g\gamma QL^{4}/(k\kappa\nu)$ is the Rayleigh number, and $\gamma$ is the thermal expansion coefficient. The fluid domain $\Omega_{f}$ is given by $\Omega_{f}=\{(x,y)\ |\ 1<r<2\}$, where $r=\sqrt{x^{2}+y^{2}}$. No-slip boundary conditions for velocity, ${\bm{u}}={\bm{0}}$, are imposed at $r=1$ and $r=2$. The boundary conditions of $\Phi$ are given by

[TABLE]

where ${\bm{n}}=(x,y)$ is the outward normal vector of the inner cylinder surface.

The no-slip conditions and Eq. (62), which are Dirichlet boundary conditions, are modeled by the classical VP method Angot1999 . In contrast, Eq. (61), which is an inhomogeneous Neumann boundary condition, can be modeled by the flux-based VP method we have proposed. We here consider flow motions in the extended square domain $\Omega=\{(x,y)|-2.56\leq x<2.56,-2.56\leq y<2.56\}$ under the periodic conditions at $x=\pm 2.56$ and $y=\pm 2.56$. The resulting penalized equations read

[TABLE]

with Eq. (59), where $\eta_{d}\,(>0)$ and $\eta_{n}\,(>0)$ are the penalization parameters, ${\bm{\beta}}({\bm{x}})={\bm{x}}/r$ with ${\bm{\beta}}({\bm{0}})={\bm{0}}$, and the mask function $\chi({\bm{x}})$ is defined by $\chi({\bm{x}})=\chi_{1}({\bm{x}})+\chi_{2}({\bm{x}})$ in which $\chi_{1}$ and $\chi_{2}$ are given by

[TABLE]

We can select another expression for ${\bm{\beta}}$ as far as ${\bm{n}}{\bm{\cdot}}{\bm{\beta}}=1$ at $r=1$. Note that $\nabla{\bm{\cdot}}(\chi_{2}{\bm{\beta}})-\chi_{2}(\nabla{\bm{\cdot}}{\bm{\beta}})=0$ for $r<1$. In Appendix C, we discuss the VP representation of mixed Dirichlet–Neumann boundary conditions for a 1D Poisson equation.

We numerically solve Eqs. (63) and (64) with (59), using the marker-and-cell method with a staggered grid at $256\times 256$ grid points. Second order central finite-differences and interpolation in Eqs. (30) and (31) are employed for the spatial discretization, and the explicit Euler method is used as the time marching. The time increment $\Delta t$ is $10^{-6}$. The computations were carried out at ${\mathrm{Pr}}=0.7$ for two ${\mathrm{Ra}}$ numbers, ${\mathrm{Ra}}=5700$ and ${\mathrm{Ra}}=5\times 10^{4}$. The penalization parameters are taken as $\eta_{d}=\eta_{n}=5\times 10^{-6}$. The Poisson equation for the pressure $p$ is solved by the use of a successive over relaxation method. The convergence of $p$ is here achieved, if the relative residual error of $p$ measured by the $\ell^{1}$-norm becomes less than $10^{-6}$. The numerical solutions are regarded as steady ones, if the relative time increments of $u$, $v$ and $\Phi$ normalized by $\Delta t$ respectively become less than $10^{-6}$, where the increments are measured by the $\ell^{1}$-norm, and ${\bm{u}}=(u,v)$.

We also performed numerical simulations at ${\mathrm{Ra}}=5\times 10^{4}$ and ${\mathrm{Pr}}=0.7$, using different number of grid points, penalization parameters $\eta_{n}=\eta_{d}$, and time increment $\Delta t$. The simulations confirmed that the numerical results for $128^{2}$ with $\eta_{d}=\eta_{n}=5\times 10^{-6}$ and $\Delta t=10^{-6}$ excellently agree with those at $256^{2}$ using the same values of $\eta_{d},\eta_{n}$, and $\Delta t$. Such an excellent agreement of the numerical results is also observed, when doubling either $\Delta t$ or $\eta_{d}(=\eta_{n})$. This convergence study justifies the choice of the parameters in the above computations. Note that for stability reasons due to the explicit time discretization, we have chosen the values of $\Delta t$ so as to satisfy $\Delta t<\eta_{d}$, $\Delta t<\eta_{n}$ and $\Delta t<1/{\mathrm{Ra}}$.

IV.2 Numerical results

Our numerical results are compared with results in Refs. Ren2013 ; Yoo2003 in which also second order finite-differences are used. Yoo Yoo2003 computed the flows in the polar coordinate system, while Ren et al. Ren2013 used a Cartesian coordinate system combined with an immersed boundary method for Neumann boundary conditions which they developed. Ren et al. Ren2013 reported that their results are in accordance with those of Ref. Yoo2003 .

Fig. 14 shows the streamlines and the isotherms of $\Phi$ at ${\mathrm{R}}=5700$. The streamlines and isotherms are symmetric with respect to $x=0$. In Fig. 14 (b), we observe that the isotherms at $0.1,\ 0.2,\ 0.3,\ 0.4,\ 0.5$ and $0.6$ excellently agree with those presented in Fig. 7 (b) of Ref. Ren2013 . We omitted the other isotherms of our results for brevity, since they also show excellent agreement with those of Ref. Ren2013 .

The spatial distributions of $\Phi$ on the inner cylinder at ${\rm Ra}=5700$ and $5\times 10^{4}$ are shown in Fig. 15, and are compared with the results in Refs. Ren2013 ; Yoo2003 . The angle $\Theta$ is measured counterclockwise from the top of the inner cylinder such that $\Theta=0^{\circ}$ at $(x,y)=(0,1)$ and $\Theta=180^{\circ}$ at $(x,y)=(0,-1)$. It can be seen that our results are in good agreement with those of Ref. Yoo2003 at both ${\rm Ra}$ numbers and the results of Ref. Ren2013 at ${\rm Ra}=5700$. At ${\rm Ra}=5\times 10^{4}$, our result shows even better agreement with Ref. Yoo2003 than with Ref. Ren2013 . The reciprocal of $\Phi$ at $r=1$ gives an average Nusselt number that characterizes the average heat transfer ratio through the inner cylinder surface Ren2013 .

Note added in proof: Guo et al. new have developed an immersed boundary method for inhomogeneous Neumman boundary conditions in the frameworks of a finite volume method. They used the same test case for validation, and showed that their numerical results likewise agree with those in Ref. Yoo2003 .

V Conclusions

We have developed a VP representation for inhomogeneous Neumann boundary conditions, generalizing the technique proposed by Kadoch et al. KaKoAnSc for homogeneous Neumann boundary conditions. It is based on the flux continuity between the fluid and solid domains. The method was first applied to a 1D Poisson equation with inhomogeneous Neumann boundary conditions. The analytical solution of the penalized Poisson equation converges to the exact solution of the non-penalized 1D Poisson equation for the penalization parameter $\eta\rightarrow 0$. The numerical solutions of the penalized equation were obtained by the use of second order central finite-differences and interpolation. In the case that the same values of the inhomogeneous Neumann boundary conditions are imposed at two interfaces between fluid and solid domains, the errors between the exact solution of the non-penalized equation and the numerical solutions of the discretized equation exhibit second order convergence in terms of the grid width $h$ for sufficiently small $\eta$. It was suggested that the $O(h^{2})$ convergence is due to the cancellation of the discretized errors of the interface between the solid and fluid domains. The eigenvalues of the discretized Laplace operator are identical for homogeneous and inhomogeneous boundary conditions and were studied in Ref. KoNgSc . We also showed that the VP representation needs a source term in the solid domain for the Poisson equation imposing two different values of inhomogeneous Neumann boundary conditions. This term allows for the convergence of the penalized solution to the non-penalized one in the limit of $\eta\rightarrow 0$. The numerical solutions of the penalized equation exhibit $O(h)$ convergence. Our results suggest that the value of $\eta$ in the range of $10^{-5}\lesssim\eta\lesssim 10^{-3}$ is optimal for $h\approx 10^{-2}$.

The developed VP method was applied to the advection-diffusion equation coupled with the Navier–Stokes equations. We simulated steady free convection of an incompressible fluid in a concentric annulus heated through the inner cylinder surface. It was found that the VP method well preserves the temperature distributions on the inner cylinder obtained in the literature Ren2013 ; Yoo2003 .

This work motivates us to develop a VP method further for Robin boundary conditions, and apply the developed method to multiphysics problems in complex geometries, e.g., Stefan problems in crystal growth and to magnetohydrodynamic problems in magnetically confined fusion.

Acknowledgments

The authors thank Dr. Dmitry Kolomenskiy for fruitful discussion about the convergence issue in the penalized Poisson equation with homogeneous Neumann boundary condition. This work was partially supported by JSPS KAKENHI Grant Numbers (S)16H06339 and (C)17K05139. K.S. acknowledges support by the French Federation for Magnetic Fusion Studies (FR-FCM) and of the Eurofusion consortium from the Euratom research and training programme 2014-2018 and 2019-2020 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

Appendix A On different values of inhomogeneous Neumann conditions

To get deeper insight into the discretization error of the VP representation for inhomogeneous Neumann boundary conditions, we consider the Poisson equation (1) with the source term $f(x)$ given by

[TABLE]

in the fluid domain $\Omega_{f}=\{x\ |\ 0\leq x<\pi\}$ imposing the inhomogeneous Neumann boundary conditions

[TABLE]

where $\epsilon$ and $\alpha$ are constant. Under the condition $\int_{\Omega_{f}}w(x)dx=0$, the exact solution is given by

[TABLE]

Applying the VP representation to Eqs. (66) and (67), we obtain

[TABLE]

where $\beta(x)=\alpha+\epsilon\cos(x)$, $\Omega_{f}$ is extended to the domain $\Omega=\{x\ |\ 0\leq x<2\pi\}$, and then $2\pi$-periodic boundary conditions are imposed. The functions $\theta(x)$ and $\chi(x)$ are respectively given by Eqs. (6) and (10). The analytical solution of Eq. (69) is obtained in the same way as in Section II.1, and is omitted here for brevity.

Fig. 16 shows that the $\ell^{2}$ error between the numerical solutions $v$ of Eq. (69) and the exact one (68) for $\eta=10^{-8}$ using different $\epsilon$. We see the $O(h^{2})$ convergence of the error at $\epsilon=0$ which degrades to an $O(h)$ convergence for $\epsilon=1$. It is observed at $\epsilon\,=0.01$ that this error decays approximately as $O(h^{2})$ for large $h\,(\gtrsim 5\times 10^{-2})$, and they show the $O(h)$ decay for smaller $h\,(\lesssim 10^{-2})$. The errors measured in the $\ell^{1}$- and $\ell^{\infty}$-norms show the similar behavior, and their figures are omitted for brevity.

Appendix B Penalized 1D thermal diffusion equation

We here consider a time-dependent problem corresponding to a 1D heat equation in the domain $\Omega_{f}=\{x\ |\ -\pi<x<\pi\}$:

[TABLE]

with the initial condition given by

[TABLE]

and the inhomogeneous time-dependent boundary conditions given by

[TABLE]

where $\partial_{t}=\partial/\partial t$. The exact solution then reads

[TABLE]

The VP representation of Eqs. (70) and (72) results in

[TABLE]

where $\alpha=-e^{-t}$, and $\chi(x)$ is a mask function defined by

[TABLE]

Here, we employ a periodic domain $\Omega=\{x\ |\ -\pi-0.2\leq x<\pi+0.2\}$. Based on Eqs. (72) and (73), we use the initial condition $\varphi_{p}(x,0)=(1-\chi)\sin{x}$.

Equation (74) was numerically solved, using second order central finite-differences and interpolation with $N=512$ grid points. The Crank-Nicolson method is employed for time integration with time increment $\Delta t=10^{-5}$. The penalization parameter $\eta$ is set to be either $10^{-2}$ or $10^{-8}$. We obtain the numerical solution $\varphi_{p}$, which is compared with the exact solution $\varphi(x,t)$ of Eq. (73). Fig. 17 shows $\varphi_{p}$ and $\varphi$ at $t=1$. In $\Omega_{f}$, $\varphi_{p}$ excellently agrees with $\varphi$. To quantify the precision, we present the $\ell^{\infty}$ and $\ell^{2}$ errors between $\varphi_{p}$ and $\varphi$ at $t=1$ in Fig. 18. In Fig. 18 (a), it can be seen that the errors decay, as the grid width $h$ decreases. Here, in this section, $h=(2\pi+0.4)/N$. For smaller $h$, we observe a $O(h)$ decay, since no grid point is on the boundary of $\Omega_{f}$. Fig. 18 (b) depicts the $\eta$-dependence of the errors. The errors decay approximately as $O(\eta^{1/2})$ for larger $\eta\,(\gtrsim 10^{-4})$ and converge for smaller $\eta\,(\lesssim 10^{-5})$. The decay of $O(\eta^{1/2})$ has been proved in theorem 2.1 of Ref. KaKoAnSc .

Appendix C Mixed Dirichlet–Neumann boundary conditions for the 1D Poisson equation

We consider here the 1D Poisson equation (1) with the source term given by

[TABLE]

in the domain $\Omega_{f}=\{x\ |\ 0<x<\pi\}$ imposing mixed Dirichlet and inhomogeneous Neumann boundary conditions, i.e.,

[TABLE]

with $\alpha$ being constant. The exact solution of Eq. (79) reads

[TABLE]

Applying the VP representation to Eqs. (79) and (80), we obtain

[TABLE]

in which $\beta(x)=\alpha\left(1+\sin x\right)$, $\theta=\left(1-\chi_{n}\right)+\eta_{n}\chi_{n}$, and $\chi(x)=\chi_{d}(x)+\chi_{n}(x)$ where

[TABLE]

and

[TABLE]

The domain $\Omega_{f}$ is extended to $\Omega=\{x\ |\ 0\leq x<2\pi\}$, and $2\pi$-periodic boundary conditions are imposed.

The analytical solution of Eq. (82) is uniquely obtained not only by using $v(\pi^{-})=v(\pi^{+})$, $v(3\pi/2^{-})=v(3\pi/2^{+})$, $v(0^{+})=v(2\pi^{-})$, and the continuity of the flux, $d_{x}v=\eta d_{x}v+\beta$, at $x=\pi$ and $3\pi/2$, but also by imposing $d_{x}v(0^{+})=d_{x}v(2\pi^{-})$. For brevity we omit the form of the analytical solution, because it is a rather complicated expression. The analytical solution converges to the exact solution (81) for $\eta_{d}\rightarrow 0$ and $\eta_{n}\rightarrow 0$.

Fig. 19 shows that the errors between the numerical solutions of Eq. (82) and the exact solution $w$, Eq. (81), for $\alpha=1$ using different $\eta_{d}=\eta_{n}=\eta$. For $\eta=10^{-3}$, the errors decay approximately as $O(h)$ with decreasing $h$, and then saturate due to the penalization error. The level of the saturation becomes smaller, as $\eta$ decreases. It can be seen that for sufficiently small $\eta$ the errors decay like $O(h)$.