Exact Riemann solver for perfect gas mixture: problem with Newton's method convergence

Question

I'm trying to solve multicomponent Euler equations for perfect gas mixture with Godunov-like scheme using exact Riemann solver. Of course, some approximate solver would probably be more cost-effective, but for now I'd like to make it work with exact solver.

Here is the problematic part of algorithm:

We have left and right conditions, for which we calculated relevant values, namely, density, velocity, pressure, heat capacity ratio and sound speed: $\rho_l, u_l, p_l, \gamma_l, a_l$ and $\rho_r, u_r, p_r, \gamma_r, a_r$.

we try to find the velocity and pressure $v, p$ in the middle state that may be reached through single shock wave or rarefaction wave from left and right states at the same time.

We construct functions $\phi_l(p), \phi_r(p)$: $$ \phi_l(p) = \begin{cases} u_l + \frac{2a_l}{\sqrt{2\gamma_l(\gamma_l-1)}}\frac{\left(1-\frac{p}{p_l}\right)}{\sqrt{1+\beta\frac{p}{p_l}}} & \text{if} \;\; p \ge p_l \;\; \text{(shock wave)} \\ u_l + \frac{2a_l}{\gamma_l-1}\left( 1-\left(\frac{p}{p_l}\right)^{\alpha} \right) & \text{if} \;\; p < p_l \;\; \text{(rarefaction)} \end{cases} $$ $$ \phi_r(p) = \begin{cases} u_r - \frac{2a_r}{\sqrt{2\gamma_r(\gamma_r-1)}}\frac{\left(1-\frac{p}{p_r}\right)}{\sqrt{1+\frac{\gamma_r+1}{\gamma_r-1}\frac{p}{p_r}}} & \text{if} \;\; p \ge p_r \;\; \text{(shock wave)} \\ u_r - \frac{2a_r}{\gamma_r-1}\left( 1-\left(\frac{p}{p_r}\right)^{\alpha} \right) & \text{if} \;\; p < p_r \;\; \text{(rarefaction)} \end{cases} $$ where $\alpha = \frac{\gamma-1}{2\gamma}$, $\beta = \frac{\gamma+1}{\gamma-1}$;

and their derivatives:

$$ \phi'_l(p) = \begin{cases} \frac{2a_l}{\sqrt{2\gamma_l(\gamma_l-1)}} \frac{ \frac{-\sqrt{1+\beta\frac{p}{p_l}}}{p_l} - \frac{ \beta \left(1-\frac{p}{p_l}\right)} {2p_l \sqrt{1+\beta\frac{p}{p_l}}}} {1+\beta\frac{p}{p_l}} & \text{if} \;\; p \ge p_l \;\; \text{(shock wave)} \\ - \alpha \frac{1}{p_l^{\alpha}} \frac{2a_l}{\gamma_l-1}\left(\frac{p}{p_l}\right)^{\alpha - 1} & \text{if} \;\; p < p_l \;\; \text{(rarefaction)} \end{cases} $$ $$ \phi'_r(p) = \begin{cases} -\frac{2a_r}{\sqrt{2\gamma_r(\gamma_r-1)}} \frac{ \frac{-\sqrt{1+\beta\frac{p}{p_r}}}{p_r} - \frac{ \beta \left(1-\frac{p}{p_r}\right)} {2p_r \sqrt{1+\beta\frac{p}{p_r}}}} {1+\beta\frac{p}{p_r}} & \text{if} \;\; p \ge p_r \;\; \text{(shock wave)} \\ \alpha \frac{1}{p_r^{\alpha}} \frac{2a_r}{\gamma_r-1}\left(\frac{p}{p_r}\right)^{\alpha - 1} & \text{if} \;\; p < p_r \;\; \text{(rarefaction)} \end{cases} $$

and try to find $p$ so that $$ \phi(p) = \phi_l(p) - \phi_r(p) = 0 $$

For that we use Newton's method: $$ p_{n+1} = p_{n} - \frac{\phi(p_n)}{\phi'(p_n)} $$ starting with $p_0 = (p_l+p_r)/2$ until it converges.

For single-component Euler equations it almost always converges extremely fast: typical Sod's problem $(\rho_l,u_l,p_l = 1,0,1$, $\rho_r,u_r,p_r = 0.125,0,0.1)$ converges in 5-6 steps. But for my multi-component problems typical values are something like $p \approx 10^5-10^7 Pa$, while $u,a \approx 10^2-10^3 m/s$ and also $\phi(p), \phi'(p) \approx 10^2-10^3$. Thus the pressure changes very little at every step of the Newton's algorithm (for example only by $\Delta p \approx 0.3 Pa$ while $p_0 = 550000 Pa$) and the convergence is extremely slow.

In some cases, I was able to achieve convergence within hundred of steps by artificially multiplying $\Delta p$ by something like $10^4$ and it gave the correct solution. But I wasn't able to construct stable algorithm for dynamic multiplier for $\Delta p$ that would give at least reasonable convergence speed in general case.

Is there some more general way to speed up the convergence here? Or may some different way to determine middle state pressure?

Answer 1

I have found very simple solution! At least, it seem to work well with my test case. Instead of explicit calculation of $\phi'(p_n)$, I use it's simple approximation: $$ \phi'(p_n) \approx \frac{\phi(p_n)-\phi(p_{n-1})}{p_n-p_{n-1}} $$ To use this algorithm, we need to have two starting values ($p_0, p_1$) instead of one ($p_0$). Currently I calculate $p_1$ via single step of original Newton's method with explicit $\phi'(p_0)$. As stated above, it gives very small difference between $p_1$ and $p_0$, but following steps with approximate $\phi'(p_n)$ produce very good increments and the method converges in usual 4-5 steps.

Update: As Praveen Chandrashekar mentioned below, I should have checked my derivatives again. There is a pretty simple error in derivative for rarefaction wave functions, it should contain $\frac{1}{p_l}$ instead of $\frac{1}{p_l^\alpha}$ (and similar thing for $p_r$). With correct derivatives, original Newton's method converges well.