Changes between Version 1 and Version 2 of u/adebrech/code

Timestamp:

03/15/17 17:13:29 (8 years ago)

Author:

adebrech

Comment:

—

u/adebrech/code

@@ -1 @@
-= Exact Riemann Solver =
+= Exact Riemann Solver (no vacuum) =
 
+The Riemann problem is the solution of a system of hyperbolic conservation laws (in our case the Euler equations). This solver assumes two one-dimensional ideal gases separated at the origin x=0, with initial densities, velocities, and pressures given as input. It gives a solution to the set of equations
+
+{{{#!latex
+$\rho + \rho u = 0$
+
+$\rho u + \rho u^2 + p = 0$
+
+$E + u(E + p) = 0$
+}}}
+
+with initial conditions
+
+
+{{{#!latex
+$\rho(x,0) = \begin{cases} \rho_L, & x < 0 \\ \rho_R, & x > 0 \end{cases}$
+
+$\rho(x,0)u(x,0) = \begin{cases} \rho_L u_L, & x < 0 \\ \rho_R u_R, & x > 0 \end{cases}$
+
+$E(x,0) = \begin{cases} E_L, & x < 0 \\ E_R, & x > 0 \end{cases}$
+}}}
+
+The solutions to this IVP at later times are split into four regions - to the left, the initial left state, unperturbed by interaction with the right; between the left and center-left regions, there will either be a rarefaction fan or a shock (assumed ideal), across which the new equilibrium left state holds; the left and right fluids will be separated by a contact discontinuity, across which the new equilibrium right state holds (note that the right and left equilibrium states differ only in their densities); the center-right region will then be separated by a rarefaction fan or (again ideal) shock from the rightmost region, which is still in the initial right state.
+
+The program reads input from a file called problem.data, and stores many of the relevant parameters in global variables.
+
+{{{#!fortran
+NAMELIST/ProblemData/ rhoL,uL,pL,rhoR,uR,pR,gamma,tol,starttime,finaltime,frames,xL,xR,res
+
+inquire(file='problem.data', exist=ex)
+if (ex) then
+ OPEN(UNIT=30, FILE='problem.data', STATUS="OLD")
+ READ(30,NML=ProblemData)                               ! read into variables in namelist
+ CLOSE(30)
+else
+ print *, 'Missing problem.data. Exiting.'
+ stop           ! exit if missing
+end if
+}}}
+
+The state vectors W,,L,, and W,,R,, are created with overloaded constructors that also calculate useful parameters A and B and the speed of sound c for use later in the program.
+
+
+These state vectors are then passed to the Newton-Raphson solver, which calls a subroutine to guess a good value for p,,0,, (currently not very effective - just uses tol) and then iterates until the tolerance for change in pressure between steps is met.
+
+
+{{{#!latex
+$f_reg(p_*) = (p_* - p_reg) [\frac{2}{(\gamma + 1) \rho_reg (p_* + \frac{\gamma - 1}{\gamma + 1} p_reg)}]$
+}}}
+
+
+See below for derivation of function f.
+
+From the shock jump conditions 
+
+{{{#!latex
+$\rho_L \hat{u_L} = \rho_{*L} \hat{u_*}$
+
+$\rho_L \hat{u_L}^2 + p_L = \rho_{*L} \hat{u_*}^2 + p_*$
+
+$\hat{u_L}(\hat{E_L} + p_L) = \hat{u_*}(\hat{E_{*L}} + p_*)$
+}}}
+
+and the energy and EOS of an ideal gas 
+
+{{{#!latex
+$e = \frac{p}{(\gamma - 1) \rho}$
+
+$p = \rho R T$
+}}}
+
+where $\hat{u} = u_L - S_L$ is the velocity of the fluid in the frame of the shock, the function $f_L(p_*) = u_L - u_*$,
+
+{{{#!latex
+$f_L(p_*) = (p_* - p_L) [\frac{2}{(\gamma + 1) \rho_L (p_* + \frac{\gamma - 1}{\gamma + 1} p_L)}]$
+}}}
+
+can be derived to connect the left state with the pressure in the central, interaction region in the case of a shock. Similarly, in the case of a rarefaction, instead of the shock jump conditions, the isentropic equation
+
+{{{#!latex
+$\frac{p_*}{\rho_*} = \frac{p_L}{\rho_L}$
+}}}
+
+and the Generalized Riemann Invariant
+
+{{{#!latex
+$u_L + \frac{2 c_L}{\gamma - 1} = u_* + \frac{2 c_{*L}}{\gamma - 1}$
+}}}
+
+where $c = \sqrt{\frac{\gamma p}{\rho}}$ is the speed of sound in the region, can be used to construct
+
+{{{#!latex
+$f_L(p_*) = \frac{2 c_L}{\gamma - 1}[(\frac{p_*}{p_L})^{\frac{\gamma - 1}{2 \gamma}} - 1]$
+}}}
+
+The same analysis can be performed on the right side of the system to yield the same results.
+
+The Riemann problem can then be solved by finding the value $p = p_*$ such that 
+
+{{{#!latex
+$f_L(p_*) + f_R(p_*) + u_R - u_L = 0$
+}}}
+
+The function and its first and second derivatives are well-behaved and monotonic, so we can use a simple Newton-Raphson iterative method, derived from a Taylor series expansion of f(p_*):
+
+If our initial guess is $p_k$ and $f(p_k) \ne 0$, we want to change $p_k$ so that $f(p_k + \delta) = f(p_{k+1}) = 0$; the Taylor expansion is given by
+
+{{{#!latex
+$f(p_k + \delta) = f(p_k) + \delta f'(p_k) + O(\delta^2)$
+}}}
+
+Dropping second order and higher terms, since we have $f(p_k + \delta) = 0$, we know $f(p_k) + \delta f'(p_k) = 0$, and solving for $\delta$,
+
+{{{#!latex
+$\delta = -\frac{f(p_k)}{f'(p_k)}$
+}}}
+
+and 
+
+{{{#!latex
+$p_{k + 1} = p_k + \delta = p_k - \frac{f(p_k)}{f'(p_k)}$
+}}}
+
+The iteration procedure is terminated at the term such that
+
+{{{#!latex
+$|\frac{2 (p_{k + 1} - p_k)}{p_{k + 1} + p_k}| < Tol$
+}}}
+
+where $Tol$ is the (generally small) tolerance passed to the solver. 
+
+Source: See Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics