[fadoobaba]

Recently, I have started developing my own CFD solver (hobby). I have had excellent success with in-compressible solver. (https://fluiddynamicscomputer.blogsp...-dynamics.html + https://fluiddynamicscomputer.blogsp...ierstokes.html + https://youtu.be/8RLIAkRKH64 + https://youtube.com/shorts/pKK5L6msmUI?feature=share etc.)

More recently, I have started reading about compressible equations and solvers. I am using continuity, convection diffusion (temperature), ideal gas (pressure), x and y momentum equations in time loop. I use MATLAB to test and then convert to C++ (open-source version) for publishing (speed difference is amazing, but this isn't the point of this thread ).

u(i, j) = (rho(i, j) .* u(i, j)...
- dt / (2 * h) * (rho(i+1, j) .* u(i+1, j) .* u(i+1, j) - rho(i-1, j) .* u(i-1, j) .* u(i+1, j)) ...
- dt / (2 * h) * (rho(i, j+1) .* u(i, j+1) .* v(i, j+1) - rho(i, j-1) .* u(i, j-1) .* v(i, j-1))...
- dt / (2 * h) * (p(i+1, j) - p(i-1, j)) ...
+ nu * (dt / h^2 * (u(i+1, j) - 2 * u(i, j) + u(i-1, j)) + dt / h^2 * (u(i, j+1) - 2 * u(i, j) + u(i, j-1)))) ./ rho(i, j); % x-momentum

This is supposed to be discretized x momentum in conservative form. Wonder what am I doing wrong? Tried different schemes but solution diverges. I get oscillations of velocity (horizontal lines in the domain). Doesn't matter the time step or mesh size.

[andy_]

Originally Posted by fadoobaba

Wonder what am I doing wrong?

You haven't introduced a mechanism to control pressure-velocity decoupling at a guess. That is, there is nothing in your scheme to oppose the oscillation you see from growing to any size. There are various schemes to control it though almost all end up doing pretty much the same thing and smoothing the pressure field albeit by different amounts.

[fadoobaba]

Originally Posted by andy_

You haven't introduced a mechanism to control pressure-velocity decoupling at a guess.

Yes, no I dont think I have done pressure velocity de coupling. How to this? I will read the literature for sure, but can you guide where to start?

Currently code inside time loop has 5 equations:

Continuity for density, convection diffusion for temperature, ideal gas for pressure, u and v momentum for velocities.

[andy_]

Originally Posted by fadoobaba

Yes, no I dont think I have done pressure velocity de coupling. How to this? I will read the literature for sure, but can you guide where to start?

Currently code inside time loop has 5 equations:

Continuity for density, convection diffusion for temperature, ideal gas for pressure, u and v momentum for velocities.

Every reasonable CFD book will discuss the range of issues that prevents a straightforward discretization of the governing transport equations for mass, momentum and energy from working. There are a range of approaches for handling the issues related to pressure-velocity decoupling, nonlineariy, compressiblity, turbulence closure and much more. Most text books will tend to focus on approaches relevant to a particular kind of flow rather than all flows. I would suggest finding a text book or course notes that is aligned with your area of interest and reading it for a few days while resisting the urge to do stuff on the computer. The links section of the website has some online sources if you don't have access to a library with a selection of CFD text books.

[FMDenaro]

I would add that a conservative FV scheme is not implemented this way! You have to compute the fluxes.

[fadoobaba]

Originally Posted by FMDenaro

I would add that a conservative FV scheme is not implemented this way! You have to compute the fluxes.

Thanks! I have read Hoffman CFD till chapter 8. I just used that knowledge to implement this. I will find another book

[Eicosane]

Originally Posted by fadoobaba

This is supposed to be discretized x momentum in conservative form. Wonder what am I doing wrong? Tried different schemes but solution diverges. I get oscillations of velocity (horizontal lines in the domain). Doesn't matter the time step or mesh size.

If I were you, in first place i would simplify solved system to Euler equations, i.e. exclude all diffusive terms from momentum and energy equations, and try to achieve stability and convergence for it. When i wrote my compressible solver, i used Toro's "Riemann Solvers and Numerical Methods for Fluid Dynamics". And only after that return to Navier-Stokes.

Second, as i see, your scheme is second order in space and first order in time. From my experience, this itself can lead to oscillations. Your scheme must be whether first or second order in space and time simultaneously. At least for convective term. For diffusion, i do not know, is difference in approximation order so important.

And third, looking at your equation, i have a question. Is this density: "./ rho(i, j)", at the end of equation the same as all other previous? If you solve compressible equations, it should be density from new time layer.

[fadoobaba]

Thanks for the book recommendation

In the in compressible flow part, I tried central and upwind for convective. To my surprise as well, central was more stable so I used that here. Same reason for viscous terms

Tried dividing by rho(i, j) here and rho_old(i, j) elsewhere and it showed no major difference. When I wrote in-compressible code, I used p_old (for pressure poisson), u_old, v_old etc. for velocities but then I deleted these _old matrices and replaced them with current time step one's. Code now runs 5 - 10% faster.

Here is compressible x-momentum I use:

u(i, j) = u(i, j) - u(i, j) * dt / (2 * h) .* (u(i+1, j) - u(i-1, j)) - v(i, j) * dt / (2 * h) .* (u(i, j+1) - u(i, j-1)) - dt / (2 * rho * h) * (p(i+1, j) - p(i-1, j)) + nu * (dt / h^2 * (u(i+1, j) - 2 * u(i, j) + u(i-1, j)) + dt / h^2 * (u(i, j+1) - 2 * u(i, j) + u(i, j-1))); % x-momentum

all u 's

[FMDenaro]

Originally Posted by fadoobaba

Thanks for the book recommendation

In the in compressible flow part, I tried central and upwind for convective. To my surprise as well, central was more stable so I used that here. Same reason for viscous terms

Tried dividing by rho(i, j) here and rho_old(i, j) elsewhere and it showed no major difference. When I wrote in-compressible code, I used p_old (for pressure poisson), u_old, v_old etc. for velocities but then I deleted these _old matrices and replaced them with current time step one's. Code now runs 5 - 10% faster.

Here is compressible x-momentum I use:

u(i, j) = u(i, j) - u(i, j) * dt / (2 * h) .* (u(i+1, j) - u(i-1, j)) - v(i, j) * dt / (2 * h) .* (u(i, j+1) - u(i, j-1)) - dt / (2 * rho * h) * (p(i+1, j) - p(i-1, j)) + nu * (dt / h^2 * (u(i+1, j) - 2 * u(i, j) + u(i-1, j)) + dt / h^2 * (u(i, j+1) - 2 * u(i, j) + u(i, j-1))); % x-momentum

all u 's

You are using a FD discretization on the quasi-linear form, I strongly discourage to use this form.

Furthermore:
- In the RHS, You are working on the same matrix u(i,j) of the update, this is wrong, you have to use two distinct arrays. What you are doing is somehow a Seidel step in time.

- in compressible flow, the diffusive term is NOT the laplacian.


I suggest to stop coding and studying the theory

- rho must be the value at (i,j)

[fadoobaba]

Correction from previous reply:

Here is in-compressible x-momentum I use:

u(i, j) = u(i, j) - u(i, j) * dt / (2 * h) .* (u(i+1, j) - u(i-1, j)) - v(i, j) * dt / (2 * h) .* (u(i, j+1) - u(i, j-1)) - dt / (2 * rho * h) * (p(i+1, j) - p(i-1, j)) + nu * (dt / h^2 * (u(i+1, j) - 2 * u(i, j) + u(i-1, j)) + dt / h^2 * (u(i, j+1) - 2 * u(i, j) + u(i, j-1))); % x-momentum

I am learning finite volume method from Malalasekera CFD book. This blog code [https://fluiddynamicscomputer.blogsp...lab-code.html] is just for my digital CV.

Indeed the current code is finite difference as I have studied Hoffman CFD till chapter 8 completely. I will read about compressible CFD and then resume coding.

[fadoobaba]

Originally Posted by FMDenaro

Furthermore:
- In the RHS, You are working on the same matrix u(i,j) of the update, this is wrong, you have to use two distinct arrays. What you are doing is somehow a Seidel step in time.

2 distinct arrays for each PDE makes code 5-10% slower with negligible difference in accuracy, hence I stopped using 2 arrays per PDE.

[FMDenaro]

Originally Posted by fadoobaba

2 distinct arrays for each PDE makes code 5-10% slower with negligible difference in accuracy, hence I stopped using 2 arrays per PDE.

And that is simply wrong for a transient computation. Your code is valid only for the steady condition.

[naffrancois]

As other suggested before, CFD is not just a matter of coding finite differences disregarding the physics, the maths and stability constraints. I recommend you, for compressible flows, to consider the conservative form of the equations and a finite volume formulation. Conservative because there are strong theorems to deal with the inherent non-linearities of compressible flows (shocks, contacts) and finite volume because it is the most natural and simplest way of ensuring conservative nature of the equations. Here are two references :

- Toro: "Riemann Solvers and Numerical Methods for Fluid Dynamics" as suggested by Eicosane, for an introduction to the Riemann solution of the Euler equations and its approximations.
- Blazek: "Computational Fluid Dynamics: Principles and Applications" for a practical introduction to finite volume methods on structured and unstructured meshes