Hello Daniel,
Thanks for the answer
An implementation of the Roe scheme would certainly go a long way in
simplifying the resolution of hyperbolic equations but I think FiPy already
has was it takes to do that - terms wise.
I know of CLAWPACK yes, but I really would like to see if FiPy can handle
fluid dynamics equations. I work at CERFACS, a lab in France centered
around the use of Computational Fluid Dynamics for academic and industrial
purpose, and although we do have a production code, I would like to let my
colleagues know about FiPy because I think it has the potential to be a
great first-approach solver for simple to mildly complex problems.
I actually already succeeded in solving (and validating with the known Sod
shock tube test case) the 1D compressible Euler equations (script attached
with an exact solver for reference).
But then in 1D, things are pretty easy - equation wise.
If you have time to check the code, can you tell me if you see any kind of
mis-use of FiPy ? Are there any lines you would have written differently ?
In 2D everything becomes harder and I can’t get it right - although I
suspect it is because I don’t know FiPy very well yet.
If you don’t mind, I have a couple questions to try and make it work.
1/ I do not get how to include the resolution of an equation of state in
the system of differential equations. I wrote a stub of code (attached,
called 7-covo2D.py) in which I solve for (rho, rhoU, rhoV, rhoE) as usual.
An easy way to write the Euler equations however is to include a notion of
pressure, which is related to the unknown by an equation of state (it is
not a differential equation though) - see for instance
http://bender.astro.sunysb.edu/hydro_by_example/compressible/Euler.pdf.
Because of this, I have to inject the EoS directly into the equations and
try and make do with that, but the Finite Volume method is not designed to
have the fluxes computed like I do in the script, so naturally everything
crashes or gives nonsensical results.
Do you see a possibility in FiPy to solve 4 differential equations and 1
algebraic equation (the EoS) in a coupled manner ?
2/ If it is not possible to include the equation of state in the system,
how would you take the -dp/dx and the -dp/dy from the momentum equations
(with p expanded as a function of rhoE, rhoU and rhoV) into account ?
Thank you very much for your help,
Best regards,
Thibault
2017-05-15 15:38 GMT+02:00 Daniel Wheeler :
> Hi Thibault,
>
> I started down the road of implementing a Riemann solver in FiPy, see
>
> https://github.com/usnistgov/fipy/blob/riemann/fipy/terms/
> baseRoeConvectionTerm.py
>
> Unfortunately, I never completed and merged this work. You might want
> to look into CLAWPACK for a code to better handle compressible flow.
>
> I've given a few answers about this in the past as well, see
>
>- http://git.net/ml/python-fipy/2012-02/msg00024.html
>
>- http://fipy.nist.narkive.com/A0gJrSl2/semilinear-wave-
> equation-in-fipy
>
>- http://fipy.nist.narkive.com/YVZTRM0G/1d-coupled-fluid-
> equations-in-fipy
>
> Cheers,
>
> Daniel
>
> On Sun, May 14, 2017 at 9:33 AM, Thibault Bridel-Bertomeu
> wrote:
> > Good afternoon,
> >
> > I was wondering if anyone had ever tried implementing the compressible
> euler
> > equations in FiPy ?
> > I have been trying for a while now but even in 1D I can’t get a proper
> shock
> > tube solution.
> >
> > If anyone has ever tried, I would be infinitely grateful if they were to
> > share their knowledge !!
> >
> > Thank you very much,
> >
> > T. BRIDEL-BERTOMEU
> >
> >
> >
> > ___
> > fipy mailing list
> > fipy@nist.gov
> > http://www.ctcms.nist.gov/fipy
> > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
> >
>
>
>
> --
> Daniel Wheeler
>
> ___
> fipy mailing list
> fipy@nist.gov
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
#!/usr/bin/env python
import fipy
from fipy import numerix
#==
nx = 1000
dx = 1./nx
Lx = nx * dx
gamma = 1.4
gm1 = gamma-1.0
#--
mesh = fipy.Grid1D(nx=nx, dx=dx)
X = mesh.cellCenters[0]
x = mesh.faceCenters[0]
#--
ro = fipy.CellVariable(mesh=mesh, name="ro", hasOld=True)
roU = fipy.CellVariable(mesh=mesh, name="roU", hasOld=True)
roE = fipy.CellVariable(mesh=mesh, name="roE", hasOld=True)
ro.setValue( 1.0*(X<=0.5)+0.125*(X>0.5) )
roU.setValue( 0.0 )
roE.setValue( 1.0/gm1*(X<=0.5)+0.1/gm1*(X>0.5) )
ro.constrain(value=1.0, where=mesh.facesLeft)
ro.constrain(value=0.125, where=mesh.facesRight)
roU.constrain(value=0.0, where=mesh.facesLeft|mesh.facesRight)
roE.constrain(value=1.0/gm1, where=mesh.facesLeft)
roE.constrain(value=0.1/gm1, where=mesh.facesRight)
vi = fipy.Viewer(vars=(ro, roU, roE), datamin=-0.1, datamax=2.6)
vi.plot()
raw_input("Initialization ...")
U = fipy.CellVa
