CFD and finite volume method: Dirichlet boundary conditions for the Euler equations

Please point me to an answer if one already exists, but after some searching, I still can't find the answer to what seems like a very simple question. There are plenty of references out there for solving the Euler equations with FVM. I am using the Toro (1997) book, Riemann Solvers and Numerical Methods for Fluid Dynamics, as a reference on Riemann solvers and the Godunov method.

I am attempting to solve the 2-dimensional Euler equations on a fixed, uniform grid of cells using the Godunov method and HLLC flux. Having implemented one of the Sod shock tube tests in the book, I feel confident that I haven't made a blatant error in my Riemann solver. The book does not address Dirichlet boundary conditions, however, which is exactly the boundary condition I'd like to implement.

Equations

I have the 2-dimensional Euler equations with a source term:

$$\frac{\partial}{\partial t} \begin{pmatrix} \rho \\ \rho \mathbf{u} \\ \rho E \end{pmatrix} + \nabla \cdot \begin{pmatrix} \rho \mathbf{u} \\ \rho \mathbf{u} \otimes \mathbf{u} + p \mathbf{I} \\ \left(\rho E + p\right) \mathbf{u} \end{pmatrix} = \mathbf{S}$$

where I have used bold notation for vectors, $$\rho$$ is the fluid density, $$\mathbf{u} = \left(u, v\right)$$ is its velocity, $$p$$ is the pressure, and $$\rho E = \frac{p}{\gamma - 1} + \frac{1}{2} \rho \left(u^2 + v^2\right)$$, where $$\gamma$$ is the (constant) ratio of specific heats.

Code

I tried to write my finite volume code in a general way so that I could easily test it and apply it to different kinds of problems. This means that each subproblem is expected to provide a function with a specific interface for each of the north, south, east, and west boundaries. To prescribe the Dirichlet boundary condition I want, there is a layer of ghost cells surrounding the mesh where I prescribe the fixed value $$\mathbf{g}_D$$. Then, I compute the flux as I normally would, specifying the fixed state in the appropriate slot for the Riemann solver. The resulting flux is then returned to the driver function, which just accumulates all the flux, weighted by the inverse grid spacing, into a big vector of derivatives.

Problem

The prescribed value on the ghost cells is easy to preserve by messing with how the source term is added, etc. Assuming the initial conditions satisfy $$\mathbf{g}_D$$, I can just set the total flux and source term to zero on the boundary so that the derivative there is zero. This made sense when I did it, but I'm no longer convinced it's correct, since I now get rather serious spikes in velocity at the boundaries of the domain where the value was forced to be constant (image below). The velocity spikes eventually work their way into every other field.

Question

Am I on the right track by using ghost cells with a prescribed value? Clearly, the exact method I've described doesn't work as expected. Is there a better way to accomplish this?

• What is the boundary supposed to be ? Solid wall boundary ? – cfdlab Nov 25 '19 at 2:33
• @cpraveen All fields are supposed to be fixed to zero for this particular problem, though there isn't a physical wall there. – emprice Nov 25 '19 at 13:28
• Well all fields fixed to 0 isn't a physical wall boundary... I tend to handle walls with a no penetration B.C. – EMP Nov 25 '19 at 21:43
• Have you taken a look at another solver such as clawpack? We tend to use two ghost cells and simply set the value in these ghost cells in the most simplistic case. This can lead to issues though so if you really want to have a value maintained at the boundary then you may need to solve a Riemann problem to guarantee the solution you want at the boundary. That being said it looks like you may have a bug in your code as I would expect your approach to work better than it is. – Kyle Mandli Nov 26 '19 at 1:54