CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > General Forums > Main CFD Forum

Near vacuum for Euler equations

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree2Likes

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   April 24, 2023, 10:17
Default Near vacuum for Euler equations
  #1
New Member
 
Davide
Join Date: Apr 2023
Posts: 10
Rep Power: 2
dappoli is on a distinguished road
Hi everyone,

I am looking for some help in solving a problem I face in solving the Euler equations (more or less).

I have looked and looked but up to now I have not found a functional solution.

The scenario I am trying to simulate is actually a little bit more complex than the usual inviscid Euler equations:


I am solving a time dependent multifluid 2D problem where I consider a plasma composed by ions and electrons. The two fluids are coupled through the electric field solved with the Poisson equation. Everything is immersed in a uniform magnetic field. To have an idea, each specie (ions and electrons) obeys the Euler equations with non zero sources on the rhs, in particular I have a mass injection, electric and magnetic field effects. I solve for the conservative variables.



To complicate things even more, the scenario I am simulating results in supersonic and subsonic regions, extremely large velocities in some points and very low densities, reaching near vacuum conditions.



I have successfully run this simulation with a Rusanov scheme with 2nd order MUSCL reconstruction and an unsplit scheme, with a third order SSP Runge Kutta integration. The Rusanov fluxes are positive preserving so that I never reach negative densities and everything works (more or less). However, since I am trying to evaluate the occurrence of turbulent effects, the Rusanov flux is too dissipative and the results are "contaminated" by numerical diffusion.



I then coded Riemann solvers utilizing more waves, in particular both the Roe and the HLLEM and I am currently using the last one, with 2nd order MUSCL reconstruction on primitive variables and MINMOD or VanLeer limiters. (the code is a home made solver written in fortran with structure meshes).



With the HLLEM however, given that the numerical diffusion is strongly reduced with respect to Rusanov, the simulation hardly manages to complete, mostly because I obtain a negative temperature under square root. I know that the issue of positivity preserving is problematic, especially dealing with flows where the total energy is close to the kinetic energy, so that computing the temperature may lead to negative values. What I am trying to do is to switch to the more diffusive Rusanov, or even 1st order Rusanov, scheme when I detect negative values, thus recomputing the time step. I even tried using such scheme in regions of the domain where I know I will have problems, near the boundaries, but I get weird effects at the interface between the two schemes and the problem arises again.



The biggest issue seems to be the computation of the temperature itself, density and energies, even pressure are smooth, but when I compute the temperature it is highly discontinuous in low density regions. My feeling is that i suffer from catastrophic cancellation in subtracting the kinetic energy from the total energy and once I divide by the density all the errors come up. I am not sure on how to proceed from here since I cannot think of a way to prevent this. I even thought about having some sort of tracking interface where in one side I have density and in the other vacuum where I don't solve anything, but is this even possible?



I know that I am stretching the fluid approximation itself, but this is exactly the point. Since I have limited experience on the subject I might have lost something.

Is there any solution or technique, switching or trick that I can apply to avoid the temperature to become negative or in general blowing up ?


thanks,

Davide
dappoli is offline   Reply With Quote

Old   April 24, 2023, 15:33
Default
  #2
Senior Member
 
Eifoehn4's Avatar
 
-
Join Date: Jul 2012
Location: Germany
Posts: 184
Rep Power: 13
Eifoehn4 is on a distinguished road
Can you write down some equations of your model? Do you have a one energy/density/velocity model or more transport equations for each species?
__________________
Check out my side project:

A multiphysics discontinuous Galerkin framework: Youtube, Gitlab.
Eifoehn4 is offline   Reply With Quote

Old   April 24, 2023, 16:32
Default
  #3
New Member
 
Davide
Join Date: Apr 2023
Posts: 10
Rep Power: 2
dappoli is on a distinguished road
Sure, here you can see the general form of the system of equations which is solved for each specie.


\frac{\partial n_\alpha}{\partial t} + \nabla \cdot \left(n_\alpha\textbf{u}_\alpha\right) = S_\alpha
\frac{\partial \left(n_\alpha \textbf{u}_\alpha\right)}{\partial t} + \nabla \cdot \left(\left(n_\alpha \textbf{u}_\alpha\right)\textbf{u}_\alpha\right) = \frac{qn_\alpha}{m_\alpha}\left[\textbf{u}_\alpha\times \textbf{B} + \nabla \phi\right] - \frac{1}{m_\alpha}\nabla P_\alpha
\frac{\partial \left(\mathcal{E}_\alpha\right)}{\partial t} + \nabla \cdot \left(\left(\mathcal{E}_\alpha + p_\alpha\right)\textbf{u}_\alpha \right) = qn_\alpha \textbf{u}_\alpha \cdot \nabla \phi + {Q}_\alpha


The system can be easily recast in the "standard" form of the Euler equations.
In this case n is the number density, E is the total energy (not specific), q is the charge (\pm e, the elementary charge), \phi is the electric potential and B = B_z z the uniform magnetic field perpendicular to the domain (x,y).

Moreover Q_\alpha=\mathcal{E}_{inj}S_\alpha is the injection energy times the injection source term.

\alpha denotes either the ions and electrons which are integrated in time separately, advanced of one step in each algorithm loop. The two systems are coupled through the Poisson equation :


\nabla^2\phi = \frac{e}{\epsilon_0}\left(n_e-n_i\right)


Just as a reference, for Xenon (as an example) m_i = 2.18e-25 kg while m_e = 9.1e-31 kg, which leads to a huge difference in timescales for ions and electrons, as well as very different sound speeds (given that electrons have temperatures 1-2 orders of magnitude larger than ions). Moreover in order to solve the non neutrality of the problem, the debye length and plasma frequency must be resolved, which usually forces the timestep to be extremely small. I am currently using between 0.5 and 10 ps.
dappoli is offline   Reply With Quote

Old   April 24, 2023, 17:06
Default
  #4
New Member
 
Maksim
Join Date: Jul 2020
Posts: 10
Rep Power: 5
Eicosane is on a distinguished road
Have You tested your code with Noh test? It's a test, where you can check, how small internal energy can be comparing to kinetic energy in your scheme.

Also what variables do You reconstruct? Density, velocity and pressure?
Eicosane is offline   Reply With Quote

Old   April 24, 2023, 17:34
Default
  #5
New Member
 
Davide
Join Date: Apr 2023
Posts: 10
Rep Power: 2
dappoli is on a distinguished road
Hi, thanks for the reply.
No I didn't try the Noh test, I tested the code with simpler tests such as the Sod test and Sedov blast. I will try.
Yes, I am reconstructing density, velocities and pressure, then from the reconstructed pressure I obtain the temperature for the eigenvalues.

I must have messed something up because I thought I sent a reply with all the equations, but apparently something went wrong. Tomorrow I'll write it again
dappoli is offline   Reply With Quote

Old   April 24, 2023, 17:57
Default
  #6
New Member
 
Maksim
Join Date: Jul 2020
Posts: 10
Rep Power: 5
Eicosane is on a distinguished road
Quote:
Originally Posted by dappoli View Post
Yes, I am reconstructing density, velocities and pressure, then from the reconstructed pressure I obtain the temperature for the eigenvalues.
And when exactly temperature becomes negative? After reconstruction?
Eifoehn4 likes this.
Eicosane is offline   Reply With Quote

Old   April 24, 2023, 18:04
Default
  #7
Senior Member
 
Eifoehn4's Avatar
 
-
Join Date: Jul 2012
Location: Germany
Posts: 184
Rep Power: 13
Eifoehn4 is on a distinguished road
Quote:
Originally Posted by Eicosane View Post
And when exactly temperature becomes negative? After reconstruction?
Could be one of many reasons. With density-pressure reconstruction your temperature is not TVD.
__________________
Check out my side project:

A multiphysics discontinuous Galerkin framework: Youtube, Gitlab.
Eifoehn4 is offline   Reply With Quote

Old   April 24, 2023, 19:23
Default
  #8
New Member
 
Maksim
Join Date: Jul 2020
Posts: 10
Rep Power: 5
Eicosane is on a distinguished road
Quote:
Originally Posted by Eifoehn4 View Post
Could be one of many reasons. With density-pressure reconstruction your temperature is not TVD.
But from my experience, on contrary, density-velocity-pressure reconstruction instead of density-momentum-total energy allows to avoid negative temperature in some cases, for example Woodward and Collela two blastwave test.
Eicosane is offline   Reply With Quote

Old   April 24, 2023, 20:04
Default
  #9
Senior Member
 
Eifoehn4's Avatar
 
-
Join Date: Jul 2012
Location: Germany
Posts: 184
Rep Power: 13
Eifoehn4 is on a distinguished road
Quote:
Originally Posted by Eicosane View Post
But from my experience, on contrary, density-velocity-pressure reconstruction instead of density-momentum-total energy allows to avoid negative temperature in some cases, for example Woodward and Collela two blastwave test.
The problem is more severe as you think and already occurs considering the spatial discretization, only. The Euler equations consist of three evolution equations for mass, momentum and energy. Since, there are four unknowns - density, velocity, energy and pressure - you need a thermodynamic law - the caloric EoS - to calculate one from the other.

One important thing you have to consider is that a thermodynamic quantity is uniquely linked to two other thermodynamic quantities. Moreover, you are using a linear interpolation during the TVD reconstruction.

You have two choices for the facial values:

Either you violate thermodynamic consistency between thermodynamic quantities or you stick to it however, violate the TVD property for all other thermodynamic quantities.

This is the price you have to pay for numerical methods higher than first order, see Godunov theorem.

It would be helpful to write down the equations, otherwise it is hard to tell you what might be the problem.
__________________
Check out my side project:

A multiphysics discontinuous Galerkin framework: Youtube, Gitlab.
Eifoehn4 is offline   Reply With Quote

Old   April 25, 2023, 02:50
Default
  #10
New Member
 
Davide
Join Date: Apr 2023
Posts: 10
Rep Power: 2
dappoli is on a distinguished road
These are the equations:

\frac{\partial n_\alpha}{\partial t} + \nabla \cdot \left(n_\alpha\textbf{u}_\alpha\right) = S_\alpha
\frac{\partial \left(n_\alpha \textbf{u}_\alpha\right)}{\partial t} + \nabla \cdot \left(\left(n_\alpha \textbf{u}_\alpha\right)\textbf{u}_\alpha\right) = \frac{qn_\alpha}{m_\alpha}\left[\textbf{u}_\alpha\times \textbf{B} + \nabla \phi\right] - \frac{1}{m_\alpha}\nabla P_\alpha
\frac{\partial \left(\mathcal{E}_\alpha\right)}{\partial t} + \nabla \cdot \left(\left(\mathcal{E}_\alpha + p_\alpha\right)\textbf{u}_\alpha \right) = qn_\alpha \textbf{u}_\alpha \cdot \nabla \phi + {Q}_\alpha
which are integrated for electrons and ions where q is the charge (\pm e, the fundamental charge), B the magnetic field, in my case perpendicular to the domain,
\mathcal{E}_\alpha = \frac{3}{2}q n_\alpha T_\alpha + \frac{1}{2}m_\alpha n_\alpha u_\alpha^2 the total energy ( Note that T is expressed in eV here), Q_\alpha = \mathcal{E}_{inj}S_\alpha the energy injection source term and \phi the elctric potential obtained from:



\nabla^2\phi = \frac{e}{\epsilon_0}\left(n_e-n_i\right)


At the moment I am trying also with a Strang splitting, so that I solve for the homogeneous system which is identical to Euler. The eigenvalues are the usual eigenvalues. In the algorithm I advance one dt the single species and then compute the new potential for the next timestep.

If you are not familiar with plasma, just as a reference, the problem is complicated by the broad range of timescales. For instance, for Xenon m_i = 2.18e-25 kg and m_e = 9.1e-31 kg, so that there is a huge difference in the timescales of ions and electrons. Thermal velocities for electrons are in the order of 1e6 m/s and the temperature is 1-2 order of magnitude larger than the one of ions. Moreover, to solve the non-neutrality of the problem, the debye length and plasma frequency must be resolved, limiting the problem to very fine meshes and timesteps in the 0.5-10 ps order.



I know that in the pressure-density-velocity the TVD is not granted, but more in general the positivity of schemes such as the Rusanov is not granted either for multidimensional system with reconstruction > order 1. The negative temperature sometimes occurs in the face reconstruction, but that is easily solvable (or at least more easily with a dedicated limiter), but sometimes it just pops up in the cell centre. After I advance in time the conservative variables, recovering the temperature in the cell centre my give rise to a negative value. Note that in this problem it is not necessary to know the temperature at cell centre, but if one has to compute collision frequencies ecc. it is mandatory.



I know that there is a work from Linde and Roe where multidimensional schemes of order >1 are made positive preserving using a dedicated limiter, which however is somehow not straightforward for the pressure. Indeed they state that the pressure must be limited iteratively to guarantee the positivity. I would try to avoid this given the fact that the timestep is already incredibly small and time can be an issue.



https://link.springer.com/chapter/10...-011-5169-6_16


By the way, does it make sense for you that i can see a discontinuity in the solution (again in the temperature, but I guess only because I am dividing by the discontinuous density) where I use 1st order rusanov instead of 2nd order HLLEM? Should the interface be smooth? Because I might have an error in the code if that is the case
dappoli is offline   Reply With Quote

Old   April 25, 2023, 04:00
Default
  #11
Senior Member
 
sbaffini's Avatar
 
Paolo Lampitella
Join Date: Mar 2009
Location: Italy
Posts: 2,152
Blog Entries: 29
Rep Power: 39
sbaffini will become famous soon enoughsbaffini will become famous soon enough
Send a message via Skype™ to sbaffini
Typically, if the error is in the roundoff, and the code has been properly coded, you can test it by temporarily switching to a higher precision and see if the problem is still there.
sbaffini is offline   Reply With Quote

Old   April 25, 2023, 04:51
Default
  #12
New Member
 
Davide
Join Date: Apr 2023
Posts: 10
Rep Power: 2
dappoli is on a distinguished road
Yes I tired doing that, but I am already in double precision and MUMPS (the parallel linear solver I am using) does not work with quadruple precision, so in the end I abandoned the test.
dappoli is offline   Reply With Quote

Old   April 25, 2023, 05:27
Default
  #13
Senior Member
 
Eifoehn4's Avatar
 
-
Join Date: Jul 2012
Location: Germany
Posts: 184
Rep Power: 13
Eifoehn4 is on a distinguished road
Quote:
Originally Posted by dappoli View Post
These are the equations:

\frac{\partial n_\alpha}{\partial t} + \nabla \cdot \left(n_\alpha\textbf{u}_\alpha\right) = S_\alpha
\frac{\partial \left(n_\alpha \textbf{u}_\alpha\right)}{\partial t} + \nabla \cdot \left(\left(n_\alpha \textbf{u}_\alpha\right)\textbf{u}_\alpha\right) = \frac{qn_\alpha}{m_\alpha}\left[\textbf{u}_\alpha\times \textbf{B} + \nabla \phi\right] - \frac{1}{m_\alpha}\nabla P_\alpha
\frac{\partial \left(\mathcal{E}_\alpha\right)}{\partial t} + \nabla \cdot \left(\left(\mathcal{E}_\alpha + p_\alpha\right)\textbf{u}_\alpha \right) = qn_\alpha \textbf{u}_\alpha \cdot \nabla \phi + {Q}_\alpha
which are integrated for electrons and ions where q is the charge (\pm e, the fundamental charge), B the magnetic field, in my case perpendicular to the domain,
\mathcal{E}_\alpha = \frac{3}{2}q n_\alpha T_\alpha + \frac{1}{2}m_\alpha n_\alpha u_\alpha^2 the total energy ( Note that T is expressed in eV here), Q_\alpha = \mathcal{E}_{inj}S_\alpha the energy injection source term and \phi the elctric potential obtained from:



\nabla^2\phi = \frac{e}{\epsilon_0}\left(n_e-n_i\right)


At the moment I am trying also with a Strang splitting, so that I solve for the homogeneous system which is identical to Euler. The eigenvalues are the usual eigenvalues. In the algorithm I advance one dt the single species and then compute the new potential for the next timestep.

If you are not familiar with plasma, just as a reference, the problem is complicated by the broad range of timescales. For instance, for Xenon m_i = 2.18e-25 kg and m_e = 9.1e-31 kg, so that there is a huge difference in the timescales of ions and electrons. Thermal velocities for electrons are in the order of 1e6 m/s and the temperature is 1-2 order of magnitude larger than the one of ions. Moreover, to solve the non-neutrality of the problem, the debye length and plasma frequency must be resolved, limiting the problem to very fine meshes and timesteps in the 0.5-10 ps order.



I know that in the pressure-density-velocity the TVD is not granted, but more in general the positivity of schemes such as the Rusanov is not granted either for multidimensional system with reconstruction > order 1. The negative temperature sometimes occurs in the face reconstruction, but that is easily solvable (or at least more easily with a dedicated limiter), but sometimes it just pops up in the cell centre. After I advance in time the conservative variables, recovering the temperature in the cell centre my give rise to a negative value. Note that in this problem it is not necessary to know the temperature at cell centre, but if one has to compute collision frequencies ecc. it is mandatory.



I know that there is a work from Linde and Roe where multidimensional schemes of order >1 are made positive preserving using a dedicated limiter, which however is somehow not straightforward for the pressure. Indeed they state that the pressure must be limited iteratively to guarantee the positivity. I would try to avoid this given the fact that the timestep is already incredibly small and time can be an issue.



https://link.springer.com/chapter/10...-011-5169-6_16


By the way, does it make sense for you that i can see a discontinuity in the solution (again in the temperature, but I guess only because I am dividing by the discontinuous density) where I use 1st order rusanov instead of 2nd order HLLEM? Should the interface be smooth? Because I might have an error in the code if that is the case

No problem, I am familiar with plasma.

Currently I am not quite sure why temperature is explictly needed in your system during time advancement. Your model is closed with a caloric EoS, only. Temperature seems to be needed only for postprocessing.

How is pressure connected to internal energy? Via standard thermodnymic relations?

Regards
__________________
Check out my side project:

A multiphysics discontinuous Galerkin framework: Youtube, Gitlab.
Eifoehn4 is offline   Reply With Quote

Old   April 25, 2023, 06:04
Default
  #14
Senior Member
 
sbaffini's Avatar
 
Paolo Lampitella
Join Date: Mar 2009
Location: Italy
Posts: 2,152
Blog Entries: 29
Rep Power: 39
sbaffini will become famous soon enoughsbaffini will become famous soon enough
Send a message via Skype™ to sbaffini
Quote:
Originally Posted by dappoli View Post
Yes I tired doing that, but I am already in double precision and MUMPS (the parallel linear solver I am using) does not work with quadruple precision, so in the end I abandoned the test.
Oh, ok... altough, on a second tought, if you already identified the offending line in your code and the reason behind it, there is little to be gained from that. I mean, if legit input then gives you a wrong output on a certain line, you already know what the problem is.
sbaffini is offline   Reply With Quote

Old   April 25, 2023, 06:38
Default
  #15
New Member
 
Davide
Join Date: Apr 2023
Posts: 10
Rep Power: 2
dappoli is on a distinguished road
Quote:
Originally Posted by Eifoehn4 View Post

Currently I am not quite sure why temperature is explictly needed in your system during time advancement. Your model is closed with a caloric EoS, only. Temperature seems to be needed only for postprocessing.

How is pressure connected to internal energy? Via standard thermodnymic relations?

Regards

You are right that in this particular problem the temperature at cell centre is not used by the code anywhere. Thus I guess the failure mechanism when I see negative temperature in the dataset is:

advance one step -> compute temperature for post processing and get T<0 -> at the next computation also the reconstructed temperature at the face is going to be negative and then it fails on c = sqrt(5/3*e*T/m)



so yeah eventually I guess it is always the face temperature that makes the code fail.

Yes I use standard relations assuming scalar pressure
dappoli is offline   Reply With Quote

Old   April 25, 2023, 07:47
Default
  #16
Senior Member
 
sbaffini's Avatar
 
Paolo Lampitella
Join Date: Mar 2009
Location: Italy
Posts: 2,152
Blog Entries: 29
Rep Power: 39
sbaffini will become famous soon enoughsbaffini will become famous soon enough
Send a message via Skype™ to sbaffini
There are simple bechmark shock tube cases that would only work for positivity preserving schemes. Have tou tried your first and second order hllem scheme against those?
sbaffini is offline   Reply With Quote

Old   April 25, 2023, 08:01
Default
  #17
Senior Member
 
Eifoehn4's Avatar
 
-
Join Date: Jul 2012
Location: Germany
Posts: 184
Rep Power: 13
Eifoehn4 is on a distinguished road
Quote:
Originally Posted by dappoli View Post
You are right that in this particular problem the temperature at cell centre is not used by the code anywhere. Thus I guess the failure mechanism when I see negative temperature in the dataset is:

advance one step -> compute temperature for post processing and get T<0 -> at the next computation also the reconstructed temperature at the face is going to be negative and then it fails on c = sqrt(5/3*e*T/m)



so yeah eventually I guess it is always the face temperature that makes the code fail.

Yes I use standard relations assuming scalar pressure
Maybe here is the problem: A proper implementation of the convective part should never generate negative temperatures in the cell center. You may violate some stability constraints including the source terms.

As Paolo already said. It may be useful to check some standard shock-tube examples. Using a TVD reconstruction + special Riemann solvers, often violate some entropy constraints.

Regards
__________________
Check out my side project:

A multiphysics discontinuous Galerkin framework: Youtube, Gitlab.
Eifoehn4 is offline   Reply With Quote

Old   April 25, 2023, 09:22
Default
  #18
New Member
 
Davide
Join Date: Apr 2023
Posts: 10
Rep Power: 2
dappoli is on a distinguished road
Quote:
Originally Posted by Eifoehn4 View Post
Maybe here is the problem: A proper implementation of the convective part should never generate negative temperatures in the cell center. You may violate some stability constraints including the source terms.

Isn't obtaining T<0 at cell centre possible because the scheme is not positive preserving? Since my conserved variable which i solve for is the total energy, if the scheme is not positive conservative or the condition for which it is are violated (for instance cfl <0.5 for Rusanov) then I can get T<0.

But this is beyond the numerical implementation of the code, but rather an issue due to the linearization of the Riemann problem performed by the approximate solver, correct me if I am wrong.


Anyway I tried doing the Noh test, very quickly so I might have an issue with the boundary condition. Since wall condition are not imposed in my code, I solved for all the cylindrical geometry and I put 1st order outflow conditions at the external boundaries (instead of enforcing the analytical solution).

With the above system of equation, with the source terms = 0 (exception made for pressure gradient which is part of the flux in my code) :
m = 1e15 kg

n = 1e-15 m^-3
p = 1e-6 Pa

domain 1x1 cm and 800x800 cells

max cfl in both directions = 0.1



attached you can see the comparison of Rusanov and HLLEM with 2nd order MUSCL. HLLEM behaves strangely at the origin and in general inside the shock so there might be some issue in the code I guess, but what could that be?

Moreover, by looking at the density plot, the square shape originating far from the shock is due to the boundary condition, being a 1st order extrapolation reflects something I guess.

In any case after a while the shock collapses in a weird shape, with both schemes ,as reported in one of the picture. Is this normal? I never saw the "long" time solution of the Noh problem.


Quote:
Originally Posted by sbaffini View Post
There are simple bechmark shock tube cases that would only work for positivity preserving schemes. Have tou tried your first and second order hllem scheme against those?
sorry I did not see this, no I did not try them.. I could try the double rarefaction wave and the hurricane-like simulation with vacuum in the centre. I will let you know
Attached Images
File Type: jpg n_e_ok.jpg (57.1 KB, 26 views)
File Type: jpg n_e_end.jpg (65.9 KB, 20 views)
File Type: jpg Density_1D_Comparison.jpg (50.6 KB, 47 views)
dappoli is offline   Reply With Quote

Old   April 25, 2023, 13:09
Default
  #19
Senior Member
 
Eifoehn4's Avatar
 
-
Join Date: Jul 2012
Location: Germany
Posts: 184
Rep Power: 13
Eifoehn4 is on a distinguished road
Quote:
Originally Posted by dappoli View Post
Isn't obtaining T<0 at cell centre possible because the scheme is not positive preserving? Since my conserved variable which i solve for is the total energy, if the scheme is not positive conservative or the condition for which it is are violated (for instance cfl <0.5 for Rusanov) then I can get T<0.

But this is beyond the numerical implementation of the code, but rather an issue due to the linearization of the Riemann problem performed by the approximate solver, correct me if I am wrong.


Anyway I tried doing the Noh test, very quickly so I might have an issue with the boundary condition. Since wall condition are not imposed in my code, I solved for all the cylindrical geometry and I put 1st order outflow conditions at the external boundaries (instead of enforcing the analytical solution).

With the above system of equation, with the source terms = 0 (exception made for pressure gradient which is part of the flux in my code) :
m = 1e15 kg

n = 1e-15 m^-3
p = 1e-6 Pa

domain 1x1 cm and 800x800 cells

max cfl in both directions = 0.1



attached you can see the comparison of Rusanov and HLLEM with 2nd order MUSCL. HLLEM behaves strangely at the origin and in general inside the shock so there might be some issue in the code I guess, but what could that be?

Moreover, by looking at the density plot, the square shape originating far from the shock is due to the boundary condition, being a 1st order extrapolation reflects something I guess.

In any case after a while the shock collapses in a weird shape, with both schemes ,as reported in one of the picture. Is this normal? I never saw the "long" time solution of the Noh problem.



sorry I did not see this, no I did not try them.. I could try the double rarefaction wave and the hurricane-like simulation with vacuum in the centre. I will let you know
I am not quite sure what you mean by positive conservative. A monoton and positive preserving scheme of course should use a limiter without anti-dissipation, MINMOD and Van Leer are ok. Moreover, you should choose a suited Riemann solver, Rusanov and HLLEM are also ok. A positive value then depends only on the CFL condition and the temporal discretization.

Both of your solutions, Rusanov and HLLEM, look quite oscillatory. My simulations do not contain these little wiggles, even if I use an HLLC instead of a HLLEM. The collapse is due to your boundary conditions and seems ok.

Regards
Attached Images
File Type: png Noh.png (59.6 KB, 41 views)
__________________
Check out my side project:

A multiphysics discontinuous Galerkin framework: Youtube, Gitlab.
Eifoehn4 is offline   Reply With Quote

Old   April 26, 2023, 04:37
Default
  #20
New Member
 
Davide
Join Date: Apr 2023
Posts: 10
Rep Power: 2
dappoli is on a distinguished road
Quote:
Originally Posted by Eifoehn4 View Post
I am not quite sure what you mean by positive conservative. A monoton and positive preserving scheme of course should use a limiter without anti-dissipation, MINMOD and Van Leer are ok. Moreover, you should choose a suited Riemann solver, Rusanov and HLLEM are also ok. A positive value then depends only on the CFL condition and the temporal discretization.

First of all thank you for the reply.

I think we are talking about the same thing, for me a positive preserving scheme is a scheme that generates "physically admissible solutions from physical data" [1], or in other words it keeps positive pressure and density if the initial condition are positive. For the Euler equations the physically admissible state in this case is defined as:

G = \{ \textbf{U} \vert \rho > 0 \quad \text{and} \quad E-(m_x^2 + m_y^2 + m_z^2)/2 \rho >0 \}
being U the conserved variable vector and m the momentum along the 3 directions. For example Roe scheme is not positive preserving and neither is Osher's. Einfeldt [2] demonstrated that the Rusanov and HLLE are positive preserving, but for example HLLEM is not. All this holds true for 1st order schemes, once 2nd order schemes are taken into account, from what I found there is no proof of satisfying the positive preserving property. That is why the Muscl reconstruction near low density fails also for Rusanov. Clearly if the timestep is strongly reduced this is in general postponed, but inevitable. Having said this, I am new to the study and application of Riemann solvers so I might have said something wrong.



Going back to Noh:

Quote:
Both of your solutions, Rusanov and HLLEM, look quite oscillatory. My simulations do not contain these little wiggles, even if I use an HLLC instead of a HLLEM. The collapse is due to your boundary conditions and seems ok.
You are right, but being focused on other problems, I forgot that I am using a variation of the MINMOD slope limiter described in Kurganov-Tadmor [3] (the scheme I am using which is equal to a MUSCL Rusanov) defined as:

u_x(x,t)\approx \text{minmod}\left(\theta \frac{u(x+a)-u(x)}{a}, \frac{u(x+a)-u(x-b)}{a+b}, \theta \frac{u(x)-u(x-b)}{b}\right)
which allows to reduce the diffusivity of the minmod by tuning the parameter \theta\in [1,2]. Setting the parameter to 1 recovers the original minmod. I had it at 1.4. Attached you can find the solution with the parameter = 1. Now the solution is sharp, although the HLLEM keeps behaving strangely at the centre of the domain (note that the simulation is the same, I just cut the plot). Do you have any idea of why it behaves like that at the origin?



Finally, I coded the "hurricane" test described in [4] where a uniform domain is initialized with a vortex like velocity. This creates a vacuum point in the centre of the domain and if the scheme is not positive preserving, eventually it will fail. Of course 2nd order MUSCL Rusanov, HLLE and HLLEM fail in this, but first order HLLE and Rusanov do not. What I did in the end is to use a posteriori approach:

I advance 1 timestep,

I check for negative pressure, temperature or NANs ,

If I detect errors, I reject the step and re do it with 1st order Rusanov in the troubled cell and its neighbours.



I already tried doing this with little success, but I was not considering the neighbouring cells and just using the Rusanov flux in the troubled one. This eventually does not work. Now I will have to deal with some special cases such as when the troubled cell is on the boundary between processors (I use domain decomposition with MPI).

Attached you can see a video of the hurricane like simulation with the HLLEM a polsteriori corrected scheme. Sorry for the low quality, the solution is smooth in reality but the forum has a very low size limit



Now I will have to try this on my problem, probably with a careful optimization of the minmod limiter and neighbour radius to avoid failure.


[1] https://link.springer.com/chapter/10...-011-5169-6_16
[2] https://www.sciencedirect.com/scienc...21999191902113
[3] https://www.math.umd.edu/~tadmor/pub...or.JCP-00I.pdf
[4] https://www.math.hkust.edu.hk/~makxu...D-inviscid.pdf
Attached Images
File Type: jpg Density_1D_Comparison.jpg (54.4 KB, 26 views)
File Type: gif ezgif.com-optimize.gif (102.1 KB, 22 views)
dappoli is offline   Reply With Quote

Reply

Tags
euler, plasma, vacuum

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Solving coupled scalar transport equations in openFoam pavaninguva OpenFOAM Running, Solving & CFD 1 December 28, 2019 08:28
Calculation of the Governing Equations Mihail CFX 7 September 7, 2014 06:27
Riemann invariants of adjoint equations of shallow water equations zqb0929 Main CFD Forum 0 March 15, 2012 00:54
CFD governing equations m.gos Main CFD Forum 0 April 30, 2011 14:21
? fluctuating equations for homogenous shear turb. ff_fan Main CFD Forum 1 September 20, 2002 07:39


All times are GMT -4. The time now is 14:59.