Other Equation of State in Laghos about laghos

[vladotomov]

Laghos can be used with different EOS, but the code should be modified at a few places, because it assumes ideal gas. Some parts of the code also assume that gamma is constant in each element.

It's difficult to say what's wrong from these pictures. Is this the initial condition or an evolved solution? Why does the energy look as if it hasn't changed?

I'd recommend to start with running the problem in 1D.

[Liu0813]

Hi, @vladotomov,
Thank you for your suggestion. To test different EOS, I try to solve the water shock tube problem in the following publication:

Chandran, J., & Salih, A. (2019). A modified equation of state for water for a wide range of pressure and the concept of water shock tube. Fluid Phase Equilibria, 483, 182-188.

The initial problem setup is asshown in the following figure, where inside a 3m long tube, a diaphragm separates liquid water under very high pressure from water at a comparatively low pressure. The liquid water is treated as compressible fluid by the relevant EOS, here the Noble-Abel Stiffened-Gas equation of state is used.

In order to solve this problem, the source file laghos_solver.cpp is modified, as follows. The calculation of EOS

Laghos/laghos_solver.cpp

Lines 1007 to 1008 in 95f2b9d

	const double P = (gamma - 1.0) * R * E;
	const double S = sqrt(gamma * (gamma - 1.0) * E);

is changed into

const double P = (gamma - 1.0) * (E + 1177788) / (1/R - 6.72125e-4) - gamma * 6.2178e8;  
const double S = sqrt(gamma * (1/R)*(1/R) * (P+6.2178e8)/(1/R-6.72125e-4));

Correspondingly, the initial condition is set in function rho0, gamma, v0, e0, and the modified code snippet is listed.

double rho0(const Vector &x)
{
   return (x(0) < 2.0) ? 1010 : 997;
}

double gamma(const Vector &x)
{
      return 1.19; 
}

void v0(const Vector &x, Vector &v)
{
      v = 0.0; 
}

double e0(const Vector &x)
{
   return 3610*300+(1/rho0(x)-6.72125e-4)*6.2178e8-1177788;
}

The result is basically reasonable, but there are many non-physical oscillations, like the following figures. Under this initial condition, the amplitude of the wave is relatively small, so the energy change is not obvious.

According to your advice, I have also run this problem in 1D, and non-physical oscillations are not very obvious, like the following figures.

By the way, if the initial density in left region is increased, the effect of non-physical oscillations will also be reduced.

[vladotomov]

Are you running with -pa or -fa? For the 1D test you're surely running full assembly, because the code forces it in 1D.

The full assembly pressure calculation is in ComputeMaterialProperties(), see here.
Update that and let me know if the results are better.

[Liu0813]

Hi, @vladotomov
The above results in 1D and 2D are obtained using -fa mode. For 1D test, pressure and sound speed is calculated in ComputeMaterialProperties(), and I have modified it into my EOS. In addition, I also try the -pa option in 2D problem, and the result is similar.

[vladotomov]

Good, then the EOS calculations are probably correct.
The 1D test looks reasonable. Some oscillations are expected, e.g., for the density and energy in the contact discontinuities. Is it converging to the reference solution under refinement?

[Liu0813]

Hi, @vladotomov
I have tested this problem at different levels of refinement. In the 1D case, as the mesh is refined, the oscillation near the contact discontinuity will decrease, and the result will converge to the reference solution. However, in the 2D case, the refinement of the mesh will cause the oscillation to be worse.

[vladotomov]

Yes surely there is some setup issue in 2D.

[vladotomov]

Are you able to run the 1D problem in 2D? This may show some inconsistencies.

[Liu0813]

I try to run the 1D problem in 2D. The calculation area is divided into 60 elements in the x direction, and there is only one elements in the y direction. It is shown that the size of element in y direction have effect on the result. When dy=0.05 (i.e. dy=dx), the numerical oscillations near the contact discontinuity are very large, as the following figure shown.

But when dy=0.1(or larger), the oscillation will decrease and have little effect on the result.

[vladotomov]

I can't guess what's wrong, but you can use these differences to debug the issue with 2D. To get equivalent 1D and 2D behavior, make sure that the boundary attributes on the 2D mesh are correct, and the initial time steps are the same (you'll have to edit the code to force this for the 2D runs).

[Liu0813]

This problem has been solved. It should be that the default time step is not suitable. When the "cfl" number is decreased, the numerical oscillations will vanish.