
[Sponsors] 
October 7, 2023, 04:47 
Compressible flow solver diverging.

#1 
New Member
Join Date: Oct 2017
Posts: 18
Rep Power: 9 
Recently, I have started developing my own CFD solver (hobby). I have had excellent success with incompressible 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++ (opensource 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(i1, j) .* u(i1, j) .* u(i+1, j)) ...  dt / (2 * h) * (rho(i, j+1) .* u(i, j+1) .* v(i, j+1)  rho(i, j1) .* u(i, j1) .* v(i, j1))...  dt / (2 * h) * (p(i+1, j)  p(i1, j)) ... + nu * (dt / h^2 * (u(i+1, j)  2 * u(i, j) + u(i1, j)) + dt / h^2 * (u(i, j+1)  2 * u(i, j) + u(i, j1)))) ./ rho(i, j); % xmomentum 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. 

October 7, 2023, 10:24 

#2 
Senior Member
andy
Join Date: May 2009
Posts: 301
Rep Power: 18 
You haven't introduced a mechanism to control pressurevelocity 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.


October 7, 2023, 10:29 

#3  
New Member
Join Date: Oct 2017
Posts: 18
Rep Power: 9 
Quote:
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. 

October 7, 2023, 12:32 

#4  
Senior Member
andy
Join Date: May 2009
Posts: 301
Rep Power: 18 
Quote:


October 7, 2023, 12:42 

#5 
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,849
Rep Power: 73 
I would add that a conservative FV scheme is not implemented this way! You have to compute the fluxes.


October 7, 2023, 12:45 

#6 
New Member
Join Date: Oct 2017
Posts: 18
Rep Power: 9 

October 9, 2023, 14:09 

#7  
New Member
Maksim
Join Date: Jul 2020
Posts: 10
Rep Power: 6 
Quote:
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. 

October 9, 2023, 15:52 

#8 
New Member
Join Date: Oct 2017
Posts: 18
Rep Power: 9 
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 incompressible 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 xmomentum I use: u(i, j) = u(i, j)  u(i, j) * dt / (2 * h) .* (u(i+1, j)  u(i1, j))  v(i, j) * dt / (2 * h) .* (u(i, j+1)  u(i, j1))  dt / (2 * rho * h) * (p(i+1, j)  p(i1, j)) + nu * (dt / h^2 * (u(i+1, j)  2 * u(i, j) + u(i1, j)) + dt / h^2 * (u(i, j+1)  2 * u(i, j) + u(i, j1))); % xmomentum all u 's 

October 9, 2023, 16:02 

#9  
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,849
Rep Power: 73 
Quote:
You are using a FD discretization on the quasilinear 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) 

October 9, 2023, 16:13 

#10 
New Member
Join Date: Oct 2017
Posts: 18
Rep Power: 9 
Correction from previous reply:
Here is incompressible xmomentum I use: u(i, j) = u(i, j)  u(i, j) * dt / (2 * h) .* (u(i+1, j)  u(i1, j))  v(i, j) * dt / (2 * h) .* (u(i, j+1)  u(i, j1))  dt / (2 * rho * h) * (p(i+1, j)  p(i1, j)) + nu * (dt / h^2 * (u(i+1, j)  2 * u(i, j) + u(i1, j)) + dt / h^2 * (u(i, j+1)  2 * u(i, j) + u(i, j1))); % xmomentum I am learning finite volume method from Malalasekera CFD book. This blog code [https://fluiddynamicscomputer.blogsp...labcode.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. 

October 9, 2023, 16:14 

#11 
New Member
Join Date: Oct 2017
Posts: 18
Rep Power: 9 
2 distinct arrays for each PDE makes code 510% slower with negligible difference in accuracy, hence I stopped using 2 arrays per PDE.


October 9, 2023, 16:24 

#12 
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,849
Rep Power: 73 

October 9, 2023, 19:24 

#13 
Senior Member
Join Date: Oct 2011
Posts: 242
Rep Power: 17 
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 nonlinearities 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 

Tags 
compreesible codes, compressible, diverging, navier stokes, solver 
Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Ncrit for a glider Xfoil. How to use it. GPT4 answer  AlanMattanó  Main CFD Forum  0  April 10, 2023 13:16 
fluent divergence for no reason  sufjanst  FLUENT  2  March 23, 2016 17:08 
compressible unsteady high Re and Temperature Flow: which solver to use?  loreasr  OpenFOAM Running, Solving & CFD  0  November 25, 2015 10:15 
Solver selection for a compressible flow problem  houkensjtu  OpenFOAM Running, Solving & CFD  1  July 31, 2015 07:53 
Solver for steady state compressible multispecies nonreacting flow?  inf.vish  OpenFOAM Running, Solving & CFD  0  October 1, 2013 02:09 